Basic Image AlgorithmS Library  2.8.0
 All Classes Namespaces Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
TextureTransformHomography.cpp
1 /*
2 This file is part of the BIAS library (Basic ImageAlgorithmS).
3 
4 Copyright (C) 2003-2009 (see file CONTACT for details)
5  Multimediale Systeme der Informationsverarbeitung
6  Institut fuer Informatik
7  Christian-Albrechts-Universitaet Kiel
8 
9 
10 BIAS is free software; you can redistribute it and/or modify
11 it under the terms of the GNU Lesser General Public License as published by
12 the Free Software Foundation; either version 2.1 of the License, or
13 (at your option) any later version.
14 
15 BIAS is distributed in the hope that it will be useful,
16 but WITHOUT ANY WARRANTY; without even the implied warranty of
17 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18 GNU Lesser General Public License for more details.
19 
20 You should have received a copy of the GNU Lesser General Public License
21 along with BIAS; if not, write to the Free Software
22 Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
23 */
24 
25 #include <Base/Common/BIASpragmaStart.hh>
26 #include <Image/TextureTransformHomography.hh>
27 
28 
29 int BIAS::TextureTransformHomography::
30 MapForward(const HomgPoint2D& src, HomgPoint2D& sink) const {
31  sink = H_ * src;
32  if (fabs(sink[0])<1e-10) return -1;
33  sink.Homogenize();
34  return 0;
35 }
36 
37 int BIAS::TextureTransformHomography::
38 MapBackward(const HomgPoint2D& sink, HomgPoint2D& src) const {
39  src = Hinv_ * sink;
40  if (fabs(src[0])<1e-10) return -1;
41  src.Homogenize();
42  return 0;
43 }
44 
45 int BIAS::TextureTransformHomography::
46 ParameterJacobianForward(Matrix<double>& Jac, const HomgPoint2D& src) {
47  Jac.newsize(2,8);
48 
49  // get direct access to H entries
50  double* pH = H_.GetData();
51  const double c11 = *pH++;
52  const double c12 = *pH++;
53  const double c13 = *pH++;
54  const double c21 = *pH++;
55  const double c22 = *pH++;
56  const double c23 = *pH++;
57  const double c31 = *pH++;
58  const double c32 = *pH;
59 
60 
61  // precompute often used values
62  const double minusc1src = -1.0*(c11*src[0]+c12*src[1]+c13);
63  const double minusc2src = -1.0*(c21*src[0]+c22*src[1]+c23);
64  const double oneoverc3src = 1.0/(c31*src[0]+c32*src[1]+1.0);
65  const double oneoverc3src2 = oneoverc3src*oneoverc3src;
66 
67  // fill in jacobian
68  double* pData = Jac.GetData();
69  *pData++ = src[0] * oneoverc3src * AFFINE_SENSITIVITY;
70  *pData++ = src[1] * oneoverc3src * AFFINE_SENSITIVITY;
71  *pData++ = oneoverc3src;
72  *pData++ = 0;
73  *pData++ = 0;
74  *pData++ = 0;
75  *pData++ = minusc1src * oneoverc3src2 * src[0]* PROJECTIVE_SENSITIVITY;
76  *pData++ = minusc1src * oneoverc3src2 * src[1]* PROJECTIVE_SENSITIVITY;
77 
78  *pData++ = 0;
79  *pData++ = 0;
80  *pData++ = 0;
81  *pData++ = src[0] * oneoverc3src * AFFINE_SENSITIVITY;
82  *pData++ = src[1] * oneoverc3src * AFFINE_SENSITIVITY;
83  *pData++ = oneoverc3src;
84  *pData++ = minusc2src * oneoverc3src2 * src[0]* PROJECTIVE_SENSITIVITY;
85  *pData++ = minusc2src * oneoverc3src2 * src[1]* PROJECTIVE_SENSITIVITY;
86  return 0;
87 }
88 
89 
90 int BIAS::TextureTransformHomography::
91 ParameterJacobianBackward(Matrix<double>& Jac, const HomgPoint2D& sink) {
92  Jac.newsize(2,8);
93 
94  // get direct access to H entries
95  double* pH = H_.GetData();
96  const double c11 = *pH++ -1.0;
97  const double c12 = *pH++;
98  const double c13 = *pH++;
99  const double c21 = *pH++;
100  const double c22 = *pH++ - 1.0;
101  const double c23 = *pH++;
102  const double c31 = *pH++;
103  const double c32 = *pH;
104 
105  // precompute often used values
106  const double x = sink[0];
107  const double xx = x*x;
108  const double y = sink[1];
109  const double yy = y*y;
110  const double longterm = (x*c21*c32-x*c31*c22-x*c31-y*c32*c11-y*c32+y*c12*c31+c11*c22+c11+c22+1-c21*c12);
111  const double longtermlongterm =longterm*longterm;
112  const double anothertermx = x*c22+x-x*c23*c32-y*c12+y*c13*c32+c12*c23-c13*c22-c13;
113  const double anothertermy = x*c21-x*c23*c31-y*c11-y+y*c31*c13+c23*c11+c23-c13*c21;
114 // fill in jacobian
115  double* pData = Jac.GetData();
116 
117 
118  *pData++ = -(-y*c32+c22+1)*(anothertermx)/longtermlongterm* AFFINE_SENSITIVITY; ;
119  *pData++ = (-y*c11*c22+c23*c11*c22+yy*c11*c32-c23*c11*y*c32-y*c11+c23*c11+x*c22*c21-x*c23*c31*c22-y*c22+y*c31*c13*c22+c23*c22-c13*c21*c22+y*c31*c13+x*c23*c31*y*c32+yy*c32-x*c23*c31+c13*c21*y*c32-yy*c31*c13*c32-x*c21*y*c32-c13*c21+c23-y-c23*y*c32+x*c21)/longtermlongterm* AFFINE_SENSITIVITY; ;
120  *pData++ = -(-y*c32+c22+1)/(longterm);
121  *pData++ = (-x*c32+c12)*(anothertermx)/longtermlongterm* AFFINE_SENSITIVITY; ;
122  *pData++ = -(c12*c23*c11-c23*c11*x*c32+x*y*c32*c11-y*c12*c11-y*c12+x*y*c32+c13*x*c21*c32+c13*y*c12*c31+c12*c23+x*c21*c12-x*c12*c23*c31-y*c31*c13*x*c32-xx*c21*c32+xx*c23*c31*c32-c13*c21*c12-x*c23*c32)/longtermlongterm* AFFINE_SENSITIVITY; ;
123  *pData++ = (-x*c32+c12)/(longterm);
124 
125  *pData++ = -(y*c12-x*c22-x)*(anothertermx)/longtermlongterm * PROJECTIVE_SENSITIVITY;
126  *pData++ = (y*c12*x*c21-xx*c22*c21-xx*c21-x*c23*y*c12*c31+xx*c23*c31*c22+xx*c23*c31-yy*c11*c12+x*c22*y*c11+x*y*c11-yy*c12+x*c22*y+x*y+yy*c13*c12*c31-y*c31*c13*x*c22-y*c31*c13*x+c23*c11*y*c12-c23*c11*x*c22-c23*c11*x+c23*y*c12-c23*x*c22-x*c23-y*c13*c21*c12+c13*c22*x*c21+c13*x*c21)/longtermlongterm * PROJECTIVE_SENSITIVITY;
127 
128 
129 
130 
131  *pData++ = (y*c31*c13*c22+x*c22*c21-c13*c21*c22-y*x*c31*c22-y*c21*c12+c23*c21*c12+y*c31*c13+c13*c21*y*c32-yy*c31*c13*c32-c13*c21+yy*c12*c31+x*c23*c31*y*c32-y*x*c31-c23*x*c21*c32-c23*y*c12*c31+x*c21)/longtermlongterm* AFFINE_SENSITIVITY; ;
132 
133  *pData++ = -(c21-y*c31)*(anothertermy)/longtermlongterm* AFFINE_SENSITIVITY;
134  *pData++ = (c21-y*c31)/longterm;
135 
136 
137  *pData++ = -(x*c11*c22+x*c11-c23*c11*x*c32-y*c12*c11+c13*c11*y*c32+c12*c23*c11-c13*c11*c22-c13*c11+x*c22-xx*c22*c31-c13*c22+c13*c22*x*c31+x*y*c12*c31-y*c12-x*c12*c23*c31-xx*c31+x+c12*c23-y*c31*c13*x*c32+c13*x*c31-x*c23*c32+y*c13*c32+xx*c23*c31*c32-c13)/longtermlongterm* AFFINE_SENSITIVITY; ;
138  *pData++ = (c11+1-x*c31)*(anothertermy)/longtermlongterm* AFFINE_SENSITIVITY;
139 
140  *pData++ = -(c11+1-x*c31)/longterm;
141  *pData++ = (x*y*c11+c13*x*c21-xx*c21-yy*c12+c13*c22*x*c21-xx*c22*c21+x*c22*y-c12*c23*x*c21+x*c22*y*c11+y*c12*x*c21-y*c13*x*c21*c32-x*c23*y*c32*c11+yy*c13*c32*c11-y*c13*c11*c22+xx*c23*c21*c32-x*c23*y*c32+c23*c11*y*c12-y*c13+yy*c13*c32+c23*y*c12-yy*c11*c12-y*c13*c22-y*c13*c11+x*y)/longtermlongterm * PROJECTIVE_SENSITIVITY;
142  *pData++ = -(-x*c21+y*c11+y)*(anothertermy)/longtermlongterm * PROJECTIVE_SENSITIVITY;
143 
144  return 0;
145 }
146 
147 void BIAS::TextureTransformHomography::
148 SetParameters(const Vector<double>& p) {
149  P_ = p;
150  H_ = ConstructHFromParameters(p);
151  H_.GetInverse(Hinv_); // should check return value !!
152  Hinv_ /= Hinv_[2][2];
153 }
154 
155 void BIAS::TextureTransformHomography::
156 ComposeWithInverseDeltaP(const Vector<double>& deltaP) {
157  HMatrix H = ConstructHFromParameters(deltaP);
158  HMatrix HInv;
159  H.GetInverse(HInv);
160  H = HInv* H_;
161  Vector<double> newp = ExtractParameters(H);
162  SetParameters(newp);
163 }
164 
165 
166 int BIAS::TextureTransformHomography::
167 ParameterInversionJacobian(Matrix<double>& Jac) const {
168  Jac.newsize(8,8);
169  const double c11 = P_[0] * AFFINE_SENSITIVITY;
170  const double c12 = P_[1] * AFFINE_SENSITIVITY;
171  const double c13 = P_[2];
172  const double c21 = P_[3] * AFFINE_SENSITIVITY;
173  const double c22 = P_[4] * AFFINE_SENSITIVITY;
174  const double c23 = P_[5];
175  const double c31 = P_[6] * PROJECTIVE_SENSITIVITY;
176  const double c32 = P_[7] * PROJECTIVE_SENSITIVITY;
177 
178 
179  const double det = c11*c22+c11+c22+1-c21*c12;
180  const double detdet = det*det;
181 
182  double *pD = Jac.GetData();
183 
184 
185 
186  *pD++ = -(c22+1-c23*c32)/detdet*(c22+1);
187  *pD++ = (c22+1-c23*c32)/detdet*c21;
188  *pD++ = 0;
189  *pD++ = (c22+1-c23*c32)/detdet*c12;
190  *pD++ = (-c21*c12+c11*c23*c32+c23*c32)/detdet;
191  *pD++ = -c32/(det);
192  *pD++ = 0;
193  *pD++ = -c23/(det);
194 
195  *pD++ = (c12-c13*c32)/detdet*(c22+1);
196  *pD++ = -(c11*c22+c11+c22+1-c21*c13*c32)/detdet;
197  *pD++ = c32/(det);
198  *pD++ = -(c12-c13*c32)/detdet*c12;
199  *pD++ = (c12-c13*c32)/detdet*(c11+1);
200  *pD++ = 0;
201  *pD++ = 0;
202  *pD++ = c13/(det);
203 
204  *pD++ = -(c12*c23-c13*c22-c13)/detdet*(c22+1);
205  *pD++ = (c23*c11*c22+c23*c11+c23*c22+c23-c13*c21*c22-c13*c21)/detdet;
206  *pD++ = -(c22+1)/(det);
207  *pD++ = (c12*c23-c13*c22-c13)/detdet*c12;
208  *pD++ = -(c23*c11+c23-c13*c21)/detdet*c12;
209  *pD++ = c12/(det);
210  *pD++ = 0;
211  *pD++ = 0;
212 
213  *pD++ = (c21-c23*c31)/detdet*(c22+1);
214  *pD++ = -(c21-c23*c31)/detdet*c21;
215  *pD++ = 0;
216  *pD++ = -(c11*c22+c11+c22+1-c31*c12*c23)/detdet;
217  *pD++ = (c21-c23*c31)/detdet*(c11+1);
218  *pD++ = c31/(det);
219  *pD++ = c23/(det);
220  *pD++ = 0;
221 
222  *pD++ = -(c21*c12-c31*c13*c22-c31*c13)/detdet;
223  *pD++ = (c11+1-c31*c13)/detdet*c21;
224  *pD++ = -c31/(det);
225  *pD++ = (c11+1-c31*c13)/detdet*c12;
226  *pD++ = -(c11+1-c31*c13)/detdet*(c11+1);
227  *pD++ = 0;
228  *pD++ = -c13/(det);
229  *pD++ = 0;
230 
231  *pD++ = (c12*c23-c13*c22-c13)/detdet*c21;
232  *pD++ = -(c23*c11+c23-c13*c21)/detdet*c21;
233  *pD++ = c21/(det);
234  *pD++ = -(-c13*c11*c22-c13*c11-c13*c22-c13+c12*c23*c11+c12*c23)/detdet;
235  *pD++ = (c23*c11+c23-c13*c21)/detdet*(c11+1);
236  *pD++ = -(c11+1)/(det);
237  *pD++ = 0;
238  *pD++ = 0;
239 
240  *pD++ = -(c21*c32-c31*c22-c31)/detdet*(c22+1);
241  *pD++ = (c21*c32-c31*c22-c31)/detdet*c21;
242  *pD++ = 0;
243  *pD++ = (c32*c11*c22+c32*c11+c32*c22+c32-c12*c31*c22-c12*c31)/detdet;
244  *pD++ = -(c32*c11+c32-c12*c31)/detdet*c21;
245  *pD++ = 0;
246  *pD++ = -(c22+1)/(det);
247  *pD++ = c21/(det);
248 
249  *pD++ = (c21*c32-c31*c22-c31)/detdet*c12;
250  *pD++ = -(-c31*c11*c22-c31*c11-c31*c22-c31+c21*c32*c11+c21*c32)/detdet;
251  *pD++ = 0;
252  *pD++ = -(c32*c11+c32-c12*c31)/detdet*c12;
253  *pD++ = (c32*c11+c32-c12*c31)/detdet*(c11+1);
254  *pD++ = 0;
255  *pD++ = c12/(det);
256  *pD++ = -(c11+1)/(det);
257 
258 
259 
260 
261  for (unsigned int i=0; i<8; i++) {
262  double factor = 1.0;
263  switch (i) {
264  case 0:
265  case 1:
266  case 3:
267  case 4: factor = AFFINE_SENSITIVITY; break;
268  case 6:
269  case 7: factor = PROJECTIVE_SENSITIVITY; break;
270  }
271  for (unsigned int j=0; j<8; j++) {
272  switch (j) {
273  case 0:
274  case 1:
275  case 3:
276  case 4: Jac[j][i] *= factor / AFFINE_SENSITIVITY; break;
277  case 6:
278  case 7: Jac[j][i] *= factor / PROJECTIVE_SENSITIVITY; break;
279  }
280  }
281  }
282 
283 
284  return 0;
285 }
286 
287 
288 BIAS::Vector<double> BIAS::TextureTransformHomography::
289 ExtractParameters(const BIAS::HMatrix& H) const{
290  Vector<double> p(8);
291  const double norminv = 1.0/H[2][2];
292  for (unsigned int i=0; i<8; i++) p[i] = H[i/3][i%3]*norminv;
293 
294  // zero warp must be identity:
295  p[0] -= 1.0;
296  p[4] -= 1.0;
297 
298  p[0] /= AFFINE_SENSITIVITY;
299  p[1] /= AFFINE_SENSITIVITY;
300  p[3] /= AFFINE_SENSITIVITY;
301  p[4] /= AFFINE_SENSITIVITY;
302 
303  p[6] /= PROJECTIVE_SENSITIVITY;
304  p[7] /= PROJECTIVE_SENSITIVITY;
305 
306  return p;
307 }
308 
309 
310 BIAS::HMatrix BIAS::TextureTransformHomography::
311 ConstructHFromParameters(const Vector<double>&p) const {
312  HMatrix H;
313  H[2][2] = 1.0;
314 
315  H[0][0] = p[0] * AFFINE_SENSITIVITY + 1.0;
316  H[0][1] = p[1] * AFFINE_SENSITIVITY;
317  H[0][2] = p[2];
318  H[1][0] = p[3] * AFFINE_SENSITIVITY;
319  H[1][1] = p[4] * AFFINE_SENSITIVITY + 1.0;
320  H[1][2] = p[5];
321  H[2][0] = p[6] * PROJECTIVE_SENSITIVITY;
322  H[2][1] = p[7] * PROJECTIVE_SENSITIVITY;
323 
324  return H;
325 }
BIAS::TextureTransformHomography::MapBackward
int MapBackward(const HomgPoint2D &sink, HomgPoint2D &src) const
map a point in image &quot;source&quot; to a point in image &quot;sink&quot;
Definition: TextureTransformHomography.cpp:38
BIAS::TextureTransformHomography::ParameterInversionJacobian
int ParameterInversionJacobian(Matrix< double > &Jac) const
compute parameters for inverse operation and obtain the jacobian of the inverse parameters with respe...
Definition: TextureTransformHomography.cpp:167
BIAS::HomgPoint2D
class HomgPoint2D describes a point with 2 degrees of freedom in projective coordinates.
Definition: HomgPoint2D.hh:67
BIAS::HMatrix
a 3x3 Matrix describing projective transformations between planes
Definition: HMatrix.hh:39
BIAS::TextureTransformHomography::MapForward
int MapForward(const HomgPoint2D &src, HomgPoint2D &sink) const
map a point in image &quot;source&quot; to a point in image &quot;sink&quot;
Definition: TextureTransformHomography.cpp:30
BIAS::HomgPoint2D::Homogenize
HomgPoint2D & Homogenize()
homogenize class data member elements to W==1 by divison by W
Definition: HomgPoint2D.hh:215
TNT::Matrix::newsize
Matrix< T > & newsize(Subscript M, Subscript N)
Definition: cmat.h:269
BIAS::TextureTransformHomography::SetParameters
void SetParameters(const Vector< double > &p)
h11-1, h12, h13, h21, h22-1, h23, h31, h32 (h33==1 assumed)
Definition: TextureTransformHomography.cpp:148
BIAS::Matrix3x3::GetInverse
int GetInverse(Matrix3x3< T > &inv) const
Matrix inversion: inverts this and stores resulty in argument inv.
Definition: Matrix3x3.cpp:373
BIAS::TextureTransformHomography::ConstructHFromParameters
BIAS::HMatrix ConstructHFromParameters(const Vector< double > &p) const
given parameters, compute H
Definition: TextureTransformHomography.cpp:311
BIAS::Matrix::GetData
T * GetData()
get the pointer to the data array of the matrix (for faster direct memeory access) ...
Definition: Matrix.hh:185
BIAS::TextureTransformHomography::ExtractParameters
Vector< double > ExtractParameters(const BIAS::HMatrix &H) const
given H, compute a parameter vector
Definition: TextureTransformHomography.cpp:289
BIAS::TextureTransformHomography::ComposeWithInverseDeltaP
void ComposeWithInverseDeltaP(const Vector< double > &deltaP)
concatenate *this and an inverse transform with param deltaP and save new parameter vector to *this...
Definition: TextureTransformHomography.cpp:156
BIAS::TextureTransformHomography::ParameterJacobianBackward
int ParameterJacobianBackward(Matrix< double > &Jac, const HomgPoint2D &src)
transformed position change when parameters change
Definition: TextureTransformHomography.cpp:91
BIAS::TextureTransformHomography::ParameterJacobianForward
int ParameterJacobianForward(Matrix< double > &Jac, const HomgPoint2D &src)
transformed position change when parameters change
Definition: TextureTransformHomography.cpp:46
Generated on Tue Feb 6 2018 02:01:40 for Basic Image AlgorithmS Library by   doxygen 1.8.5