1 /****************************************************************************
2 *
3 * ViSP, open source Visual Servoing Platform software.
4 * Copyright (C) 2005 - 2019 by Inria. All rights reserved.
5 *
6 * This software is free software; you can redistribute it and/or modify
7 * it under the terms of the GNU General Public License as published by
8 * the Free Software Foundation; either version 2 of the License, or
9 * (at your option) any later version.
10 * See the file LICENSE.txt at the root directory of this source
11 * distribution for additional information about the GNU GPL.
12 *
13 * For using ViSP with software that can not be combined with the GNU
14 * GPL, please contact Inria about acquiring a ViSP Professional
15 * Edition License.
16 *
17 * See http://visp.inria.fr for more information.
18 *
19 * This software was developed at:
20 * Inria Rennes - Bretagne Atlantique
21 * Campus Universitaire de Beaulieu
22 * 35042 Rennes Cedex
23 * France
24 *
25 * If you have questions regarding the use of this file, please contact
26 * Inria at visp@inria.fr
27 *
28 * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
29 * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
30 *
31 * Description:
32 * Template tracker.
33 *
34 * Authors:
35 * Amaury Dame
36 * Aurelien Yol
37 * Fabien Spindler
38 *
39 *****************************************************************************/
40 #include <visp3/core/vpTrackingException.h>
41 #include <visp3/tt/vpTemplateTrackerWarpHomography.h>
42
43 /*!
44 * Construct an homography model with 8 parameters initialized to zero.
45 */
vpTemplateTrackerWarpHomography()46 vpTemplateTrackerWarpHomography::vpTemplateTrackerWarpHomography()
47 {
48 nbParam = 8;
49 }
50
51 /*!
52 * Get the parameters of the warping function one level down
53 * where image size is divided by two along the lines and the columns.
54 * \param p : 8-dim vector that contains the current parameters of the warping function.
55 * \param p_down : 8-dim vector that contains the resulting parameters one level down.
56 */
getParamPyramidDown(const vpColVector & p,vpColVector & p_down)57 void vpTemplateTrackerWarpHomography::getParamPyramidDown(const vpColVector &p, vpColVector &p_down)
58 {
59 p_down[0] = p[0];
60 p_down[1] = p[1];
61 p_down[2] = p[2] * 2.;
62 p_down[3] = p[3];
63 p_down[4] = p[4];
64 p_down[5] = p[5] * 2.;
65 p_down[6] = p[6] / 2.;
66 p_down[7] = p[7] / 2.;
67 }
68
69 /*!
70 * Get the parameters of the warping function one level up
71 * where image size is multiplied by two along the lines and the columns.
72 * \param p : 8-dim vector that contains the current parameters of the warping function.
73 * \param p_up : 8-dim vector that contains the resulting parameters one level up.
74 */
getParamPyramidUp(const vpColVector & p,vpColVector & p_up)75 void vpTemplateTrackerWarpHomography::getParamPyramidUp(const vpColVector &p, vpColVector &p_up)
76 {
77 p_up[0] = p[0];
78 p_up[1] = p[1];
79 p_up[2] = p[2] / 2.;
80 p_up[3] = p[3];
81 p_up[4] = p[4];
82 p_up[5] = p[5] / 2.;
83 p_up[6] = p[6] * 2.;
84 p_up[7] = p[7] * 2.;
85 }
86
87 /*!
88 * Compute the derivative of the image with relation to the warping function parameters.
89 * \param v : Coordinate (along the image rows axis) of the point to consider in the image.
90 * \param u : Coordinate (along the image columns axis) of the point to consider in the image.
91 * \param dv : Derivative on the v-axis (along the rows) of the point (u,v).
92 * \param du : Derivative on the u-axis (along the columns) of the point (u,v).
93 * \param dIdW : Resulting derivative matrix (image according to the warping function).
94 */
getdW0(const int & v,const int & u,const double & dv,const double & du,double * dIdW)95 void vpTemplateTrackerWarpHomography::getdW0(const int &v, const int &u, const double &dv, const double &du, double *dIdW)
96 {
97 double u_du_ = u * du;
98 double v_dv_ = v * dv;
99 dIdW[0] = u_du_;
100 dIdW[1] = u * dv;
101 dIdW[2] = -u * (u_du_ + v_dv_);
102 dIdW[3] = v * du;
103 dIdW[4] = v_dv_;
104 dIdW[5] = -v * (u_du_ + v_dv_);
105 dIdW[6] = du;
106 dIdW[7] = dv;
107 }
108
109 /*!
110 * Compute the derivative of the warping model \f$M\f$ according to the initial parameters \f$p_0\f$
111 * at point \f$X=(u,v)\f$:
112 * \f[
113 * \frac{\partial M}{\partial p}(X, p_0)
114 * \f]
115 *
116 * \param v : Coordinate (along the image rows axis) of the point X(u,v) to consider in the image.
117 * \param u : Coordinate (along the image columns axis) of the point X(u,v) to consider in the image.
118 * \param dIdW : Resulting 2-by-8 derivative matrix.
119 */
getdWdp0(const int & v,const int & u,double * dIdW)120 void vpTemplateTrackerWarpHomography::getdWdp0(const int &v, const int &u, double *dIdW)
121 {
122 double uv_ = u * v;
123 dIdW[0] = u;
124 dIdW[1] = 0;
125 dIdW[2] = -u * u;
126 dIdW[3] = v;
127 dIdW[4] = 0;
128 dIdW[5] = - uv_;
129 dIdW[6] = 1.;
130 dIdW[7] = 0;
131
132 dIdW[8] = 0;
133 dIdW[9] = u;
134 dIdW[10] = - uv_;
135 dIdW[11] = 0;
136 dIdW[12] = v;
137 dIdW[13] = -v * v;
138 dIdW[14] = 0;
139 dIdW[15] = 1.;
140 }
141
142 /*!
143 * Compute the projection denominator (Z) used in x = X/Z and y = Y/Z.
144 * \param X : Point with coordinates (u, v) to consider.
145 * \param p : 8-dim vector that contains the current parameters of the warping function.
146 *
147 * \sa warpX(const vpColVector &, vpColVector &, const vpColVector &)
148 * \sa warpX(const int &, const int &, double &, double &, const vpColVector &)
149 * \sa dWarp(), dWarpCompo()
150 */
computeDenom(vpColVector & X,const vpColVector & p)151 void vpTemplateTrackerWarpHomography::computeDenom(vpColVector &X, const vpColVector &p)
152 {
153 double value = (p[2] * X[0] + p[5] * X[1] + 1.);
154
155 if (std::fabs(value) > std::numeric_limits<double>::epsilon()) {
156 denom = (1. / value);
157 } else {
158 throw(vpTrackingException(vpTrackingException::fatalError,
159 "Division by zero in vpTemplateTrackerWarpHomography::computeDenom()"));
160 }
161 }
162
163 /*!
164 * Warp point \f$X_1=(u_1,v_1)\f$ using the transformation model with parameters \f$p\f$.
165 * \f[X_2 = {^2}M_1(p) * X_1\f]
166 * \param v1 : Coordinate (along the image rows axis) of the point \f$X_1=(u_1,v_1)\f$ to warp.
167 * \param u1 : Coordinate (along the image columns axis) of the point \f$X_1=(u_1,v_1)\f$ to warp.
168 * \param v2 : Coordinate of the warped point \f$X_2=(u_2,v_2)\f$ along the image rows axis.
169 * \param u2 : Coordinate of the warped point \f$X_2=(u_2,v_2)\f$ along the image column axis.
170 * \param p : 8-dim vector that contains the parameters of the transformation.
171 */
warpX(const int & v1,const int & u1,double & v2,double & u2,const vpColVector & p)172 void vpTemplateTrackerWarpHomography::warpX(const int &v1, const int &u1, double &v2, double &u2, const vpColVector &p)
173 {
174 u2 = ((1. + p[0]) * u1 + p[3] * v1 + p[6]) * denom;
175 v2 = (p[1] * u1 + (1. + p[4]) * v1 + p[7]) * denom;
176 }
177
178 /*!
179 * Warp point \f$X_1=(u_1,v_1)\f$ using the transformation model.
180 * \f[X_2 = {^2}M_1(p) * X_1\f]
181 * \param X1 : 2-dim vector corresponding to the coordinates \f$(u_1, v_1)\f$ of the point to warp.
182 * \param X2 : 2-dim vector corresponding to the coordinates \f$(u_2, v_2)\f$ of the warped point.
183 * \param p : 8-dim vector that contains the parameters of the transformation.
184 *
185 * \sa computeDenom()
186 */
warpX(const vpColVector & X1,vpColVector & X2,const vpColVector & p)187 void vpTemplateTrackerWarpHomography::warpX(const vpColVector &X1, vpColVector &X2, const vpColVector &p)
188 {
189 X2[0] = ((1 + p[0]) * X1[0] + p[3] * X1[1] + p[6]) * denom;
190 X2[1] = (p[1] * X1[0] + (1 + p[4]) * X1[1] + p[7]) * denom;
191 }
192
193 /*!
194 * Compute the derivative matrix of the warping function at point \f$X=(u,v)\f$ according to the model parameters:
195 * \f[
196 * \frac{\partial M}{\partial p}(X, p)
197 * \f]
198 * \param X : 2-dim vector corresponding to the coordinates \f$(u_1, v_1)\f$ of the point to
199 * consider in the derivative computation.
200 * \param dM : Resulting warping model derivative returned as a 2-by-8 matrix.
201 *
202 * \sa computeDenom()
203 */
dWarp(const vpColVector & X,const vpColVector &,const vpColVector &,vpMatrix & dM)204 void vpTemplateTrackerWarpHomography::dWarp(const vpColVector &X, const vpColVector &, const vpColVector &,
205 vpMatrix &dM)
206 {
207 double u = X[0];
208 double v = X[1];
209 dM[0][0] = u * denom;
210 dM[0][1] = 0;
211 dM[0][2] = -u * X[0] * denom;
212 dM[0][3] = v * denom;
213 dM[0][4] = 0;
214 dM[0][5] = -v * X[0] * denom;
215 dM[0][6] = denom;
216
217 dM[1][1] = u * denom;
218 dM[1][2] = -u * X[1] * denom;
219 dM[1][4] = v * denom;
220 dM[1][5] = -v * X[1] * denom;
221 dM[1][7] = denom;
222 }
223
224 /*!
225 * Compute the compositionnal derivative matrix of the warping function according to the model parameters.
226 * \param X : 2-dim vector corresponding to the coordinates \f$(u_1, v_1)\f$ of the point to
227 * consider in the derivative computation.
228 * \param p : 8-dim vector that contains the parameters of the warping function.
229 * \param dwdp0 : 2-by-8 derivative matrix of the warping function according to
230 * the initial warping function parameters (p=0).
231 * \param dM : Resulting warping model compositionnal derivative returned as a 2-by-8 matrix.
232 *
233 * \sa computeDenom()
234 */
dWarpCompo(const vpColVector &,const vpColVector & X,const vpColVector & p,const double * dwdp0,vpMatrix & dM)235 void vpTemplateTrackerWarpHomography::dWarpCompo(const vpColVector &, const vpColVector &X,
236 const vpColVector &p, const double *dwdp0, vpMatrix &dM)
237 {
238 double dwdx0, dwdx1;
239 double dwdy0, dwdy1;
240
241 dwdx0 = ((1. + p[0]) - X[0] * p[2]) * denom;
242 dwdx1 = (p[1] - X[1] * p[2]) * denom;
243 dwdy0 = (p[3] - X[0] * p[5]) * denom;
244 dwdy1 = ((1. + p[4]) - X[1] * p[5]) * denom;
245 for (unsigned int i = 0; i < nbParam; i++) {
246 dM[0][i] = dwdx0 * dwdp0[i] + dwdy0 * dwdp0[i + nbParam];
247 dM[1][i] = dwdx1 * dwdp0[i] + dwdy1 * dwdp0[i + nbParam];
248 }
249 }
250
251 /*!
252 * Warp a point X1 with the inverse transformation \f$M\f$.
253 * \f[ X_2 = {\left({^1}M_2\right)}^{-1} \; X_1\f]
254 * \param X1 : 2-dim vector corresponding to the coordinates (u,v) of the point to warp.
255 * \param X2 : 2-dim vector corresponding to the coordinates (u,v) of the warped point.
256 * \param p : Parameters corresponding to the warping model \f${^1}M_2\f$.
257 */
warpXInv(const vpColVector & X1,vpColVector & X2,const vpColVector & p)258 void vpTemplateTrackerWarpHomography::warpXInv(const vpColVector &X1, vpColVector &X2, const vpColVector &p)
259 {
260 double value = (p[2] * X1[0] + p[5] * X1[1] + 1.);
261
262 if (std::fabs(value) > std::numeric_limits<double>::epsilon()) {
263 X2[0] = ((1 + p[0]) * X1[0] + p[3] * X1[1] + p[6]) / value;
264 X2[1] = (p[1] * X1[0] + (1 + p[4]) * X1[1] + p[7]) / value;
265 } else {
266 throw(vpTrackingException(vpTrackingException::fatalError, "Division by zero in "
267 "vpTemplateTrackerWarpHomography::"
268 "warpXInv()"));
269 }
270 }
271
272 /*!
273 * Compute inverse of the warping transformation.
274 * \param p : 8-dim vector that contains the parameters corresponding
275 * to the transformation to inverse.
276 * \param p_inv : 8-dim vector that contains the parameters of the inverse transformation \f$ {M(p)}^{-1}\f$.
277 */
getParamInverse(const vpColVector & p,vpColVector & p_inv) const278 void vpTemplateTrackerWarpHomography::getParamInverse(const vpColVector &p, vpColVector &p_inv) const
279 {
280 double h_00 = 1. + p[0];
281 double h_10 = p[1];
282 double h_20 = p[2];
283 double h_01 = p[3];
284 double h_11 = 1. + p[4];
285 double h_21 = p[5];
286 double h_02 = p[6];
287 double h_12 = p[7];
288
289 double h_inv_22 = (h_00 * h_11 - h_01 * h_10);
290
291 if (std::fabs(h_inv_22) < std::numeric_limits<double>::epsilon()) {
292 throw(vpException(vpException::fatalError, "Cannot get homography inverse parameters. Matrix determinant is 0."));
293 }
294
295 p_inv[0] = (h_11 - h_12 * h_21) / h_inv_22 - 1.;
296 p_inv[3] = (h_02 * h_21 - h_01) / h_inv_22;
297 p_inv[6] = (h_01 * h_12 - h_02 * h_11) / h_inv_22;
298
299 p_inv[1] = (h_12 * h_20 - h_10) / h_inv_22;
300 p_inv[4] = (h_00 - h_02 * h_20) / h_inv_22 - 1.;
301 p_inv[7] = (h_02 * h_10 - h_00 * h_12) / h_inv_22;
302
303 p_inv[2] = (h_10 * h_21 - h_11 * h_20) / h_inv_22;
304 p_inv[5] = (h_01 * h_20 - h_00 * h_21) / h_inv_22;
305 }
306
307 /*!
308 * Return the homography corresponding to the parameters.
309 * \param p : 8-dim vector that contains the parameters corresponding
310 * to the transformation to inverse.
311 * \return Corresponding homography.
312 */
getHomography(const vpColVector & p) const313 vpHomography vpTemplateTrackerWarpHomography::getHomography(const vpColVector &p) const
314 {
315 vpHomography H;
316 H[0][0] = 1. + p[0];
317 H[1][0] = p[1];
318 H[2][0] = p[2];
319 H[0][1] = p[3];
320 H[1][1] = 1. + p[4];
321 H[2][1] = p[5];
322 H[0][2] = p[6];
323 H[1][2] = p[7];
324 H[2][2] = 1.;
325
326 return H;
327 }
328
329 /*!
330 * Compute the parameters corresponding to an homography.
331 * \param[in] H : Homography
332 * \param[out] p : 8-dim vector corresponding to the homography.
333 */
getParam(const vpHomography & H,vpColVector & p) const334 void vpTemplateTrackerWarpHomography::getParam(const vpHomography &H, vpColVector &p) const
335 {
336 p.resize(getNbParam(), false);
337 p[0] = H[0][0] / H[2][2] - 1.;
338 p[1] = H[1][0] / H[2][2];
339 p[2] = H[2][0] / H[2][2];
340 p[3] = H[0][1] / H[2][2];
341 p[4] = H[1][1] / H[2][2] - 1.;
342 p[5] = H[2][1] / H[2][2];
343 p[6] = H[0][2] / H[2][2];
344 p[7] = H[1][2] / H[2][2];
345 }
346
347 /*!
348 * Compute the parameters corresponding to an homography as a 3-by-3 matrix.
349 * \param[in] H : 3-by-3 matrix corresponding to an homography
350 * \param[out] p : 8-dim vector corresponding to the homography.
351 */
getParam(const vpMatrix & H,vpColVector & p) const352 void vpTemplateTrackerWarpHomography::getParam(const vpMatrix &H, vpColVector &p) const
353 {
354 p.resize(getNbParam(), false);
355 p[0] = H[0][0] / H[2][2] - 1.;
356 p[1] = H[1][0] / H[2][2];
357 p[2] = H[2][0] / H[2][2];
358 p[3] = H[0][1] / H[2][2];
359 p[4] = H[1][1] / H[2][2] - 1.;
360 p[5] = H[2][1] / H[2][2];
361 p[6] = H[0][2] / H[2][2];
362 p[7] = H[1][2] / H[2][2];
363 }
364
365 /*!
366 * Compute the transformation resulting from the composition of two other transformations.
367 * \param p1 : 8-dim vector that contains the parameters corresponding
368 * to first transformation.
369 * \param p2 : 8-dim vector that contains the parameters corresponding
370 * to second transformation.
371 * \param p12 : 8-dim vector that contains the resulting transformation \f$ p_{12} = p_1 \circ p_2\f$.
372 */
pRondp(const vpColVector & p1,const vpColVector & p2,vpColVector & p12) const373 void vpTemplateTrackerWarpHomography::pRondp(const vpColVector &p1, const vpColVector &p2, vpColVector &p12) const
374 {
375 double h1_00 = 1. + p1[0];
376 double h1_10 = p1[1];
377 double h1_20 = p1[2];
378 double h1_01 = p1[3];
379 double h1_11 = 1. + p1[4];
380 double h1_21 = p1[5];
381 double h1_02 = p1[6];
382 double h1_12 = p1[7];
383
384 double h2_00 = 1. + p2[0];
385 double h2_10 = p2[1];
386 double h2_20 = p2[2];
387 double h2_01 = p2[3];
388 double h2_11 = 1. + p2[4];
389 double h2_21 = p2[5];
390 double h2_02 = p2[6];
391 double h2_12 = p2[7];
392
393 double h12_22 = h1_20 * h2_02 + h1_21 * h2_12 + 1.;
394
395 p12[0] = (h1_00 * h2_00 + h1_01 * h2_10 + h1_02 * h2_20) / h12_22 - 1.;
396 p12[3] = (h1_00 * h2_01 + h1_01 * h2_11 + h1_02 * h2_21) / h12_22;
397 p12[6] = (h1_00 * h2_02 + h1_01 * h2_12 + h1_02) / h12_22;
398
399 p12[1] = (h1_10 * h2_00 + h1_11 * h2_10 + h1_12 * h2_20) / h12_22;
400 p12[4] = (h1_10 * h2_01 + h1_11 * h2_11 + h1_12 * h2_21) / h12_22 - 1.;
401 p12[7] = (h1_10 * h2_02 + h1_11 * h2_12 + h1_12) / h12_22;
402
403 p12[2] = (h1_20 * h2_00 + h1_21 * h2_10 + h2_20) / h12_22;
404 p12[5] = (h1_20 * h2_01 + h1_21 * h2_11 + h2_21) / h12_22;
405 }
406