diff --git a/include/mcc/mcc_pcm_construct.h b/include/mcc/mcc_pcm_construct.h index 1f02371..645df3d 100644 --- a/include/mcc/mcc_pcm_construct.h +++ b/include/mcc/mcc_pcm_construct.h @@ -21,6 +21,7 @@ namespace mcc::impl enum class MccDefaultPCMConstructorErrorCode : int { ERROR_OK, ERROR_NOT_ENOUGH_DATA, + ERROR_INVALID_INDEX, #ifdef USE_BSPLINE_PCM ERROR_INVALID_KNOTS_NUMBER, ERROR_BSPLINE_FIT @@ -46,7 +47,8 @@ struct MccDefaultPCMConstructorErrorCategory : std::error_category { return "OK"; case MccDefaultPCMConstructorErrorCode::ERROR_NOT_ENOUGH_DATA: return "not enough data point"; - + case MccDefaultPCMConstructorErrorCode::ERROR_INVALID_INDEX: + return "invalid index of data point"; #ifdef USE_BSPLINE_PCM case MccDefaultPCMConstructorErrorCode::ERROR_INVALID_KNOTS_NUMBER: return "invalid number of B-spline knots"; @@ -106,16 +108,16 @@ public: struct table_t { - std::vector target_colon, target_colat; - std::vector hw_colon, hw_colat; - std::vector colon_res, colat_res; // target - hw + 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 model_colon, model_colat; // fitting model values std::vector colon_res, colat_res; // target - model #ifdef USE_BSPLINE_PCM @@ -127,12 +129,10 @@ public: #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) + MccError addPoint(mcc_skypoint_c auto target_coords, mcc_coord_pair_c auto const& hw_counts) requires(decltype(hw_counts)::pairKind == MccCoordPairKind::COORDS_KIND_XY) { - auto err = target_coords.to(ref_coordpair_t::pairKind, ref_epoch); + auto err = target_coords.to(ref_coordpair_t::pairKind, hw_counts.epoch()); if (err) { return mcc_deduced_err(err, MccDefaultPCMErrorCode::ERROR_CCTE); } @@ -140,6 +140,8 @@ public: // _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()}); + _tableEpoch.emplace_back(hw_counts.epoch()); + _table.target_colon.emplace_back(target_coords.co_lon()); _table.target_colat.emplace_back(target_coords.co_lon()); @@ -153,11 +155,25 @@ public: } - MccError addPoint(mcc_skypoint_c auto target_coords, mcc_coord_pair_c auto const& hw_counts) + template + MccError getPoint(size_t idx, std::tuple& point) + requires(std::same_as && + std::same_as && + HW_T::pairKind == MccCoordPairKind::COORDS_KIND_XY) { - auto ref_epoch = hw_counts.epoch(); + if (idx > _table.target_colon.size()) { + return MccDefaultPCMConstructorErrorCode::ERROR_INVALID_INDEX; + } - return addPoint(std::move(target_coords), hw_counts, ref_epoch); + std::get<0>(point).setX(_table.target_colon[idx]); + std::get<0>(point).setY(_table.target_colat[idx]); + std::get<0>(point).setEpoch(_tableEpoch[idx]); + + std::get<1>(point).setX(_table.hw_colon[idx]); + std::get<1>(point).setY(_table.hw_colat[idx]); + std::get<1>(point).setEpoch(_tableEpoch[idx]); + + return MccDefaultPCMConstructorErrorCode::ERROR_OK; } @@ -168,6 +184,8 @@ public: void deletePoints() { + _tableEpoch.clear(); + _table.target_colon.clear(); _table.target_colat.clear(); @@ -285,13 +303,13 @@ public: // 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()); + std::vector xres(numberOfPoints()), yres(numberOfPoints()); // 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]; + for (size_t i = 0; i < numberOfPoints(); ++i) { + xres = _table.target_colon[i] - result.; yres = _table.target_colat[i]; } } @@ -301,7 +319,8 @@ public: protected: // std::vector _table; - table_t _table; + table_t _table{}; + std::vector _tableEpoch{}; compute_result_t bsplineFitting(MccDefaultPCM::pcm_data_t& pcm_data) { @@ -320,30 +339,54 @@ protected: 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, + 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; + } + + + 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(pcm_data.bspline.knotsY, pcm_data.bspline.knotsX, pcm_data.bspline.coeffsY, + _table.hw_colat, _table.hw_colon, result.model_colat); // get fitted residuals!!! + + + result.colon_res.resize(numberOfPoints()); + result.colat_res.resize(numberOfPoints()); + 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 + } + if (pcm_data.type == MccDefaultPCMType::PCM_TYPE_BSPLINE) { // 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) { @@ -366,8 +409,30 @@ protected: 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, + 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, + result.inv_model_colat); // get fitted residuals!!! + + result.inv_colon_res.resize(numberOfPoints()); + result.inv_colat_res.resize(numberOfPoints()); + 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 + } } + + return result; } + + + void robustLinearRegress() {} }; } // namespace mcc::impl \ No newline at end of file