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-03-27 17:29:46 +03:00
parent d4721d7b67
commit 8116658e0d
1 changed files with 52 additions and 28 deletions
Download Patch File Download Diff File Expand all files Collapse all files

80
include/mcc/mcc_pcm_construct.h
View File

@@ -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<double> tx(pcm_data.bspline.knotsX.size() + 6), ty(pcm_data.bspline.knotsY.size() + 6);
std::vector<double> 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<double> theta_tag(numberOfPoints());
// inverse (encoder = celestial + pcm)
std::vector<double> 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
}
}