timur/mcc
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

...

Browse Source
This commit is contained in:
Timur A. Fatkhullin
2026-04-16 18:42:11 +03:00
parent a8a1d9a763
commit 3373bedb39
4 changed files with 268 additions and 84 deletions
Download Patch File Download Diff File Expand all files Collapse all files

11
include/mcc/mcc_pcm.h
View File

@@ -485,6 +485,17 @@ private:
}
}
if constexpr (mcc_is_equatorial_mount<MOUNT_TYPE>) { // HA - DEC system (needs normalization)
if (!inverse) {
x = MccAngleHA_OBS(x);
y = MccAngleDEC_OBS(y);
}
// in the case of inverse PCM it is assumed that angles are already normalized (inputs are celestial
// coordiinates)
}
pcm_geom_coeffs_t* geom_coeffs = &_pcmData.geomCoefficients;
#ifdef USE_BSPLINE_PCM

226
include/mcc/mcc_pcm_fit.h
View File

@@ -112,9 +112,9 @@ public:
size_t final_iter{}; // number of final iteration
std::vector<double> model_colon{}, model_colat{}; // fitting model values
std::vector<double> colon_res{}, colat_res{}; // target - model
std::vector<double> colon_weight{}, colat_weight; // Tukey's weights
std::vector<double> model_colonRES{}, model_colatRES{}; // fitting model values: fit(target-hw)
std::vector<double> colon_err{}, colat_err{}; // fitting residuals: (target-hw) - fit
std::vector<double> colon_weight{}, colat_weight; // Tukey's weights
#ifdef USE_BSPLINE_PCM
std::array<int, 4> bspline_fit_err{
@@ -122,8 +122,8 @@ public:
// 0 - pcmX(X,Y), 1 - pcmY(X,Y), 2 - revPcmX(COLON, COLAT), 3 - revPcmY(COLON, COLAT)
// quantities below are computed only fo pcm_type == MccDefaultPCMType::PCM_TYPE_BSPLINE
std::vector<double> inv_model_colon{}, inv_model_colat{}; // fitted inverse model values
std::vector<double> inv_colon_res{}, inv_colat_res{}; // encoder - model
std::vector<double> inv_model_colonRES{}, inv_model_colatRES{}; // fitted inverse model values
std::vector<double> inv_colon_err{}, inv_colat_err{}; // encoder - model
#endif
};
@@ -134,6 +134,15 @@ public:
};
struct pcm_table_elem_t {
ref_coordpair_t target{}; // celestial
MccGenXY hw{}; // hardware (encoder)
MccGenXY res{}; // residuals: celestial - hardware
};
typedef std::vector<pcm_table_elem_t> pcm_table_t;
size_t numberOfPoints() const
{
return _targetCOLON.size();
@@ -206,21 +215,31 @@ public:
}
// something like:
// auto [tag, hw, err] = getPoint(7);
template <typename COLON_T, typename COLAT_T, typename HW_X_T, typename HW_Y_T>
MccError addPoint(COLON_T const& colon,
COLAT_T const& colat,
HW_X_T const& x,
HW_Y_T const& y,
mcc_coord_epoch_c auto const& epoch)
requires(std::is_arithmetic_v<COLON_T> && std::is_arithmetic_v<COLAT_T> && std::is_arithmetic_v<HW_X_T> &&
std::is_arithmetic_v<HW_Y_T>)
{
ref_coordpair_t cp{typename ref_coordpair_t::x_t{colon}, typename ref_coordpair_t::y_t{colat}, epoch};
MccGenXY hw{x, y, epoch};
return addPoint(cp, hw);
}
template <mcc_coord_pair_c TAG_T, mcc_coord_pair_c HW_T>
std::tuple<TAG_T, HW_T, MccError> getPoint(size_t idx)
MccError getPoint(size_t idx, std::tuple<TAG_T, HW_T>& point)
requires(std::same_as<typename TAG_T::x_t, typename ref_coordpair_t::x_t> &&
std::same_as<typename TAG_T::y_t, typename ref_coordpair_t::y_t> &&
HW_T::pairKind == MccCoordPairKind::COORDS_KIND_XY)
{
std::tuple<TAG_T, HW_T, MccError> point;
if (idx > numberOfPoints()) {
std::get<2>(point) = MccPCMFitterErrorCode::ERROR_INVALID_INDEX;
return MccPCMFitterErrorCode::ERROR_INVALID_INDEX;
} else {
std::get<2>(point) = MccPCMFitterErrorCode::ERROR_OK;
std::get<0>(point).setX(_targetCOLON[idx]);
std::get<0>(point).setY(_targetCOLAT[idx]);
std::get<0>(point).setEpoch(_epoch[idx]);
@@ -230,7 +249,103 @@ public:
std::get<1>(point).setEpoch(_epoch[idx]);
}
return point;
return MccPCMFitterErrorCode::ERROR_OK;
}
pcm_table_t getPCMTable() const
{
pcm_table_t tab;
if (numberOfPoints()) {
for (size_t i = 0; i < numberOfPoints(); ++i) {
tab.emplace_back(pcm_table_elem_t{.target{typename ref_coordpair_t::x_t{_targetCOLON[i]},
typename ref_coordpair_t::y_t{_targetCOLAT[i]}, _epoch[i]},
.hw{MccAngleX{_hwX[i]}, MccAngleY{_hwY[i]}, _epoch[i]},
.res{MccAngleX{_colonRES[i]}, MccAngleY{_colatRES[i]}, _epoch[i]}});
}
}
return tab;
}
//
// 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
// it are interpretated as border knots and final full knots set is:
// knots = [input_knots[0], input_knots[0], input_knots[0], input_knots[0], input_knots[1], input_knots[2],
// ..., input_knots[N-1], input_knots[N-1], input_knots[N-1], input_knots[N-1]], where N = input_knots.size()
//
// WARNING: the input knots for inverse B-spline are ignored so the direct and inverse B-spline coefficients are
// calculated on the same mesh!
compute_result_t computeModel(MccDefaultPCM<MOUNT_TYPE>::pcm_data_t& pcm_data,
compute_params_t const& comp_params = {})
{
compute_result_t result{.pcm_type = pcm_data.type, .error = MccPCMFitterErrorCode::ERROR_OK};
size_t min_data_size = 2; // 2 is for BSPLINE
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 {
min_data_size = 8;
}
if (numberOfPoints() < min_data_size) {
result.error = MccPCMFitterErrorCode::ERROR_NOT_ENOUGH_DATA;
return result;
}
// robust linear regression with Tukey's loss function
result = robustLinearRegress(pcm_data, comp_params);
} else {
if (numberOfPoints() < min_data_size) {
result.error = MccPCMFitterErrorCode::ERROR_NOT_ENOUGH_DATA;
return result;
}
}
#ifdef USE_BSPLINE_PCM
if (pcm_data.type == MccDefaultPCMType::PCM_TYPE_BSPLINE) {
return bsplineFitting(pcm_data);
} 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
std::vector<double> xres = _colonRES, yres = _colatRES;
for (size_t i = 0; i < numberOfPoints(); ++i) {
_colonRES[i] = _targetCOLON[i] - result.model_colonRES[i];
_colatRES[i] = _targetCOLAT[i] - result.model_colatRES[i];
}
auto r = bsplineFitting(pcm_data);
result.error = r.error;
result.bspline_fit_err = r.bspline_fit_err;
if (!r.error) {
for (size_t i = 0; i < numberOfPoints(); ++i) {
result.model_colonRES[i] += r.model_colonRES[i];
result.model_colatRES[i] += r.model_colatRES[i];
result.colon_err[i] = _targetCOLON[i] - result.model_colonRES[i];
result.colat_err[i] = _targetCOLAT[i] - result.model_colatRES[i];
}
}
// restore original residuals
_colonRES = xres;
_colatRES = yres;
}
#endif
return result;
}
protected:
@@ -242,21 +357,21 @@ protected:
void computeModelResi(compute_result_t& result)
{
result.colon_res.resize(numberOfPoints());
result.colat_res.resize(numberOfPoints());
result.colon_err.resize(numberOfPoints());
result.colat_err.resize(numberOfPoints());
for (size_t i = 0; i < numberOfPoints(); ++i) {
result.colon_res[i] = _colonRES[i] - result.model_colon[i]; // = target - model
result.colat_res[i] = _colatRES[i] - result.model_colat[i]; // = target - model
result.colon_err[i] = _colonRES[i] - result.model_colonRES[i]; // = target - model
result.colat_err[i] = _colatRES[i] - result.model_colatRES[i]; // = target - model
}
#ifdef USE_BSPLINE_PCM
if (result.pcm_type == MccDefaultPCMType::PCM_TYPE_BSPLINE) {
result.inv_colon_res.resize(numberOfPoints());
result.inv_colat_res.resize(numberOfPoints());
result.inv_colon_err.resize(numberOfPoints());
result.inv_colat_err.resize(numberOfPoints());
for (size_t i = 0; i < numberOfPoints(); ++i) {
result.inv_colon_res[i] = _colonRES[i] - result.inv_model_colon[i]; // = hw - model
result.inv_colat_res[i] = _colatRES[i] - result.inv_model_colat[i]; // = hw - model
result.inv_colon_err[i] = _colonRES[i] - result.inv_model_colonRES[i]; // = hw - model
result.inv_colat_err[i] = _colatRES[i] - result.inv_model_colatRES[i]; // = hw - model
}
}
#endif
@@ -279,9 +394,25 @@ protected:
// 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<double> tx(pcm_data.bspline.knotsX.size() + 6), ty(pcm_data.bspline.knotsY.size() + 6);
std::vector<double>&theta_hw = _hwY, theta_tag = _targetCOLAT, phi_hw = _hwX, phi_tag = _targetCOLON;
std::vector<double>&theta_hw = _hwY, &theta_tag = _targetCOLAT, &phi_hw = _hwX, &phi_tag = _targetCOLON;
std::vector<double> theta, phi;
// full set of B-spline knots
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];
ty[ty.size() - 1] = ty[ty.size() - 2] = ty[ty.size() - 3] =
pcm_data.bspline.knotsY[pcm_data.bspline.knotsY.size() - 1];
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];
}
/*
WARNING:
FITPACK B-spline on sphere: in the fitting routines the first angle argument is THETA - co-latitude
@@ -303,21 +434,14 @@ protected:
theta_hw[i] = _hwY[i] + MCC_HALF_PI;
phi_hw[i] = _hwX[i] + std::numbers::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];
for (size_t i = 0; i < tx.size(); ++i) {
tx[i] += std::numbers::pi;
}
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;
for (size_t i = 0; i < ty.size(); ++i) {
ty[i] += MCC_HALF_PI;
}
}
size_t Ncoeffs = (tx.size() - 4) * (ty.size() - 4);
@@ -343,10 +467,10 @@ protected:
bsplines::fitpack_eval_spl2d(ty, tx, pcm_data.bspline.coeffsX, theta_hw, phi_hw,
result.model_colon); // get fitted residuals!!!
result.model_colonRES); // get fitted residuals!!!
bsplines::fitpack_eval_spl2d(ty, tx, pcm_data.bspline.coeffsY, theta_hw, phi_hw,
result.model_colat); // get fitted residuals!!!
result.model_colatRES); // get fitted residuals!!!
if (pcm_data.type == MccDefaultPCMType::PCM_TYPE_BSPLINE) {
@@ -358,16 +482,16 @@ protected:
// inverse (encoder = celestial + pcm)
std::vector<double> colon_res = _colonRES;
std::vector<double> colat_res = _colatRES;
std::vector<double> colon_err = _colonRES;
std::vector<double> colat_err = _colatRES;
if constexpr (std::same_as<ref_coordpair_t, MccSkyHADEC_OBS>) {
theta_tag = theta;
phi_tag = phi;
}
for (size_t i = 0; i < colat_res.size(); ++i) {
colon_res[i] = -colon_res[i];
colat_res[i] = -colat_res[i];
for (size_t i = 0; i < colat_err.size(); ++i) {
colon_err[i] = -colon_err[i];
colat_err[i] = -colat_err[i];
if constexpr (std::same_as<ref_coordpair_t, MccSkyHADEC_OBS>) {
theta_tag[i] = _targetCOLAT[i] + MCC_HALF_PI;
@@ -375,14 +499,14 @@ protected:
}
}
result.bspline_fit_err[2] = bsplines::fitpack_sphere_fit(theta_tag, phi_tag, colon_res, 1.0, ty, tx,
result.bspline_fit_err[2] = bsplines::fitpack_sphere_fit(theta_tag, phi_tag, colon_err, 1.0, ty, tx,
pcm_data.bspline.inverseCoeffsX, resi2x);
if (result.bspline_fit_err[2] > 0) {
result.error = MccPCMFitterErrorCode::ERROR_BSPLINE_FIT;
return result;
}
result.bspline_fit_err[3] = bsplines::fitpack_sphere_fit(theta_tag, phi_tag, colat_res, 1.0, ty, tx,
result.bspline_fit_err[3] = bsplines::fitpack_sphere_fit(theta_tag, phi_tag, colat_err, 1.0, ty, tx,
pcm_data.bspline.inverseCoeffsY, resi2y);
if (result.bspline_fit_err[3] > 0) {
result.error = MccPCMFitterErrorCode::ERROR_BSPLINE_FIT;
@@ -390,10 +514,10 @@ protected:
}
bsplines::fitpack_eval_spl2d(ty, tx, pcm_data.bspline.inverseCoeffsX, theta_tag, phi_tag,
result.inv_model_colon); // get fitted residuals!!!
result.inv_model_colonRES); // get fitted residuals!!!
bsplines::fitpack_eval_spl2d(ty, tx, pcm_data.bspline.inverseCoeffsY, theta_tag, phi_tag,
result.inv_model_colat); // get fitted residuals!!!
result.inv_model_colatRES); // get fitted residuals!!!
}
computeModelResi(result);
@@ -546,10 +670,10 @@ protected:
pcm_data.geomCoefficients.forkFlexure = new_beta(8);
}
result.model_colon = {model.begin(), model.begin() + numberOfPoints()};
result.model_colat = {model.begin() + numberOfPoints(), model.end()};
result.colon_res.resize(numberOfPoints());
result.colat_res.resize(numberOfPoints());
result.model_colonRES = {model.begin(), model.begin() + numberOfPoints()};
result.model_colatRES = {model.begin() + numberOfPoints(), model.end()};
result.colon_err.resize(numberOfPoints());
result.colat_err.resize(numberOfPoints());
result.colon_weight = {weights.begin(), weights.begin() + numberOfPoints()};
result.colat_weight = {weights.begin() + numberOfPoints(), weights.end()};