36#ifndef VIGRA_POLYNOMIAL_REGISTRATION_HXX
37#define VIGRA_POLYNOMIAL_REGISTRATION_HXX
39#include "mathutil.hxx"
41#include "linear_solve.hxx"
42#include "tinyvector.hxx"
43#include "splineimageview.hxx"
69 unsigned int poly_count = (polynom_order+1)*(polynom_order+2)/2;
71 std::vector<double> weights(poly_count);
73 unsigned int weight_idx=0;
75 for (
unsigned int order=0; order<=polynom_order; order++)
77 for(
unsigned int i=0; i<=order; i++, weight_idx++)
79 weights[weight_idx] = pow(x,(
double)order-i)*pow(y,(
double)i);
111template <
int PolynomOrder,
112 class SrcPointIterator,
113 class DestPointIterator>
114linalg::TemporaryMatrix<double>
118 int point_count = s_end - s;
119 int poly_count = (PolynomOrder+1)*(PolynomOrder+2)/2;
121 vigra::Matrix<double> A(point_count,poly_count), b1(point_count,1), res1(poly_count,1), b2(point_count,1), res2(poly_count,1);
122 std::vector<double> weights;
124 for (
int i =0; i<point_count; ++i, ++s, ++d)
128 for(
int c=0; c<poly_count; c++)
133 b1(i,0)=(*s)[0];b2(i,0)=(*s)[1];
137 vigra_fail(
"polynomialMatrix2DFromCorrespondingPoints(): singular solution matrix.");
141 for(
int c=0; c<poly_count; c++)
143 res(c,0) = res1(c,0);
144 res(c,1) = res2(c,0);
207template <
int PolynomOrder,
209 class DestIterator,
class DestAccessor,
212 DestIterator dul, DestIterator dlr, DestAccessor dest,
215 int poly_count = (PolynomOrder+1)*(PolynomOrder+2)/2;
217 vigra_precondition(
rowCount(polynomialMatrix) == poly_count &&
columnCount(polynomialMatrix) == 2,
218 "polynomialWarpImage(): matrix doesn't represent a polynomial transformation of given degreee in 2D coordinates.");
220 double w = dlr.x - dul.x;
221 double h = dlr.y - dul.y;
223 std::vector<double> weights(poly_count);
225 for(
double y = 0.0; y < h; ++y, ++dul.y)
227 typename DestIterator::row_iterator rd = dul.rowIterator();
228 for(
double x=0.0; x < w; ++x, ++rd)
235 for(
int c=0; c<poly_count; c++)
237 sx += weights[c]*polynomialMatrix(c,0);
238 sy += weights[c]*polynomialMatrix(c,1);
242 dest.set(src(sx, sy), rd);
247template <
int PolynomOrder,
249 class DestIterator,
class DestAccessor,
253 triple<DestIterator, DestIterator, DestAccessor> dest,
254 MultiArrayView<2, double, C>
const & polynomialMatrix)
260template <
int PolynomOrder,
Base class for, and view to, MultiArray.
Definition multi_array.hxx:705
Create a continuous view onto a discrete image using splines.
Definition splineimageview.hxx:100
bool isInside(double x, double y) const
Definition splineimageview.hxx:487
Definition matrix.hxx:125
MultiArrayIndex columnCount(const MultiArrayView< 2, T, C > &x)
Definition matrix.hxx:684
linalg::TemporaryMatrix< double > polynomialMatrix2DFromCorrespondingPoints(SrcPointIterator s, SrcPointIterator s_end, DestPointIterator d)
Create polynomial matrix of a certain degree that maps corresponding points onto each other.
Definition polynomial_registration.hxx:115
std::vector< double > polynomialWarpWeights(double x, double y, unsigned int polynom_order)
Definition polynomial_registration.hxx:67
MultiArrayIndex rowCount(const MultiArrayView< 2, T, C > &x)
Definition matrix.hxx:671
void polynomialWarpImage(...)
Warp an image according to an polynomial transformation.