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-18 23:34:46 +03:00
parent 40d36545f3
commit 9ad5bfb09b
2 changed files with 287 additions and 41 deletions
Download Patch File Download Diff File Expand all files Collapse all files

33
include/mcc/mcc_pcm.h
View File

@@ -166,6 +166,9 @@ public:
std::vector<coeff_t> coeffsX{};
std::vector<coeff_t> coeffsY{};
std::vector<coeff_t> inverseCoeffsX{};
std::vector<coeff_t> 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<DXT>) {
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<DYT>) {
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<double>::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<double>::quiet_NaN();

295
include/mcc/mcc_pcm_construct.h
View File

@@ -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<MccDefaultPCMConstructorErrorCode>(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<int>(ec), MccDefaultPCMConstructorErrorCategory::get());
}
} // namespace mcc::impl
namespace std
{
template <>
class is_error_code_enum<mcc::impl::MccDefaultPCMConstructorErrorCode> : public true_type
{
};
} // namespace std
namespace mcc::impl
{
template <MccMountType MOUNT_TYPE>
class MccPcmRefPointTable
class MccPCMConstructor
{
public:
using ref_coordpair_t = std::conditional_t<mcc_is_equatorial_mount<MOUNT_TYPE>,
@@ -34,6 +104,29 @@ public:
double colon_res, colat_res; // target - hw
};
struct table_t {
std::vector<double> target_colon, target_colat;
std::vector<double> hw_colon, hw_colat;
std::vector<double> colon_res, colat_res; // target - hw
};
struct compute_result_t {
MccDefaultPCMType pcm_type;
MccError error{}; // final model computation error
std::vector<double> model_colon, model_colat; // fitted model values
std::vector<double> 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<double> inv_model_colon, inv_model_colat; // fitted inverse model values
std::vector<double> 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<MOUNT_TYPE>::pcm_data_t& pcm_data)
compute_result_t computeModel(MccDefaultPCM<MOUNT_TYPE>::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<double> 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<double> colon_res = _table.colon_res;
std::vector<double> 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<double> 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_elem_t> _table;
table_t _table;
compute_result_t bsplineFitting(MccDefaultPCM<MOUNT_TYPE>::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<double> 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<double> 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<double> colon_res = _table.colon_res;
std::vector<double> 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_elem_t> _table;
};
} // namespace mcc::impl