From 9ad5bfb09b810ca74c98818ac984b8a737ae8c6a Mon Sep 17 00:00:00 2001 From: "Timur A. Fatkhullin" Date: Wed, 18 Mar 2026 23:34:46 +0300 Subject: [PATCH] ... --- include/mcc/mcc_pcm.h | 33 ++-- include/mcc/mcc_pcm_construct.h | 295 +++++++++++++++++++++++++++++--- 2 files changed, 287 insertions(+), 41 deletions(-) diff --git a/include/mcc/mcc_pcm.h b/include/mcc/mcc_pcm.h index 2dcd2c7..8ee6ef1 100644 --- a/include/mcc/mcc_pcm.h +++ b/include/mcc/mcc_pcm.h @@ -166,6 +166,9 @@ public: std::vector coeffsX{}; std::vector coeffsY{}; + + std::vector inverseCoeffsX{}; + std::vector inverseCoeffsY{}; }; #endif @@ -466,7 +469,7 @@ private: #ifdef USE_BSPLINE_PCM pcm_bspline_t* bspline = &_pcmData.bspline; - pcm_bspline_t* inv_bspline = &_pcmData.inverseBspline; + // pcm_bspline_t* inv_bspline = &_pcmData.inverseBspline; #endif #ifdef USE_BSPLINE_PCM @@ -567,9 +570,8 @@ private: int ret = 0; if constexpr (!std::is_null_pointer_v) { - ret = - bsplines::fitpack_parder_spl2d(inv_bspline->knotsX, inv_bspline->knotsY, inv_bspline->coeffsX, x, y, - dspl_valX, 1, 0, inv_bspline->bsplDegreeX, inv_bspline->bsplDegreeY); + ret = bsplines::fitpack_parder_spl2d(bspline->knotsX, bspline->knotsY, bspline->coeffsX, x, y, + dspl_valX, 1, 0, bspline->bsplDegreeX, bspline->bsplDegreeY); if (ret) { return MccDefaultPCMErrorCode::ERROR_INVALID_INPUTS_BISPLEV; @@ -577,9 +579,8 @@ private: derivX->pcmX += dspl_valX; // dpcmX/dX - ret = - bsplines::fitpack_parder_spl2d(inv_bspline->knotsX, inv_bspline->knotsY, inv_bspline->coeffsX, x, y, - dspl_valX, 0, 1, inv_bspline->bsplDegreeX, inv_bspline->bsplDegreeY); + ret = bsplines::fitpack_parder_spl2d(bspline->knotsX, bspline->knotsY, bspline->coeffsX, x, y, + dspl_valX, 0, 1, bspline->bsplDegreeX, bspline->bsplDegreeY); if (ret) { return MccDefaultPCMErrorCode::ERROR_INVALID_INPUTS_BISPLEV; @@ -589,9 +590,8 @@ private: } if constexpr (!std::is_null_pointer_v) { - ret = - bsplines::fitpack_parder_spl2d(inv_bspline->knotsX, inv_bspline->knotsY, inv_bspline->coeffsY, x, y, - dspl_valY, 1, 0, inv_bspline->bsplDegreeX, inv_bspline->bsplDegreeY); + ret = bsplines::fitpack_parder_spl2d(bspline->knotsX, bspline->knotsY, bspline->coeffsY, x, y, + dspl_valY, 1, 0, bspline->bsplDegreeX, bspline->bsplDegreeY); if (ret) { return MccDefaultPCMErrorCode::ERROR_INVALID_INPUTS_BISPLEV; @@ -599,9 +599,8 @@ private: derivY->pcmX += dspl_valY; // dpcmY/dX - ret = - bsplines::fitpack_parder_spl2d(inv_bspline->knotsX, inv_bspline->knotsY, inv_bspline->coeffsY, x, y, - dspl_valY, 0, 1, inv_bspline->bsplDegreeX, inv_bspline->bsplDegreeY); + ret = bsplines::fitpack_parder_spl2d(bspline->knotsX, bspline->knotsY, bspline->coeffsY, x, y, + dspl_valY, 0, 1, bspline->bsplDegreeX, bspline->bsplDegreeY); if (ret) { return MccDefaultPCMErrorCode::ERROR_INVALID_INPUTS_BISPLEV; @@ -615,8 +614,8 @@ private: if (_pcmData.type == MccDefaultPCMType::PCM_TYPE_BSPLINE && inverse) { double spl_valX, spl_valY; - int ret = bsplines::fitpack_eval_spl2d(inv_bspline->knotsX, inv_bspline->knotsY, inv_bspline->coeffsX, x, y, - spl_valX, inv_bspline->bsplDegreeX, inv_bspline->bsplDegreeY); + int ret = bsplines::fitpack_eval_spl2d(bspline->knotsX, bspline->knotsY, bspline->inverseCoeffsX, x, y, + spl_valX, bspline->bsplDegreeX, bspline->bsplDegreeY); if (ret) { res->pcmX = std::numeric_limits::quiet_NaN(); @@ -626,8 +625,8 @@ private: } - ret = bsplines::fitpack_eval_spl2d(inv_bspline->knotsX, inv_bspline->knotsY, inv_bspline->coeffsY, x, y, - spl_valY, inv_bspline->bsplDegreeX, inv_bspline->bsplDegreeY); + ret = bsplines::fitpack_eval_spl2d(bspline->knotsX, bspline->knotsY, bspline->inverseCoeffsY, x, y, + spl_valY, bspline->bsplDegreeX, bspline->bsplDegreeY); if (ret) { res->pcmX = std::numeric_limits::quiet_NaN(); diff --git a/include/mcc/mcc_pcm_construct.h b/include/mcc/mcc_pcm_construct.h index e9cec93..1f02371 100644 --- a/include/mcc/mcc_pcm_construct.h +++ b/include/mcc/mcc_pcm_construct.h @@ -18,8 +18,78 @@ namespace mcc::impl { +enum class MccDefaultPCMConstructorErrorCode : int { + ERROR_OK, + ERROR_NOT_ENOUGH_DATA, +#ifdef USE_BSPLINE_PCM + ERROR_INVALID_KNOTS_NUMBER, + ERROR_BSPLINE_FIT +#endif +}; + +// error category + +struct MccDefaultPCMConstructorErrorCategory : std::error_category { + MccDefaultPCMConstructorErrorCategory() = default; + + const char* name() const noexcept + { + return "MCC-DEFAULT-PCM-CONSTRUCTOR-ERROR-CATEGORY"; + } + + std::string message(int ec) const + { + MccDefaultPCMConstructorErrorCode err = static_cast(ec); + + switch (err) { + case MccDefaultPCMConstructorErrorCode::ERROR_OK: + return "OK"; + case MccDefaultPCMConstructorErrorCode::ERROR_NOT_ENOUGH_DATA: + return "not enough data point"; + +#ifdef USE_BSPLINE_PCM + case MccDefaultPCMConstructorErrorCode::ERROR_INVALID_KNOTS_NUMBER: + return "invalid number of B-spline knots"; + case MccDefaultPCMConstructorErrorCode::ERROR_BSPLINE_FIT: + return "B-spline fitting error"; +#endif + default: + return "UNKNOWN"; + } + } + + static const MccDefaultPCMConstructorErrorCategory& get() + { + static const MccDefaultPCMConstructorErrorCategory constInst; + return constInst; + } +}; + + +inline std::error_code make_error_code(MccDefaultPCMConstructorErrorCode ec) +{ + return std::error_code(static_cast(ec), MccDefaultPCMConstructorErrorCategory::get()); +} + +} // namespace mcc::impl + + +namespace std +{ + +template <> +class is_error_code_enum : public true_type +{ +}; + + +} // namespace std + +namespace mcc::impl +{ + template -class MccPcmRefPointTable +class MccPCMConstructor { public: using ref_coordpair_t = std::conditional_t, @@ -34,6 +104,29 @@ public: double colon_res, colat_res; // target - hw }; + + struct table_t { + std::vector target_colon, target_colat; + std::vector hw_colon, hw_colat; + std::vector colon_res, colat_res; // target - hw + }; + + struct compute_result_t { + MccDefaultPCMType pcm_type; + MccError error{}; // final model computation error + + std::vector model_colon, model_colat; // fitted model values + std::vector colon_res, colat_res; // target - model + +#ifdef USE_BSPLINE_PCM + int bspline_fit_err{}; // bivariate B-spline fitting exit code (see FITPACK) + + // quantities below are computed only fo pcm_type == MccDefaultPCMType::PCM_TYPE_BSPLINE + std::vector inv_model_colon, inv_model_colat; // fitted inverse model values + std::vector inv_colon_res, inv_colat_res; // encoder - model +#endif + }; + MccError addPoint(mcc_skypoint_c auto target_coords, mcc_coord_pair_c auto const& hw_counts, mcc_coord_epoch_c auto const& ref_epoch) @@ -44,8 +137,17 @@ public: return mcc_deduced_err(err, MccDefaultPCMErrorCode::ERROR_CCTE); } - _table.push_back({target_coords.co_lon(), target_coords.co_lat(), hw_counts.x(), hw_counts.y(), - target_coords.co_lon() - hw_counts.x(), target_coords.co_lat() - hw_counts.y()}); + // _table.push_back({target_coords.co_lon(), target_coords.co_lat(), hw_counts.x(), hw_counts.y(), + // target_coords.co_lon() - hw_counts.x(), target_coords.co_lat() - hw_counts.y()}); + + _table.target_colon.emplace_back(target_coords.co_lon()); + _table.target_colat.emplace_back(target_coords.co_lon()); + + _table.hw_colon.emplace_back(hw_counts.x()); + _table.hw_colat.emplace_back(hw_counts.y()); + + _table.colon_res.emplace_back(target_coords.co_lon() - hw_counts.x()); + _table.colat_res.emplace_back(target_coords.co_lat() - hw_counts.y()); return {}; } @@ -59,6 +161,24 @@ public: } + size_t numberOfPoints() const + { + return _table.colon_res.size(); + } + + void deletePoints() + { + _table.target_colon.clear(); + _table.target_colat.clear(); + + _table.hw_colon.clear(); + _table.hw_colat.clear(); + + _table.colon_res.clear(); + _table.colat_res.clear(); + } + + // // for B-splines interior knots along axes must be given in 'pcm_data' // NOTE: the size of the interior knots array must be at least 2 as @@ -68,11 +188,17 @@ public: // // WARNING: the input knots for inverse B-spline are ignored so the direct and inverse B-spline coefficients are // calculated on the same mesh! - MccError computeModel(MccDefaultPCM::pcm_data_t& pcm_data) + compute_result_t computeModel(MccDefaultPCM::pcm_data_t& pcm_data) { + compute_result_t result{.pcm_type = pcm_data.type, .error = MccDefaultPCMConstructorErrorCode::ERROR_OK}; + size_t min_data_size = 2; // 2 is for BSPLINE - if (pcm_data.type == MccDefaultPCMType::PCM_TYPE_GEOMETRY) { + if (pcm_data.type == MccDefaultPCMType::PCM_TYPE_GEOMETRY +#ifdef USE_BSPLINE_PCM + || pcm_data.type == MccDefaultPCMType::PCM_TYPE_GEOMETRY_BSPLINE +#endif + ) { if constexpr (MOUNT_TYPE == MccMountType::FORK_TYPE) { min_data_size = 9; } else { @@ -80,15 +206,112 @@ public: } if (_table.size() < min_data_size) { - // return error + result.error = MccDefaultPCMConstructorErrorCode::ERROR_NOT_ENOUGH_DATA; + return result; } + + // robust linear regression with Tukey's loss function } #ifdef USE_BSPLINE_PCM - if (pcm_data.bspline.knotsX.size() < 2 || pcm_data.bspline.knotsY.size() < 2) { - // TODO: return error + if (pcm_data.type == MccDefaultPCMType::PCM_TYPE_BSPLINE || + pcm_data.type == MccDefaultPCMType::PCM_TYPE_GEOMETRY_BSPLINE) { + if (pcm_data.bspline.knotsX.size() < 2 || pcm_data.bspline.knotsY.size() < 2) { + return MccDefaultPCMConstructorErrorCode::ERROR_INVALID_KNOTS_NUMBER; + } + + double resi2x, resi2y; // fitting residuals + + std::vector tx(pcm_data.bspline.knotsX.size() + 6), ty(pcm_data.bspline.knotsY.size() + 6); + + size_t Ncoeffs = (tx.size() - 4) * (ty.size() - 4); + + pcm_data.bspline.coeffsX.resize(Ncoeffs); + pcm_data.bspline.coeffsY.resize(Ncoeffs); + + if (pcm_data.type == MccDefaultPCMType::PCM_TYPE_BSPLINE) { + if (_table.size() < min_data_size) { + result.error = MccDefaultPCMConstructorErrorCode::ERROR_NOT_ENOUGH_DATA; + return result; + } + + // here both direct and inverse coefficients will be calculated + pcm_data.inverseBspline.coeffsX.resize(Ncoeffs); + pcm_data.inverseBspline.coeffsY.resize(Ncoeffs); + + // 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, 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, pcm_data.bspline.coeffsY, resi2y); + if (result.bspline_fit_err > 0) { + result.error = MccDefaultPCMConstructorErrorCode::ERROR_BSPLINE_FIT; + return result; + } + + + // inverse (encoder = celestial + pcm) + std::vector colon_res = _table.colon_res; + std::vector colat_res = _table.colat_res; + for (size_t i = 0; i < colat_res.size(); ++i) { + colon_res[i] = -colon_res[i]; + colat_res[i] = -colat_res[i]; + } + + 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); + 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); + if (result.bspline_fit_err > 0) { + result.error = MccDefaultPCMConstructorErrorCode::ERROR_BSPLINE_FIT; + return result; + } + } else { // geometry + B-spline + // the fitting for geometrical coefficients is already done above so + // one must fit residuals by bivariate B-splines + + std::vector xres(_table.size()), yres(_table.size()); + // for (size_t i = 0; i < _table.size(); ++i) { + // xres = _table[i].target_colon; + // yres = _table[i].target_colat; + // } + for (size_t i = 0; i < _table.size(); ++i) { + xres = _table.target_colon[i]; + yres = _table.target_colat[i]; + } + } } + } +#endif + +protected: + // std::vector _table; + table_t _table; + + compute_result_t bsplineFitting(MccDefaultPCM::pcm_data_t& pcm_data) + { + compute_result_t result{.pcm_type = pcm_data.type, .error = MccDefaultPCMConstructorErrorCode::ERROR_OK}; + + if (pcm_data.bspline.knotsX.size() < 2 || pcm_data.bspline.knotsY.size() < 2) { + return MccDefaultPCMConstructorErrorCode::ERROR_INVALID_KNOTS_NUMBER; + } + + double resi2x, resi2y; // fitting residuals std::vector tx(pcm_data.bspline.knotsX.size() + 6), ty(pcm_data.bspline.knotsY.size() + 6); @@ -98,29 +321,53 @@ public: pcm_data.bspline.coeffsY.resize(Ncoeffs); if (pcm_data.type == MccDefaultPCMType::PCM_TYPE_BSPLINE) { - if (_table.size() < min_data_size) { - // return error - } - + // here both direct and inverse coefficients will be calculated pcm_data.inverseBspline.coeffsX.resize(Ncoeffs); pcm_data.inverseBspline.coeffsY.resize(Ncoeffs); - } else if (pcm_data.type == MccDefaultPCMType::PCM_TYPE_GEOMETRY_BSPLINE) { - // the fitting for geometrical coefficients is already done above so - // one must fit residuals by bivariate B-splines + // 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, + pcm_data.bspline.coeffsX, resi2x); + if (result.bspline_fit_err > 0) { + result.error = MccDefaultPCMConstructorErrorCode::ERROR_BSPLINE_FIT; + return result; + } - std::vector xres(_table.size()), yres(_table.size()); - for (size_t i = 0; i < _table.size(); ++i) { - xres = _table[i].target_colon; - yres = _table[i].target_colat; + 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, + pcm_data.bspline.coeffsY, resi2y); + if (result.bspline_fit_err > 0) { + result.error = MccDefaultPCMConstructorErrorCode::ERROR_BSPLINE_FIT; + return result; + } + + + // inverse (encoder = celestial + pcm) + std::vector colon_res = _table.colon_res; + std::vector colat_res = _table.colat_res; + for (size_t i = 0; i < colat_res.size(); ++i) { + colon_res[i] = -colon_res[i]; + colat_res[i] = -colat_res[i]; + } + + 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); + 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); + if (result.bspline_fit_err > 0) { + result.error = MccDefaultPCMConstructorErrorCode::ERROR_BSPLINE_FIT; + return result; } } } -#endif - -protected: - std::vector _table; }; - } // namespace mcc::impl \ No newline at end of file