From a4b3a32d94dee547b1c83e3e9557e9648cf4ab6a Mon Sep 17 00:00:00 2001 From: "Timur A. Fatkhullin" Date: Sun, 29 Mar 2026 22:18:46 +0300 Subject: [PATCH] ... --- include/mcc/mcc_bsplines.h | 213 +++++++++++++++++++++++++++++--- include/mcc/mcc_pcm_construct.h | 2 +- 2 files changed, 194 insertions(+), 21 deletions(-) diff --git a/include/mcc/mcc_bsplines.h b/include/mcc/mcc_bsplines.h index 0b50dd5..5781628 100644 --- a/include/mcc/mcc_bsplines.h +++ b/include/mcc/mcc_bsplines.h @@ -7,6 +7,8 @@ #include extern "C" { + +/* fitting on XY-grid */ void surfit(int* iopt, int* m, double* x, @@ -38,6 +40,7 @@ void surfit(int* iopt, int* kwrk, int* ier); +/* XY-grid */ void bispev(double* tx, int* nx, double* ty, @@ -56,6 +59,7 @@ void bispev(double* tx, int* kwrk, int* ier); +/* XY-grid */ void parder(double* tx, int* nx, double* ty, @@ -76,7 +80,7 @@ void parder(double* tx, int* kwrk, int* ier); - +/* fitting on sphere */ void sphere(int* iopt, int* m, double* teta, @@ -100,8 +104,44 @@ void sphere(int* iopt, int* iwrk, int* kwrk, int* ier); -} +/* calculation for set of points (not grid) */ +void bispeu(double* tx, + int* nx, + double* ty, + int* ny, + double* c, + int* kx, + int* ky, + double* x, + double* y, + double* z, + int* m, + double* wrk, + int* lwrk, + int* ier); + + +/* calculation for set of points (not grid) */ +void pardeu(double* tx, + int* nx, + double* ty, + int* ny, + double* c, + int* kx, + int* ky, + int* nux, + int* nuy, + double* x, + double* y, + double* z, + int* m, + double* wrk, + int* lwrk, + int* iwrk, + int* kwrk, + int* ier); +} namespace mcc::bsplines @@ -275,14 +315,14 @@ template -int fitpack_eval_spl2d(const TXT& tx, - const TYT& ty, - const CoeffT& coeffs, - const XT& x, - const YT& y, - FuncT& func, - int kx = 3, - int ky = 3) +int fitpack_eval_spl2d_grid(const TXT& tx, + const TYT& ty, + const CoeffT& coeffs, + const XT& x, + const YT& y, + FuncT& func, + int kx = 3, + int ky = 3) { static_assert(std::same_as, double> && std::same_as, double> && @@ -334,6 +374,70 @@ int fitpack_eval_spl2d(const TXT& tx, } +/* for set of points */ +template +int fitpack_eval_spl2d(const TXT& tx, + const TYT& ty, + const CoeffT& coeffs, + const XT& x, + const YT& y, + FuncT& func, + int kx = 3, + int ky = 3) +{ + static_assert(std::same_as, double> && + std::same_as, double> && + std::same_as, double> && + std::same_as, double> && + std::same_as, double> && + std::same_as, double>, + "Input ranges elements type must be double!"); + + if (kx < 0 || ky < 0) { + return 10; + } + + if (std::ranges::size(x) != std::ranges::size(y)) { + return 10; + } + + int m = std::ranges::size(x); + + int ier = 0; + + int ntx = std::ranges::size(tx); + int nty = std::ranges::size(ty); + + auto n_coeffs = (ntx - kx - 1) * (nty - ky - 1); + if (std::ranges::size(coeffs) < n_coeffs) { + return 10; + } + + if (std::ranges::size(func) < m) { + std::ranges::fill_n(std::back_inserter(func), m - std::ranges::size(func), 0.0); + } + + auto tx_ptr = const_cast(std::ranges::data(tx)); + auto ty_ptr = const_cast(std::ranges::data(ty)); + auto coeffs_ptr = const_cast(std::ranges::data(coeffs)); + auto x_ptr = const_cast(std::ranges::data(x)); + auto y_ptr = const_cast(std::ranges::data(y)); + + int lwrk = kx + ky + 2; + std::vector wrk(lwrk); + + bispeu(tx_ptr, &ntx, ty_ptr, &nty, coeffs_ptr, &kx, &ky, x_ptr, y_ptr, std::ranges::data(func), &m, wrk.data(), + &lwrk, &ier); + + return ier; +} + + // scalar x,y version template -int fitpack_parder_spl2d(const TXT& tx, - const TYT& ty, - const CoeffT& coeffs, - const XT& x, - const YT& y, - PderT& pder, - int dx, // partial derivatives order along X - int dy, // partial derivatives order along Y - int kx = 3, - int ky = 3) +int fitpack_parder_spl2d_grid(const TXT& tx, + const TYT& ty, + const CoeffT& coeffs, + const XT& x, + const YT& y, + PderT& pder, + int dx, // partial derivatives order along X + int dy, // partial derivatives order along Y + int kx = 3, + int ky = 3) { static_assert(std::same_as, double> && std::same_as, double> && @@ -438,6 +542,75 @@ int fitpack_parder_spl2d(const TXT& tx, return ier; } + +template +int fitpack_parder_spl2d(const TXT& tx, + const TYT& ty, + const CoeffT& coeffs, + const XT& x, + const YT& y, + PderT& pder, + int dx, // partial derivatives order along X + int dy, // partial derivatives order along Y + int kx = 3, + int ky = 3) +{ + static_assert(std::same_as, double> && + std::same_as, double> && + std::same_as, double> && + std::same_as, double> && + std::same_as, double> && + std::same_as, double>, + "Input ranges elements type must be double!"); + + if (kx < 0 || ky < 0 || dx < 0 || dy < 0) { + return 10; + } + + int ntx = std::ranges::size(tx); + int nty = std::ranges::size(ty); + + auto n_coeffs = (ntx - kx - 1) * (nty - ky - 1); + if (std::ranges::size(coeffs) < n_coeffs) { + return 10; + } + + int m = std::ranges::size(x); + if (std::ranges::size(x) != std::ranges::size(y)) { + return 10; + } + + if (std::ranges::size(pder) < m) { + std::ranges::fill_n(std::back_inserter(pder), m - std::ranges::size(pder), 0.0); + } + + + auto tx_ptr = const_cast(std::ranges::data(tx)); + auto ty_ptr = const_cast(std::ranges::data(ty)); + auto coeffs_ptr = const_cast(std::ranges::data(coeffs)); + auto x_ptr = const_cast(std::ranges::data(x)); + auto y_ptr = const_cast(std::ranges::data(y)); + + int ier = 0; + + int lwrk = m * (kx + 1 - dx) + m * (ky + 1 - dy) + (ntx - kx - 1) * (nty - ky - 1); + std::vector wrk(lwrk); + + int kwrk = 2 * m; + std::vector iwrk(kwrk); + + pardeu(tx_ptr, &ntx, ty_ptr, &nty, coeffs_ptr, &kx, &ky, &dx, &dy, x_ptr, y_ptr, std::ranges::data(pder), &m, + wrk.data(), &lwrk, iwrk.data(), &kwrk, &ier); + + return ier; +} + + // scalar version template 0) { result.error = MccDefaultPCMConstructorErrorCode::ERROR_BSPLINE_FIT;