From 8116658e0d2993ea8750d2bed7e3dc7f4df4a09e Mon Sep 17 00:00:00 2001 From: "Timur A. Fatkhullin" Date: Fri, 27 Mar 2026 17:29:46 +0300 Subject: [PATCH] ... --- include/mcc/mcc_pcm_construct.h | 80 +++++++++++++++++++++------------ 1 file changed, 52 insertions(+), 28 deletions(-) diff --git a/include/mcc/mcc_pcm_construct.h b/include/mcc/mcc_pcm_construct.h index d65beb3..12037db 100644 --- a/include/mcc/mcc_pcm_construct.h +++ b/include/mcc/mcc_pcm_construct.h @@ -320,6 +320,7 @@ protected: { compute_result_t result{.pcm_type = pcm_data.type, .error = MccDefaultPCMConstructorErrorCode::ERROR_OK}; + // at least border knots must be given ([t_min, t_max]) if (pcm_data.bspline.knotsX.size() < 2 || pcm_data.bspline.knotsY.size() < 2) { result.error = MccDefaultPCMConstructorErrorCode::ERROR_INVALID_KNOTS_NUMBER; @@ -328,31 +329,53 @@ protected: double resi2x, resi2y; // sum of fitting residuals squares + // full knot vectors: [t_min, t_min, t_min, t_min, ..., t_1, t_2, ..., t_max, t_max, t_max, t_max] std::vector tx(pcm_data.bspline.knotsX.size() + 6), ty(pcm_data.bspline.knotsY.size() + 6); + std::vector theta_hw(numberOfPoints()), theta_tag(numberOfPoints()); + + /* + WARNING: + FITPACK B-spline on sphere: in the fitting routines the first angle argument is THETA - co-latitude + coordinate and the second one is PHI - co-longitude!!! + + In the fitting routines the first angle argument THETA must be in the range of [0.0, PI] + (while celestial declination is in the range of [-PI/2, PI/2])!!! + */ + + for (size_t i = 0; i < numberOfPoints(); ++i) { + theta_hw[i] = _table.hw_colat[i] + MCC_HALF_PI; + } + + tx[0] = tx[1] = tx[2] = pcm_data.bspline.knotsX[0]; + tx[tx.size() - 1] = tx[tx.size() - 2] = tx[tx.size() - 3] = + pcm_data.bspline.knotsX[pcm_data.bspline.knotsX.size() - 1]; + + ty[0] = ty[1] = ty[2] = pcm_data.bspline.knotsY[0] + MCC_HALF_PI; + ty[ty.size() - 1] = ty[ty.size() - 2] = ty[ty.size() - 3] = + pcm_data.bspline.knotsY[pcm_data.bspline.knotsY.size() - 1] + MCC_HALF_PI; + + for (size_t i = 0; i < pcm_data.bspline.knotsX.size(); ++i) { + tx[3 + i] = pcm_data.bspline.knotsX[i]; + } + for (size_t i = 0; i < pcm_data.bspline.knotsY.size(); ++i) { + ty[3 + i] = pcm_data.bspline.knotsY[i] + MCC_HALF_PI; + } size_t Ncoeffs = (tx.size() - 4) * (ty.size() - 4); pcm_data.bspline.coeffsX.resize(Ncoeffs); pcm_data.bspline.coeffsY.resize(Ncoeffs); - /* - WARNING: - FITPACK B-spline on sphere: in the fitting routines the first angle argument is THETA - co-latitude - coordinate and the second one is PHI - co-longitude - - */ // direct (celestial = encoder + pcm) - result.bspline_fit_err = bsplines::fitpack_sphere_fit(_table.hw_colat, _table.hw_colon, _table.colon_res, 1.0, - pcm_data.bspline.knotsY, pcm_data.bspline.knotsX, + result.bspline_fit_err = bsplines::fitpack_sphere_fit(theta_hw, _table.hw_colon, _table.colon_res, 1.0, ty, tx, pcm_data.bspline.coeffsX, resi2x); if (result.bspline_fit_err > 0) { result.error = MccDefaultPCMConstructorErrorCode::ERROR_BSPLINE_FIT; return result; } - result.bspline_fit_err = bsplines::fitpack_sphere_fit(_table.hw_colat, _table.hw_colon, _table.colat_res, 1.0, - pcm_data.bspline.knotsY, pcm_data.bspline.knotsX, + result.bspline_fit_err = bsplines::fitpack_sphere_fit(theta_hw, _table.hw_colon, _table.colon_res, 1.0, ty, tx, pcm_data.bspline.coeffsY, resi2y); if (result.bspline_fit_err > 0) { result.error = MccDefaultPCMConstructorErrorCode::ERROR_BSPLINE_FIT; @@ -360,11 +383,11 @@ protected: } - bsplines::fitpack_eval_spl2d(pcm_data.bspline.knotsY, pcm_data.bspline.knotsX, pcm_data.bspline.coeffsX, - _table.hw_colat, _table.hw_colon, result.model_colon); // get fitted residuals!!! + bsplines::fitpack_eval_spl2d(ty, tx, pcm_data.bspline.coeffsX, theta_hw, _table.hw_colon, + result.model_colon); // get fitted residuals!!! - bsplines::fitpack_eval_spl2d(pcm_data.bspline.knotsY, pcm_data.bspline.knotsX, pcm_data.bspline.coeffsY, - _table.hw_colat, _table.hw_colon, result.model_colat); // get fitted residuals!!! + bsplines::fitpack_eval_spl2d(ty, tx, pcm_data.bspline.coeffsY, theta_hw, _table.hw_colon, + result.model_colat); // get fitted residuals!!! result.colon_res.resize(numberOfPoints()); @@ -372,8 +395,8 @@ protected: for (size_t i = 0; i < numberOfPoints(); ++i) { result.colon_res[i] = _table.colon_res[i] - result.model_colon[i]; // = target - model result.colat_res[i] = _table.colat_res[i] - result.model_colat[i]; // = target - model - result.model_colon[i] += _table.hw_colon[i]; // == hw + fitted_pcmX - result.model_colat[i] += _table.hw_colat[i]; // == hw + fitted_pcmY + // result.model_colon[i] += _table.hw_colon[i]; // == hw + fitted_pcmX + // result.model_colat[i] += _table.hw_colat[i]; // == hw + fitted_pcmY } if (pcm_data.type == MccDefaultPCMType::PCM_TYPE_BSPLINE) { @@ -381,6 +404,8 @@ protected: pcm_data.inverseBspline.coeffsX.resize(Ncoeffs); pcm_data.inverseBspline.coeffsY.resize(Ncoeffs); + std::vector theta_tag(numberOfPoints()); + // inverse (encoder = celestial + pcm) std::vector colon_res = _table.colon_res; @@ -388,30 +413,28 @@ protected: for (size_t i = 0; i < colat_res.size(); ++i) { colon_res[i] = -colon_res[i]; colat_res[i] = -colat_res[i]; + + theta_tag[i] = _table.target_colat[i] + MCC_HALF_PI; } - result.bspline_fit_err = bsplines::fitpack_sphere_fit(_table.target_colon, _table.target_colat, colon_res, - 1.0, pcm_data.bspline.knotsY, pcm_data.bspline.knotsX, - pcm_data.bspline.inverseCoeffsX, resi2x); + result.bspline_fit_err = bsplines::fitpack_sphere_fit(theta_tag, _table.target_colon, colon_res, 1.0, ty, + tx, pcm_data.bspline.inverseCoeffsX, resi2x); if (result.bspline_fit_err > 0) { result.error = MccDefaultPCMConstructorErrorCode::ERROR_BSPLINE_FIT; return result; } - result.bspline_fit_err = bsplines::fitpack_sphere_fit(_table.target_colon, _table.target_colat, colat_res, - 1.0, pcm_data.bspline.knotsY, pcm_data.bspline.knotsX, - pcm_data.bspline.inverseCoeffsY, resi2y); + result.bspline_fit_err = bsplines::fitpack_sphere_fit(theta_tag, _table.target_colon, colat_res, 1.0, ty, + tx, pcm_data.bspline.inverseCoeffsY, resi2y); if (result.bspline_fit_err > 0) { result.error = MccDefaultPCMConstructorErrorCode::ERROR_BSPLINE_FIT; return result; } - bsplines::fitpack_eval_spl2d(pcm_data.bspline.knotsY, pcm_data.bspline.knotsX, - pcm_data.bspline.inverseCoeffsX, _table.hw_colat, _table.hw_colon, + bsplines::fitpack_eval_spl2d(ty, tx, pcm_data.bspline.inverseCoeffsX, theta_tag, _table.target_colon, result.inv_model_colon); // get fitted residuals!!! - bsplines::fitpack_eval_spl2d(pcm_data.bspline.knotsY, pcm_data.bspline.knotsX, - pcm_data.bspline.inverseCoeffsY, _table.hw_colat, _table.hw_colon, + bsplines::fitpack_eval_spl2d(ty, tx, pcm_data.bspline.inverseCoeffsY, theta_tag, _table.target_colon, result.inv_model_colat); // get fitted residuals!!! result.inv_colon_res.resize(numberOfPoints()); @@ -419,8 +442,9 @@ protected: for (size_t i = 0; i < numberOfPoints(); ++i) { result.inv_colon_res[i] = _table.colon_res[i] - result.inv_model_colon[i]; // = hw - model result.inv_colat_res[i] = _table.colat_res[i] - result.inv_model_colat[i]; // = hw - model - result.inv_model_colon[i] += _table.target_colon[i]; // == target + fitted_pcmX - result.inv_model_colat[i] += _table.target_colat[i]; // == target + fitted_pcmY + // result.inv_model_colon[i] += _table.target_colon[i]; // == target + + // fitted_pcmX result.inv_model_colat[i] += _table.target_colat[i]; // == target + // + fitted_pcmY } }