[ 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
46namespace vigra {
47
48/** \addtogroup Registration Image Registration
49 */
50//@{
51
52namespace detail {
53
54// Small and hopefully fast helper function to compute the squared-distance
55// between two points (2d)
56template <class SrcPoint, class DestPoint>
57inline 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 */
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 */
90template<int N>
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 */
123template <class RadialBasisFunctor,
125linalg::TemporaryMatrix<double>
127{
128 int point_count = s_end - s;
129
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 */
232template <int ORDER, class T,
233 class DestIterator, class DestAccessor,
234 class DestPointIterator,
235 class C,
236 class RadialBasisFunctor>
237void
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
276template <int ORDER, class T,
277 class DestIterator, class DestAccessor,
278 class DestPointIterator,
279 class C,
280 class RadialBasisFunctor>
281inline
282void 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
291template <int ORDER, class T,
292 class T2, class S2,
293 class DestPointIterator,
294 class C,
295 class RadialBasisFunctor>
296inline
297void 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 */
Two dimensional difference vector.
Definition diff2d.hxx:186
Class for a single RGB value.
Definition rgbvalue.hxx:128
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
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
Definition rbf_registration.hxx:92
Definition rbf_registration.hxx:70

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