[ VIGRA Homepage | Function Index | Class Index | Namespaces | File List | Main Page ]

polynomial_registration.hxx
1/************************************************************************/
2/* */
3/* Copyright 2007-2013 by Benjamin Seppke */
4/* */
5/* This file is part of the VIGRA computer vision library. */
6/* The VIGRA Website is */
7/* http://hci.iwr.uni-heidelberg.de/vigra/ */
8/* Please direct questions, bug reports, and contributions to */
9/* ullrich.koethe@iwr.uni-heidelberg.de or */
10/* vigra@informatik.uni-hamburg.de */
11/* */
12/* Permission is hereby granted, free of charge, to any person */
13/* obtaining a copy of this software and associated documentation */
14/* files (the "Software"), to deal in the Software without */
15/* restriction, including without limitation the rights to use, */
16/* copy, modify, merge, publish, distribute, sublicense, and/or */
17/* sell copies of the Software, and to permit persons to whom the */
18/* Software is furnished to do so, subject to the following */
19/* conditions: */
20/* */
21/* The above copyright notice and this permission notice shall be */
22/* included in all copies or substantial portions of the */
23/* Software. */
24/* */
25/* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND */
26/* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES */
27/* OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND */
28/* NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT */
29/* HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, */
30/* WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING */
31/* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR */
32/* OTHER DEALINGS IN THE SOFTWARE. */
33/* */
34/************************************************************************/
35
36#ifndef VIGRA_POLYNOMIAL_REGISTRATION_HXX
37#define VIGRA_POLYNOMIAL_REGISTRATION_HXX
38
39#include "mathutil.hxx"
40#include "matrix.hxx"
41#include "linear_solve.hxx"
42#include "tinyvector.hxx"
43#include "splineimageview.hxx"
44
45namespace vigra
46{
47
48/** \addtogroup Registration Image Registration
49 */
50//@{
51
52/**
53 * Iterative function for determinination of the polynom weights:
54 *
55 * Example: order=2, x, y
56 * ----->
57 * [1,
58 * x, y,
59 * x^2, x*y, y^2]
60 *
61 * This function is needed, because the polynomial transformation Matrix
62 * has the the same number of rows. the target position is then determined
63 * by multiplying each x- and y-transformation result value with the
64 * corresponding weight for the current x- and y-coordinate, given by this
65 * function.
66 */
67std::vector<double> polynomialWarpWeights(double x, double y, unsigned int polynom_order)
68{
69 unsigned int poly_count = (polynom_order+1)*(polynom_order+2)/2;
70
71 std::vector<double> weights(poly_count);
72
73 unsigned int weight_idx=0;
74
75 for (unsigned int order=0; order<=polynom_order; order++)
76 {
77 for(unsigned int i=0; i<=order; i++, weight_idx++)
78 {
79 weights[weight_idx] = pow(x,(double)order-i)*pow(y,(double)i);
80 }
81 }
82 return weights;
83}
84
85/********************************************************/
86/* */
87/* polynomialMatrix2DFromCorrespondingPoints */
88/* */
89/********************************************************/
90
91/** \brief Create polynomial matrix of a certain degree that maps corresponding points onto each other.
92
93 For use with \ref polynomialWarpImage() of same degree.
94
95 Since polynoms are usually non-linear functions, a special semantics is embedded to define
96 a matrix here. Each matrix consist of two rows, containing x- and y-factors of the polynom.
97
98 The meaning of the matrix is explained at the example of a polynom of 2nd order:
99
100 First Row = [a_x b_x c_x d_x e_x f_x]
101 Second Row = [a_y b_y c_y d_y e_y f_y]
102
103 The transformed coordinate p'=[x' y'] of a position p=[x y] is then:
104
105 x' = a_x + b_x*x + c_x*y + d_x*x^2 + e_x*x*y + f_x*y^2
106 y' = a_y + b_y*x + c_y*y + d_y*x^2 + e_y*x*y + f_y*y^2
107
108 Note that the order of the polynom's factors is directly influenced by the
109 \ref polynomialWarpWeights() function and follows the intuitive scheme.
110*/
111template <int PolynomOrder,
112 class SrcPointIterator,
113 class DestPointIterator>
114linalg::TemporaryMatrix<double>
117{
118 int point_count = s_end - s;
119 int poly_count = (PolynomOrder+1)*(PolynomOrder+2)/2;
120
122 std::vector<double> weights;
123
124 for (int i =0; i<point_count; ++i, ++s, ++d)
125 {
126 weights = polynomialWarpWeights((*d)[0], (*d)[1], PolynomOrder);
127
128 for(int c=0; c<poly_count; c++)
129 {
130 A(i,c) = weights[c];
131 }
132
133 b1(i,0)=(*s)[0];b2(i,0)=(*s)[1];
134 }
135
136 if(!vigra::linearSolve( A, b1, res1 ) || !vigra::linearSolve( A, b2, res2 ))
137 vigra_fail("polynomialMatrix2DFromCorrespondingPoints(): singular solution matrix.");
138
140
141 for(int c=0; c<poly_count; c++)
142 {
143 res(c,0) = res1(c,0);
144 res(c,1) = res2(c,0);
145 }
146
147 return res;
148}
149
150
151/********************************************************/
152/* */
153/* polynomialWarpImage */
154/* */
155/********************************************************/
156
157/** \brief Warp an image according to an polynomial transformation.
158
159 To get more information about the structure of the matrix,
160 see \ref polynomialMatrix2DFromCorrespondingPoints().
161
162 <b>\#include</b> <vigra/polynomial_registration.hxx><br>
163 Namespace: vigra
164
165 pass 2D array views:
166 \code
167 namespace vigra {
168 template <int ORDER, class T,
169 class T2, class S2,
170 class C>
171 void
172 polynomialWarpImage(SplineImageView<ORDER, T> const & src,
173 MultiArrayView<2, T2, S2> dest,
174 MultiArrayView<2, double, C> const & polynomialMatrix);
175 }
176 \endcode
177
178 \deprecatedAPI{polynomialWarpImage}
179
180 pass arguments explicitly:
181 \code
182 namespace vigra {
183 template <int ORDER, class T,
184 class DestIterator, class DestAccessor,
185 class C>
186 void polynomialWarpImage(SplineImageView<ORDER, T> const & src,
187 DestIterator dul, DestIterator dlr, DestAccessor dest,
188 MultiArrayView<2, double, C> const & polynomialMatrix);
189 }
190 \endcode
191
192 use argument objects in conjunction with \ref ArgumentObjectFactories :
193 \code
194 namespace vigra {
195 template <int ORDER, class T,
196 class DestIterator, class DestAccessor,
197 class C>
198 void polynomialWarpImage(SplineImageView<ORDER, T> const & src,
199 triple<DestIterator, DestIterator, DestAccessor> dest,
200 MultiArrayView<2, double, C> const & polynomialMatrix);
201 }
202 \endcode
203 \deprecatedEnd
204 */
205doxygen_overloaded_function(template <...> void polynomialWarpImage)
206
207template <int PolynomOrder,
208 int ORDER, class T,
209 class DestIterator, class DestAccessor,
210 class C>
214{
215 int poly_count = (PolynomOrder+1)*(PolynomOrder+2)/2;
216
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.");
219
220 double w = dlr.x - dul.x;
221 double h = dlr.y - dul.y;
222
223 std::vector<double> weights(poly_count);
224
225 for(double y = 0.0; y < h; ++y, ++dul.y)
226 {
227 typename DestIterator::row_iterator rd = dul.rowIterator();
228 for(double x=0.0; x < w; ++x, ++rd)
229 {
230 weights = polynomialWarpWeights(x,y, PolynomOrder);
231
232 double sx=0;
233 double sy=0;
234
235 for(int c=0; c<poly_count; c++)
236 {
237 sx += weights[c]*polynomialMatrix(c,0);
238 sy += weights[c]*polynomialMatrix(c,1);
239 }
240
241 if(src.isInside(sx, sy))
242 dest.set(src(sx, sy), rd);
243 }
244 }
245}
246
247template <int PolynomOrder,
248 int ORDER, class T,
249 class DestIterator, class DestAccessor,
250 class C>
251inline
252void polynomialWarpImage(SplineImageView<ORDER, T> const & src,
253 triple<DestIterator, DestIterator, DestAccessor> dest,
254 MultiArrayView<2, double, C> const & polynomialMatrix)
255{
256 polynomialWarpImage<PolynomOrder>(src, dest.first, dest.second, dest.third, polynomialMatrix);
257}
258
259
260template <int PolynomOrder,
261 int ORDER, class T,
262 class T2, class S2,
263 class C>
264inline
265void polynomialWarpImage(SplineImageView<ORDER, T> const & src,
266 MultiArrayView<2, T2, S2> dest,
267 MultiArrayView<2, double, C> const & polynomialMatrix)
268{
269 polynomialWarpImage<PolynomOrder>(src, destImageRange(dest), polynomialMatrix);
270}
271
272
273//@}
274
275} // namespace vigra
276
277
278#endif /* VIGRA_POLYNOMIAL_REGISTRATION_HXX */
Class for a single RGB value.
Definition rgbvalue.hxx:128
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
void polynomialWarpImage(...)
Warp an image according to an polynomial transformation.

© Ullrich Köthe (ullrich.koethe@iwr.uni-heidelberg.de)
Heidelberg Collaboratory for Image Processing, University of Heidelberg, Germany

html generated using doxygen and Python
vigra 1.12.2