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

rbf_registration.hxx |

1 /************************************************************************/
2 /* */
3 /* Copyright 2007-2014 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_RBF_REGISTRATION_HXX
37 #define VIGRA_RBF_REGISTRATION_HXX
38
39 #include <vigra/mathutil.hxx>
40 #include <vigra/matrix.hxx>
41 #include <vigra/linear_solve.hxx>
42 #include <vigra/tinyvector.hxx>
43 #include <vigra/splineimageview.hxx>
44 #include <vigra/affine_registration.hxx>
45
46 namespace vigra {
47
48 /** \addtogroup Registration Image Registration
49 */
50 //@{
51
52 namespace detail {
53
54 // Small and hopefully fast helper function to compute the squared-distance
55 // between two points (2d)
56 template <class SrcPoint, class DestPoint>
57 inline double distance2(SrcPoint const & p1, DestPoint const & p2)
58 {
59 return (double)(p1[0]-p2[0])*(p1[0]-p2[0]) + (p1[1]-p2[1])*(p1[1]-p2[1]);
60 }
61
62 } //end of namespace vigra::detail
63
64
65 /**
66 * The famous thin plate spline functor [weight(d) = d^2*log(d^2)]
67 * only affects up to max distance...
68 */
69 struct ThinPlateSplineFunctor
70 {
71 template <class SrcPoint, class DestPoint>
72 inline double operator()(SrcPoint const & p1, DestPoint const & p2) const
73 {
74 double dist2 = detail::distance2(p1, p2);
75
76 if(dist2 == 0 )
77 {
78 return 0;
79 }
80 else
81 {
82 return dist2*log(dist2);
83 }
84 }
85 };
86
87 /**
88 * A distance power based radial basis functor [weight = dist^N]
89 */
90 template<int N>
91 struct DistancePowerFunctor
92 {
93 template <class SrcPoint, class DestPoint>
94 inline double operator()(SrcPoint const & p1, DestPoint const & p2) const
95 {
96 double dist2 = detail::distance2(p1, p2);
97
98 if(dist2 == 0)
99 {
100 return 0;
101 }
102 else
103 {
104 return pow(dist2, N/2.0);
105 }
106 }
107 };
108
109 /********************************************************/
110 /* */
111 /* rbfMatrix2DFromCorrespondingPoints */
112 /* */
113 /********************************************************/
114
115 /** \brief Create a matrix that maps corresponding points onto each other using a given RBF.
116
117 For use with \ref rbfWarpImage(). For n given (corresponding) points,
118 the matrix will be of size (n+3,2). Note that the representation of this matrix is exactly
119 the same as the "W" matrix of Bookstein. More information can be found in the following article:
120
121 Fred L. Bookstein. Principal Warps: Thin-Plate Splines and the Decomposition of Deformations. IEEE PAMI, Vol 11, No 8. 1989
122 */
123 template <class RadialBasisFunctor,
124 class SrcPointIterator, class DestPointIterator>
125 linalg::TemporaryMatrix<double>
126 rbfMatrix2DFromCorrespondingPoints(SrcPointIterator s, SrcPointIterator s_end, DestPointIterator d, RadialBasisFunctor const & rbf)
127 {
128 int point_count = s_end - s;
129
130 Matrix<double> L(point_count+3, point_count+3, 0.0);
131 Matrix<double> Y(point_count+3, 2, 0.0);
132 Matrix<double> W(point_count+3, 2, 0.0);
133
134 //fill P (directly into K) and V (directly into Y)
135 for(int i=0; i<point_count; ++i)
136 {
137 L(i, point_count ) = L(point_count, i) = 1;
138 L(i, point_count+1) = L(point_count+1, i) = (d[i])[0];
139 L(i, point_count+2) = L(point_count+2, i) = (d[i])[1];
140
141 Y(i,0)= (s[i])[0];
142 Y(i,1)= (s[i])[1];
143 }
144
145 //fill K (directly into L)
146 for(int j=0; j<point_count; j++)
147 {
148 for(int i=0; i<j; i++)
149 {
150 L(i,j) = L(j,i) = rbf(d[i], d[j]);
151 }
152 }
153
154 linearSolve(L, Y, W);
155 //Results are okay, even if vigra reports failure...
156 //so I commented this out
157 // if(!linearSolve(L, Y, W))
158 // vigra_fail("radialBasisMatrix2DFromCorrespondingPoints(): singular solution matrix.");
159
160 return W;
161 };
162
163
164 /********************************************************/
165 /* */
166 /* rbfWarpImage */
167 /* */
168 /********************************************************/
169
170 /** \brief Warp an image according to an radial basis function based transformation.
171
172 To get more information about the structure of the matrix, see \ref rbfMatrix2DFromCorrespondingPoints()
173
174 <b>\#include</b> <vigra/rbf_registration.hxx><br>
175 Namespace: vigra
176
177 pass 2D array views:
178 \code
179 namespace vigra {
180 template <int ORDER, class T,
181 class T2, class S2,
182 class DestPointIterator,
183 class C,
184 class RadialBasisFunctor>
185 void
186 rbfWarpImage(SplineImageView<ORDER, T> const & src,
187 MultiArrayView<2, T2, S2> dest,
188 DestPointIterator d, DestPointIterator d_end,
189 MultiArrayView<2, double, C> const & W,
190 RadialBasisFunctor rbf);
191 }
192 \endcode
193
194 \deprecatedAPI{rbfWarpImage}
195
196 pass arguments explicitly:
197 \code
198 namespace vigra {
199 template <int ORDER, class T,
200 class DestIterator, class DestAccessor,
201 class DestPointIterator,
202 class C,
203 class RadialBasisFunctor>
204 void
205 rbfWarpImage(SplineImageView<ORDER, T> const & src,
206 DestIterator dul, DestIterator dlr, DestAccessor dest,
207 DestPointIterator d, DestPointIterator d_end,
208 MultiArrayView<2, double, C> const & W,
209 RadialBasisFunctor rbf);
210 }
211 \endcode
212
213 use argument objects in conjunction with \ref ArgumentObjectFactories :
214 \code
215 namespace vigra {
216 template <int ORDER, class T,
217 class DestIterator, class DestAccessor,
218 class DestPointIterator,
219 class C,
220 class RadialBasisFunctor>
221 void
222 rbfWarpImage(SplineImageView<ORDER, T> const & src,
223 triple<DestIterator, DestIterator, DestAccessor> dest,
224 DestPointIterator d, DestPointIterator d_end,
225 MultiArrayView<2, double, C> const & W,
226 RadialBasisFunctor rbf);
227 }
228 }
229 \endcode
230 \deprecatedEnd
231 */
232 template <int ORDER, class T,
233 class DestIterator, class DestAccessor,
234 class DestPointIterator,
235 class C,
236 class RadialBasisFunctor>
237 void
238 rbfWarpImage(SplineImageView<ORDER, T> const & src,
239 DestIterator dul, DestIterator dlr, DestAccessor dest,
240 DestPointIterator d, DestPointIterator d_end,
241 MultiArrayView<2, double, C> const & W,
242 RadialBasisFunctor rbf)
243 {
244 int point_count = d_end - d;
245
246 vigra_precondition(rowCount(W) == point_count+3 && columnCount(W) == 2,
247 "vigra::rbfWarpImage(): matrix doesn't represent a proper transformation of given point size in 2D coordinates.");
248
249 double w = dlr.x - dul.x;
250 double h = dlr.y - dul.y;
251
252 for(double y = 0.0; y < h; ++y, ++dul.y)
253 {
254 typename DestIterator::row_iterator rd = dul.rowIterator();
255 for(double x=0.0; x < w; ++x, ++rd)
256 {
257 //Affine part
258 double sx = W(point_count,0)+W(point_count+1,0)*x+ W(point_count+2,0)*y,
259 sy = W(point_count,1)+W(point_count+1,1)*x+ W(point_count+2,1)*y;
260
261 //RBS part
262 for(int i=0; i<point_count; i++)
263 {
264 double weight = rbf(d[i], Diff2D(x,y));
265 sx += W(i,0)*weight;
266 sy += W(i,1)*weight;
267
268 }
269
270 if(src.isInside(sx, sy))
271 dest.set(src(sx, sy), rd);
272 }
273 }
274 };
275
276 template <int ORDER, class T,
277 class DestIterator, class DestAccessor,
278 class DestPointIterator,
279 class C,
280 class RadialBasisFunctor>
281 inline
282 void rbfWarpImage(SplineImageView<ORDER, T> const & src,
283 triple<DestIterator, DestIterator, DestAccessor> dest,
284 DestPointIterator d, DestPointIterator d_end,
285 MultiArrayView<2, double, C> const & W,
286 RadialBasisFunctor rbf)
287 {
288 rbfWarpImage(src, dest.first, dest.second, dest.third, d, d_end, W, rbf);
289 }
290
291 template <int ORDER, class T,
292 class T2, class S2,
293 class DestPointIterator,
294 class C,
295 class RadialBasisFunctor>
296 inline
297 void rbfWarpImage(SplineImageView<ORDER, T> const & src,
298 MultiArrayView<2, T2, S2> dest,
299 DestPointIterator d, DestPointIterator d_end,
300 MultiArrayView<2, double, C> const & W,
301 RadialBasisFunctor rbf)
302 {
303 rbfWarpImage(src, destImageRange(dest), d, d_end, W, rbf);
304 }
305
306
307 //@}
308
309 } // namespace vigra
310
311
312 #endif /* VIGRA_RBF_REGISTRATION_HXX */
vigra::DistancePowerFunctor
Definition: rbf_registration.hxx:91
vigra::linalg::rowCount
MultiArrayIndex rowCount(const MultiArrayView< 2, T, C > &x)
Definition: matrix.hxx:669
vigra::rbfWarpImage
void rbfWarpImage(SplineImageView< ORDER, T > const &src, DestIterator dul, DestIterator dlr, DestAccessor dest, DestPointIterator d, DestPointIterator d_end, MultiArrayView< 2, double, C > const &W, RadialBasisFunctor rbf)
Warp an image according to an radial basis function based transformation.
Definition: rbf_registration.hxx:238
vigra::linalg::Matrix
Definition: matrix.hxx:121
vigra::ThinPlateSplineFunctor
Definition: rbf_registration.hxx:69
vigra::Diff2D
Two dimensional difference vector.
Definition: diff2d.hxx:185
vigra
Definition: accessor.hxx:43
vigra::SplineImageView::isInside
bool isInside(double x, double y) const
Definition: splineimageview.hxx:487
vigra::linalg::log
linalg::TemporaryMatrix< T > log(MultiArrayView< 2, T, C > const &v)
vigra::linalg::columnCount
MultiArrayIndex columnCount(const MultiArrayView< 2, T, C > &x)
Definition: matrix.hxx:682
vigra::MultiArrayView
Base class for, and view to, vigra::MultiArray.
Definition: multi_array.hxx:652
vigra::SplineImageView
Create a continuous view onto a discrete image using splines.
Definition: splineimageview.hxx:99
vigra::rbfMatrix2DFromCorrespondingPoints
linalg::TemporaryMatrix< double > rbfMatrix2DFromCorrespondingPoints(SrcPointIterator s, SrcPointIterator s_end, DestPointIterator d, RadialBasisFunctor const &rbf)
Create a matrix that maps corresponding points onto each other using a given RBF. ...
Definition: rbf_registration.hxx:126
vigra::linalg::linearSolve
bool linearSolve(...)

© 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.11.0