Basic Image AlgorithmS Library  2.8.0
 All Classes Namespaces Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
RegionMatcher.hh
1 /* This file is part of the BIAS library (Basic ImageAlgorithmS).
2 
3  Copyright (C) 2003-2009 (see file CONTACT for details)
4  Multimediale Systeme der Informationsverarbeitung
5  Institut fuer Informatik
6  Christian-Albrechts-Universitaet Kiel
7 
8  BIAS is free software; you can redistribute it and/or modify
9  it under the terms of the GNU Lesser General Public License as published by
10  the Free Software Foundation; either version 2.1 of the License, or
11  (at your option) any later version.
12 
13  BIAS is distributed in the hope that it will be useful,
14  but WITHOUT ANY WARRANTY; without even the implied warranty of
15  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16  GNU Lesser General Public License for more details.
17 
18  You should have received a copy of the GNU Lesser General Public License
19  along with BIAS; if not, write to the Free Software
20  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA*/
21 #ifndef __RegionMatcher_hh__
22 #define __RegionMatcher_hh__
23 
24 #include <bias_config.h>
25 
26 #include <Base/Debug/Debug.hh>
27 #include <Base/Debug/Error.hh>
28 #include <MathAlgo/Lapack.hh>
29 #include <Base/Math/Matrix2x2.hh>
30 #include <Base/Math/Vector2.hh>
31 
32 #define D_REGION_MATCHER_SEARCHNCC 0x0001
33 #define D_REGION_MATCHER_KLT 0x0002
34 #define D_REGION_MATCHER_LINE 0x0004
35 #define D_REGION_MATCHER_PARABOLA 0x0008
36 #define D_REGION_MATCHER_COLORNCC 0x0010
37 #define D_REGION_MATCHER_COLORCCV 0x0020
38 #define D_REGION_MATCHER_KLT_EDGE 0x0040
39 #define D_REGION_MATCHER_SEARCH_EVEN 0x0080
40 #define D_REGION_MATCHER_SEARCH_ODD 0x0100
41 
42 // data type for internal calculations in KLT
43 #define KLT_TYPE double
44 //#define KLT_TYPE float
45 
46 namespace BIAS
47 {
48 
49  // forward declaration
50  template<class T>
51  class BIASImageBase_EXPORT Image;
52 
53  /** some fast SAD functions via defines
54  @author woelk */
55 #define ABSDIFF(a, b) ((a>b)?(a-b):(b-a))
56 #define SAD3x3(x1, y1, x2, y2, ida1, ida2, re)\
57  unsigned x1m=x1-1, x1p=x1+1, y1m=y1-1, y1p=y1+1;\
58  unsigned x2m=x2-1, x2p=x2+1, y2m=y2-1, y2p=y2+1;\
59  re = (unsigned)ABSDIFF(ida1[y1m][x1m], ida2[y2m][x2m]) + \
60  (unsigned)ABSDIFF(ida1[y1m][x1], ida2[y2m][x2]) + \
61  (unsigned)ABSDIFF(ida1[y1m][x1p], ida2[y2m][x2p]) + \
62  (unsigned)ABSDIFF(ida1[y1][x1m], ida2[y2][x2m]) + \
63  (unsigned)ABSDIFF(ida1[y1][x1], ida2[y2][x2]) + \
64  (unsigned)ABSDIFF(ida1[y1][x1p], ida2[y2][x2p]) + \
65  (unsigned)ABSDIFF(ida1[y1p][x1m], ida2[y2p][x2m]) + \
66  (unsigned)ABSDIFF(ida1[y1p][x1], ida2[y2p][x2]) + \
67  (unsigned)ABSDIFF(ida1[y1p][x1p], ida2[y2p][x2p]);
68 
69 #define SAD5x5(x1, y1, x2, y2, ida1, ida2, re)\
70  unsigned x1m=x1-1, x1p=x1+1, y1m=y1-1, y1p=y1+1, x1m2=x1-2, x1p2=x2+2, y1m2=y2-2, y1p2=y2+2;\
71  unsigned x2m=x2-1, x2p=x2+1, y2m=y2-1, y2p=y2+1, x2m2=x1-2, x2p2=x2+2, y2m2=y2-2, y2p2=y2+2;\
72  re = (unsigned)ABSDIFF(ida1[y1m2][x1m2], ida2[y2m2][x2m2]) + \
73  (unsigned)ABSDIFF(ida1[y1m2][x1m], ida2[y2m2][x2m]) + \
74  (unsigned)ABSDIFF(ida1[y1m2][x1], ida2[y2m2][x2]) + \
75  (unsigned)ABSDIFF(ida1[y1m2][x1p], ida2[y2m2][x2p]) + \
76  (unsigned)ABSDIFF(ida1[y1m2][x1p2], ida2[y2m2][x2p2]) + \
77  (unsigned)ABSDIFF(ida1[y1m][x1m2], ida2[y2m][x2m2]) + \
78  (unsigned)ABSDIFF(ida1[y1m][x1m], ida2[y2m][x2m]) + \
79  (unsigned)ABSDIFF(ida1[y1m][x1], ida2[y2m][x2]) + \
80  (unsigned)ABSDIFF(ida1[y1m][x1p], ida2[y2m][x2p]) + \
81  (unsigned)ABSDIFF(ida1[y1m][x1p2], ida2[y2m][x2p2]) + \
82  (unsigned)ABSDIFF(ida1[y1][x1m2], ida2[y2][x2m2]) + \
83  (unsigned)ABSDIFF(ida1[y1][x1m], ida2[y2][x2m]) + \
84  (unsigned)ABSDIFF(ida1[y1][x1], ida2[y2][x2]) + \
85  (unsigned)ABSDIFF(ida1[y1][x1p], ida2[y2][x2p]) + \
86  (unsigned)ABSDIFF(ida1[y1][x1p2], ida2[y2][x2p2]) + \
87  (unsigned)ABSDIFF(ida1[y1p][x1m2], ida2[y2p][x2m2]) + \
88  (unsigned)ABSDIFF(ida1[y1p][x1m], ida2[y2p][x2m]) + \
89  (unsigned)ABSDIFF(ida1[y1p][x1], ida2[y2p][x2]) + \
90  (unsigned)ABSDIFF(ida1[y1p][x1p], ida2[y2p][x2p]) + \
91  (unsigned)ABSDIFF(ida1[y1p][x1p2], ida2[y2p][x2p2]) + \
92  (unsigned)ABSDIFF(ida1[y1p2][x1m2], ida2[y2p2][x2m2]) + \
93  (unsigned)ABSDIFF(ida1[y1p2][x1m], ida2[y2p2][x2m]) + \
94  (unsigned)ABSDIFF(ida1[y1p2][x1], ida2[y2p2][x2]) + \
95  (unsigned)ABSDIFF(ida1[y1p2][x1p], ida2[y2p2][x2p]) + \
96  (unsigned)ABSDIFF(ida1[y1p2][x1p2], ida2[y2p2][x2p2]);
97 
98 #define NCC5x5(x1, y1, x2, y2, ida1, ida2, re)\
99  unsigned x1m=x1-1, x1p=x1+1, y1m=y1-1, y1p=y1+1, x1m2=x1-2, x1p2=x2+2, y1m2=y2-2, y1p2=y2+2;\
100  unsigned x2m=x2-1, x2p=x2+1, y2m=y2-1, y2p=y2+1, x2m2=x1-2, x2p2=x2+2, y2m2=y2-2, y2p2=y2+2;\
101  int im1, im2; unsigned long S11=0, S12=0, S22=0, N=0, S1=0, S2=0;\
102  im1=ida1[y1m2][x1m2]; im2=ida2[y2m2][x2m2]; N++; S1+=im1; S2+=im2; S11+=im1*im1; S12+=im1*im2; S22+=im2*im2; \
103  im1=ida1[y1m2][x1m]; im2=ida2[y2m2][x2m]; N++; S1+=im1; S2+=im2; S11+=im1*im1; S12+=im1*im2; S22+=im2*im2; \
104  im1=ida1[y1m2][x1]; im2=ida2[y2m2][x2]; N++; S1+=im1; S2+=im2; S11+=im1*im1; S12+=im1*im2; S22+=im2*im2; \
105  im1=ida1[y1m2][x1p]; im2=ida2[y2m2][x2p]; N++; S1+=im1; S2+=im2; S11+=im1*im1; S12+=im1*im2; S22+=im2*im2; \
106  im1=ida1[y1m2][x1p2]; im2=ida2[y2m2][x2p2]; N++; S1+=im1; S2+=im2; S11+=im1*im1; S12+=im1*im2; S22+=im2*im2; \
107  im1=ida1[y1m][x1m2]; im2=ida2[y2m][x2m2]; N++; S1+=im1; S2+=im2; S11+=im1*im1; S12+=im1*im2; S22+=im2*im2; \
108  im1=ida1[y1m][x1m]; im2=ida2[y2m][x2m]; N++; S1+=im1; S2+=im2; S11+=im1*im1; S12+=im1*im2; S22+=im2*im2; \
109  im1=ida1[y1m][x1]; im2=ida2[y2m][x2]; N++; S1+=im1; S2+=im2; S11+=im1*im1; S12+=im1*im2; S22+=im2*im2; \
110  im1=ida1[y1m][x1p]; im2=ida2[y2m][x2p]; N++; S1+=im1; S2+=im2; S11+=im1*im1; S12+=im1*im2; S22+=im2*im2; \
111  im1=ida1[y1m][x1p2]; im2=ida2[y2m][x2p2]; N++; S1+=im1; S2+=im2; S11+=im1*im1; S12+=im1*im2; S22+=im2*im2; \
112  im1=ida1[y1][x1m2]; im2=ida2[y2][x2m2]; N++; S1+=im1; S2+=im2; S11+=im1*im1; S12+=im1*im2; S22+=im2*im2; \
113  im1=ida1[y1][x1m]; im2=ida2[y2][x2m]; N++; S1+=im1; S2+=im2; S11+=im1*im1; S12+=im1*im2; S22+=im2*im2; \
114  im1=ida1[y1][x1]; im2=ida2[y2][x2]; N++; S1+=im1; S2+=im2; S11+=im1*im1; S12+=im1*im2; S22+=im2*im2; \
115  im1=ida1[y1][x1p]; im2=ida2[y2][x2p]; N++; S1+=im1; S2+=im2; S11+=im1*im1; S12+=im1*im2; S22+=im2*im2; \
116  im1=ida1[y1][x1p2]; im2=ida2[y2][x2p2]; N++; S1+=im1; S2+=im2; S11+=im1*im1; S12+=im1*im2; S22+=im2*im2; \
117  im1=ida1[y1p][x1m2]; im2=ida2[y2p][x2m2]; N++; S1+=im1; S2+=im2; S11+=im1*im1; S12+=im1*im2; S22+=im2*im2; \
118  im1=ida1[y1p][x1m]; im2=ida2[y2p][x2m]; N++; S1+=im1; S2+=im2; S11+=im1*im1; S12+=im1*im2; S22+=im2*im2; \
119  im1=ida1[y1p][x1]; im2=ida2[y2p][x2]; N++; S1+=im1; S2+=im2; S11+=im1*im1; S12+=im1*im2; S22+=im2*im2; \
120  im1=ida1[y1p][x1p]; im2=ida2[y2p][x2p]; N++; S1+=im1; S2+=im2; S11+=im1*im1; S12+=im1*im2; S22+=im2*im2; \
121  im1=ida1[y1p][x1p2]; im2=ida2[y2p][x2p2]; N++; S1+=im1; S2+=im2; S11+=im1*im1; S12+=im1*im2; S22+=im2*im2; \
122  im1=ida1[y1p2][x1m2]; im2=ida2[y2p2][x2m2]; N++; S1+=im1; S2+=im2; S11+=im1*im1; S12+=im1*im2; S22+=im2*im2; \
123  im1=ida1[y1p2][x1m]; im2=ida2[y2p2][x2m]; N++; S1+=im1; S2+=im2; S11+=im1*im1; S12+=im1*im2; S22+=im2*im2; \
124  im1=ida1[y1p2][x1]; im2=ida2[y2p2][x2]; N++; S1+=im1; S2+=im2; S11+=im1*im1; S12+=im1*im2; S22+=im2*im2; \
125  im1=ida1[y1p2][x1p]; im2=ida2[y2p2][x2p]; N++; S1+=im1; S2+=im2; S11+=im1*im1; S12+=im1*im2; S22+=im2*im2; \
126  im1=ida1[y1p2][x1p2]; im2=ida2[y2p2][x2p2]; N++; S1+=im1; S2+=im2; S11+=im1*im1; S12+=im1*im2; S22+=im2*im2; \
127  re =((double)N*(double)S12 - (double)S1*(double)S2)/ \
128  sqrt(((double)N*(double)S11 - (double)S1*(double)S1) * \
129  ((double)N*(double)S22 - (double)S2*(double)S2)); \
130  if (re!=re) re=0.0;
131 
132 #define hamming_distance(a, b, re) \
133  bool ab, bb; re=0; \
134  for(int i = 0 ; i < 8; i++ ) { \
135  ab = (a & ( 1 << i )) ? 1 : 0; \
136  bb = (b & ( 1 << i )) ? 1 : 0; \
137  if( (ab && !bb) || (!ab && bb) ) { \
138  re++; \
139  } \
140  }
141 
142  /** @class RegionMatcher
143  @ingroup g_corresp
144  @brief Basic functions for CornerMatcher.
145  class that holds sum of square diffeences (ssd),
146  normalized cross correlation (ncc) and KLT and parabola interpolation
147  for a convenient interface working with HomgPoint2D
148  doing boundary checks use CornerMatcher
149  @author Felix Woelk
150  */
151  class BIASMatcher2D_EXPORT RegionMatcher : public Debug
152  {
153  public:
154  RegionMatcher();
155 
156  ~RegionMatcher();
157 
158  /** Fast bilinear interpolation of a region in grey image
159  @author woelk 12/2007 (c) www.vision-n.de */
160  template<class StorageType>
161  static bool
162  BilinearRegion(const Image<StorageType>& im, const double &x,
163  const double &y, const unsigned hws, Matrix<double>& res);
164 
165  /** Fast bilinear interpolation of a region in interleaved color image
166  @author woelk 12/2007 (c) www.vision-n.de */
167  template<class StorageType>
168  static bool
169  BilinearRegionColor(const Image<StorageType>& im, const double &x,
170  const double &y, const unsigned hws, Matrix<double>& res);
171 
172  /** Fast bilinear interpolation of a region in interleaved color image
173  @author woelk 12/2007 (c) www.vision-n.de */
174  template<class StorageType>
175  static bool
176  BilinearRegionColor3(const Image<StorageType>& im, const double &x,
177  const double &y, const unsigned hws, Matrix<double>& res);
178 
179  /** Fast linear interpolation of a region in the image, untested
180  @author woelk 01/2008 (c) www.vision-n.de */
181  template<class StorageType>
182  static bool
183  LinearRegionX(const Image<StorageType>& im, const double &x,
184  const unsigned &y, const unsigned hws, Matrix<double>& res);
185 
186  /** fast sum of square differences (SSD) sim. measurment
187  between p1 and p2 without boundary checking
188  specialized version for unsigned char exist
189  for halfwindowsizes <= 127 if sizeof(unsigned char)=1 and
190  sizeof(unsigned long int)=4
191  @author woelk ? untested */
192  template<class StorageType>
193  void
194  SSD(const unsigned int p1[2], const unsigned int p2[2],
195  const StorageType **ida1, const StorageType **ida2,
196  const unsigned int halfwinsize, double& result) const;
197 
198  /** fast sum of absolut differences (SAD) sim. measurement
199  between p1 and p2 without boundary checking
200  specialized version for unsigned char exist only
201  for halfwindowsizes < 2050
202  @author woelk 02 2003 untested */
203  template<class StorageType>
204  void
205  SAD(const unsigned int p1[2], const unsigned int p2[2],
206  const StorageType **ida1, const StorageType **ida2,
207  const unsigned int halfwinsize, double& result) const;
208 
209  /** fast sum of absolut differences (SAD) sim. measurement
210  between p1 and p2 without boundary checking
211  Only for interleavd color images
212  @author woelk 02 2008 untested */
213  template<class StorageType>
214  void
215  SAD(const unsigned int p1[2], const unsigned int p2[2],
216  const StorageType **ida1, const StorageType **ida2,
217  const unsigned int halfwinsize, const unsigned channel_count,
218  double& result) const;
219 
220  /** 'normalized' version of SAD, returns 1.0 incase of no differences
221  0.0 in case of max difference
222  specialized version for unsigned char exist only
223  for halfwindowsizes < 2050
224  @author woelk 05 2003 untested */
225  template<class StorageType>
226  void
227  SADN(const unsigned int p1[2], const unsigned int p2[2],
228  const StorageType **ida1, const StorageType **ida2,
229  const unsigned int halfwinsize, double& result) const;
230 
231  /** fast ncc between p1 and p2 without boundary checking
232  Faster specialisation for unsigned char exists
233  because of integer arithmetic
234  Specialisation for unsigned char does work without overflow with
235  halfwindowsizes <= 127 if sizeof(unsigned char)=1 and
236  sizeof(unsigned long int)=4
237  @author woelk */
238  template<class StorageType>
239  void
240  NCC(const unsigned int p1[2], const unsigned int p2[2],
241  const StorageType **ida1, const StorageType **ida2,
242  const unsigned int halfwinsize, double& result) const;
243 
244  /** fast interpolated ncc
245  @author woelk */
246  template<class StorageType>
247  void
248  NCC(const double p1[2], const double p2[2], const StorageType **ida1,
249  const StorageType **ida2, const unsigned halfwinsize, double &result);
250 
251  /** fast ncc between p1 in ida1 and p2 in ida2
252  the correlation is calculated in a window of width=2*hww+1
253  and height = 2*hwh+1
254  only those pixels within the window are used, where the according roi
255  is not zero
256  @author woelk 09 2003 , tested */
257  template<class StorageType>
258  void
259  NCC(const unsigned p1[2], const unsigned p2[2], const StorageType **ida1,
260  const StorageType **roi1, const StorageType **ida2,
261  const StorageType **roi2, const unsigned hww, const unsigned hwh,
262  double& result) const;
263 
264  /** ColorNCC for interleaved HSV images
265  @author Daniel Grest, March 2003 */
266  template<class StorageType>
267  void
268  ColorNCC(const unsigned int p1[2], const unsigned int p2[2],
269  const StorageType **ida1, const StorageType **ida2,
270  const unsigned int halfwinsize, double &result) const;
271 
272  /** computes the Cross-CoVariance ,which is equal to NCC
273  but does no variance(contrast) normailzation
274  @author Daniel Grest, April 2003 */
275  template<class StorageType>
276  void
277  ColorCCV(const unsigned int p1[2], const unsigned int p2[2],
278  const StorageType **ida1, const StorageType **ida2,
279  const unsigned int halfwinsize, double &result) const;
280 
281  /** ColorNCC for interleaved hsL images
282  @author Daniel Grest, June 2003 */
283  template<class StorageType>
284  void
285  CNCC(const unsigned int p1[2], const unsigned int p2[2],
286  const StorageType **ida1, const StorageType **ida2,
287  const unsigned int halfwinsize, double &result) const;
288 
289  /** Searches in a window around p2 [p2[i]+-halfsearchwinsize] for a
290  point with correlation > mincorrelation.
291  If point with correlation > mincorrelation is found, the search is
292  aborted and resultpoint is set.
293  If no point with correlation > mincorrelation is found, a negative
294  value is returned. resultpoint is set to the point with biggest
295  correlation within search window, result is set tot this correlation.
296  Searching follows an outbound spiral around p2, either until
297  searchwinsize is completly searched or until a point with
298  correlation greater than mincorrelation is found
299  If whole search area should be used, use mincorrelation>=1.0
300  Searchwindow is always odd-sized (p2+/-halfsearchwinsize)
301  returns 0 on success, negative value if no point found with
302  correlation > mincorrelation
303  !!! warning, no boundary checking is performed,
304  i.e. p1[i]/p2[i] +-(halfnccwinsize+halfsearchwinsize) must lie
305  within StorageType-array
306  input: p1, p2, ida1, ida2, halfnccwinsize, halfsearchwinsize,
307  mincorrelation
308  output: resultpoint, result
309  */
310  template<class StorageType>
311  int
312  NCCSearchOdd(const unsigned int p1[2], const unsigned int p2[2],
313  const StorageType **ida1, const StorageType **ida2,
314  const unsigned int halfnccwinsize,
315  const unsigned int halfsearchwinsize, unsigned int resultpoint[2],
316  const double mincorrelation, double& result) const;
317 
318  /** As above but middle point of searchwindow is p1[0]+0.5, p1[1]+0.5
319  Searchwindow is always even-sized floor(p2[i]+0.5+-halfsearchwinsize)
320  */
321  template<class StorageType>
322  int
323  NCCSearchEven(const unsigned int p1[2], const unsigned int p2[2],
324  const StorageType **ida1, const StorageType **ida2,
325  const unsigned int halfnccwinsize,
326  const unsigned int halfsearchwinsize, unsigned int resultpoint[2],
327  const double mincorrelation, double& result) const;
328  /** fast hamming distance between p1 and p2 without boundary checking.
329  parameter censussize is number of channels in which census information is stored.
330  @author fkellner */
331  template<class StorageType>
332  void
333  HammingDistance(const unsigned int p1[2], const unsigned int p2[2],
334  const StorageType **ida1, const StorageType **ida2,
335  const unsigned int halfwinsize, const unsigned int censussize,
336  double& result) const;
337  template <class StorageType>
338  void
339  SADSamplingInsensitive(const unsigned int p1[2], const unsigned int p2[2],
340  const StorageType **ida1, const StorageType **ida2,
341  const unsigned int halfwinsize, double& result) const;
342 
343  // /** find correspondence to p1 from ida1 in ida2
344  // at aproximate starting position p2 using KLT algorithm:
345 
346  // Solve G * d = e with
347 
348  // G = | sum(gx*gx) sum(gx*gy) | , and
349  // | sum(gx*gy) sum(gy*gy) |
350 
351  // e = | sum(gx*(I2-I1)) sum(gy*(I2-I1)) |^T
352 
353  // sums are over the window, gx is x-gradient and gy the y-gradient in
354  // first image, I1 and I2 are the grey values in the first
355  // resp. the second image
356  // iterate numiterations times
357  // !!! no initial boundary checks is performed in this algorithm !!!
358  // returns 0 on success
359  // returns -1 if point slids out of [0, with][0, height] modulo halfwinsize
360  // output: result - image coo of p1 in im2
361  // error - norm(e)
362  // input : p1 - location of feature point in image1
363  // p2 - aproximate location of same feature point in new image 2
364  // could be set to p1 featurepoint is not moving a lot
365  // ida1 - pointer to the image data (row major) of first image
366  // ida2 - pointer to the image data (row major) of new image
367  // gradx - pointer to the x-gradient data of image1
368  // grady - pointer to the y-gradient data of image1
369  // halfwinsize - the window to use
370  // numiterations - number of iterations to do
371  // width - width of StorageType Arrays
372  // height - height of StorageType Array ida1[height-1][with-1]
373  // is the maximal allowed acces to the arrays
374  // @author fw */
375  // template <class StorageType>
376  // int KLT(KLT_TYPE p1[2], KLT_TYPE p2[2], StorageType **ida1,
377  // StorageType **ida2, StorageType **gradx, StorageType **grady,
378  // unsigned int width, unsigned int height, unsigned int halfwinsize,
379  // unsigned int numiterations, KLT_TYPE result[2], KLT_TYPE& error);
380 
381  // /** same as KLT above but do it until either maxiterations is reached
382  // or the error is dropping below maxerror */
383  // template <class StorageType>
384  // int KLT(KLT_TYPE p1[2], KLT_TYPE p2[2], StorageType **ida1,
385  // StorageType **ida2, StorageType **gradx, StorageType **grady,
386  // unsigned int width, unsigned int height, unsigned int halfwinsize,
387  // unsigned int maxiterations, KLT_TYPE maxerror, KLT_TYPE result[2],
388  // KLT_TYPE& error);
389 
390  // /** as above but uses the averaged gradient of the 2 images */
391  // template <class StorageType>
392  // int KLT2(KLT_TYPE p1[2], KLT_TYPE p2[2], StorageType **ida1,
393  // StorageType **ida2, StorageType **gradx1, StorageType **grady1,
394  // StorageType **gradx2, StorageType **grady2,
395  // unsigned int width, unsigned int height, unsigned int halfwinsize,
396  // unsigned int maxiterations, KLT_TYPE maxerror, KLT_TYPE result[2],
397  // KLT_TYPE& error);
398 
399  // /** same as above but also returns the number of iterations used */
400  // template <class StorageType>
401  // int KLT2(KLT_TYPE p1[2], KLT_TYPE p2[2], StorageType **ida1,
402  // StorageType **ida2, StorageType **gradx1, StorageType **grady1,
403  // StorageType **gradx2, StorageType **grady2,
404  // unsigned int width, unsigned int height, unsigned int halfwinsize,
405  // unsigned int maxiterations, KLT_TYPE maxerror, KLT_TYPE result[2],
406  // KLT_TYPE& error, unsigned &iter);
407 
408  // /** as above but uses a 5x5 window
409  // @author woelk 08 2003 */
410  // template <class StorageType>
411  // int KLT5x5(KLT_TYPE p1[2], KLT_TYPE p2[2], StorageType **ida1,
412  // StorageType **ida2, StorageType **gradx1, StorageType **grady1,
413  // StorageType **gradx2, StorageType **grady2,
414  // unsigned width, unsigned height, unsigned maxiterations,
415  // KLT_TYPE maxerror, KLT_TYPE result[2], KLT_TYPE& error);
416 
417  // /** searches only perpendicular to the dominating edge, i.e.
418  // parallel to the dominating gradient
419  // the search direction is discretized in 4 different directions
420  // the discretization-error of direction is max 22 deg which
421  // accounts to 7% of length error
422  // decides which KLT_? to use
423  // @author woelk 08 2003 */
424 
425  // template <class StorageType>
426  // int KLTEdge(KLT_TYPE p1[2], KLT_TYPE p2[2], StorageType **ida1,
427  // StorageType **ida2, StorageType **gradx, StorageType **grady,
428  // StorageType **gradNWSE, StorageType **gradNESW,
429  // unsigned width, unsigned height, unsigned maxiterations,
430  // KLT_TYPE maxerror, KLT_TYPE result[2], KLT_TYPE& error);
431 
432  // template <class StorageType>
433  // int KLTEdge_NS(KLT_TYPE p1[2], KLT_TYPE p2[2], StorageType **ida1,
434  // StorageType **ida2, StorageType **grad,
435  // unsigned width, unsigned height, unsigned maxiterations,
436  // KLT_TYPE maxerror, KLT_TYPE result[2], KLT_TYPE& error);
437 
438  // template <class StorageType>
439  // int KLTEdge_WE(KLT_TYPE p1[2], KLT_TYPE p2[2], StorageType **ida1,
440  // StorageType **ida2, StorageType **grad,
441  // unsigned width, unsigned height, unsigned maxiterations,
442  // KLT_TYPE maxerror, KLT_TYPE result[2], KLT_TYPE& error);
443 
444  // template <class StorageType>
445  // int KLTEdge_NESW(KLT_TYPE p1[2], KLT_TYPE p2[2], StorageType **ida1,
446  // StorageType **ida2, StorageType **grad,
447  // unsigned width, unsigned height, unsigned maxiterations,
448  // KLT_TYPE maxerror, KLT_TYPE result[2], KLT_TYPE& error);
449 
450  // template <class StorageType>
451  // int KLTEdge_NWSE(KLT_TYPE p1[2], KLT_TYPE p2[2], StorageType **ida1,
452  // StorageType **ida2, StorageType **grad,
453  // unsigned width, unsigned height, unsigned maxiterations,
454  // KLT_TYPE maxerror, KLT_TYPE result[2], KLT_TYPE& error);
455 
456 
457  /** searches along line from start2 to end2 within +/- epsilon pixels
458  away from the line in ida2 for a match to p1
459  using ncc with halfwinsize
460  boundary checking is not done, i.e. start2 and end2
461  must lie more than halfwindowsize+epsilon away from the
462  image border
463  p1 must lie more than halfwinsize away from image border
464  only for those points in im2 a correlation is calculated for which
465  the gradient is > grad1[p1[1]][p1[0]]*gradientscale
466  input: ida1 - pointer to the image data array of im1, row mayor
467 
468 
469  output: result - image coo of point with greatest corelation
470  correlation - the correlation betwen p1 and result
471  returns: 0 on success
472  -1 no matching point found
473  <-1 error */
474  template<class StorageType>
475  int
476  LineMatcherNCC(const unsigned int p1[2], const unsigned int start2[2],
477  const unsigned int end2[2], const StorageType **ida1,
478  const StorageType **ida2, const StorageType **grad1,
479  const StorageType **grad2, const unsigned int halfwinsize,
480  const unsigned int epsilon, const double gradientscale,
481  unsigned int result[2], double& correlation) const;
482 
483  /** as above but uses SAD */
484  template<class StorageType>
485  int
486  LineMatcherSAD(const unsigned int p1[2], const unsigned int start2[2],
487  const unsigned int end2[2], const StorageType **ida1,
488  const StorageType **ida2, const StorageType **grad1,
489  const StorageType **grad2, const unsigned int halfwinsize,
490  const unsigned int epsilon, const double gradientscale,
491  unsigned int result[2], double& sad) const;
492 
493  /** as above but uses SSD */
494  template<class StorageType>
495  int
496  LineMatcherSSD(const unsigned int p1[2], const unsigned int start2[2],
497  const unsigned int end2[2], const StorageType **ida1,
498  const StorageType **ida2, const StorageType **grad1,
499  const StorageType **grad2, const unsigned int halfwinsize,
500  const unsigned int epsilon, const double gradientscale,
501  unsigned int result[2], double& ssd) const;
502 
503  /** function to achieve sub-pixel accuracy if p1 and p2
504  are corresponding image points from ida1 resp ida2
505  exact parabola approximation with 3 values is used
506  no boundary checks, p1 and p2 must be further than
507  halfwinsize+1 from array border
508  correlation between p1 and p2 must be a local discrete maximum
509  i.e. NCC(p1, {p2[0]+/-1, p2[1]+/-1})<NCC(p1, p2)
510  @author woelk 03 2003 */
511  template<class StorageType>
512  int
513  ParabolaNCC(const unsigned int p1[2], const unsigned int p2[2],
514  const StorageType **ida1, const StorageType **ida2,
515  const unsigned int halfwinsize, KLT_TYPE result[2]) const;
516 
517  /** function to achieve sub-pixel accuracy if p1 and p2
518  are corresponding image points from ida1 resp ida2
519  over-determined parabola approximation with 5 values is used
520  no boundary checks, p1 and p2 must be further than
521  halfwinsize+2 from array border
522  correlation between p1 and p2 must be a local discrete maximum
523  i.e. NCC(p1, {p2[0]+/-2, p2[1]+/-2})<NCC(p1, p2)
524  @author woelk 04 2003 */
525  template<class StorageType>
526  int
527  ParabolaNCC5(const unsigned int p1[2], const unsigned int p2[2],
528  const StorageType **ida1, const StorageType **ida2,
529  const unsigned int halfwinsize, KLT_TYPE result[2]) const;
530 
531  /*
532  template <class StorageType>
533  int GetTrackableRegions(StorageType **gx, StorageType **gy,
534  unsigned int halfwinsize, StorageType **tmap,
535  StorageType &min, StorageType&max);
536  */
537 
538  /**
539  @author frick
540  @brief calculates absolute difference between pixel (x1,y1) in image img1 and (x2,y2) in img2
541  */
542  template<class StorageType>
543  static float
544  AD(const unsigned int x1, const unsigned int y1,
545  BIAS::Image<StorageType>& img1, const unsigned int x2,
546  const unsigned int y2, BIAS::Image<StorageType>& img2);
547 
548  /**
549  @author fkellner
550  @brief calculates sampling insensitive dissimilarity measure (birchfield, tomasi)
551  */
552  template<class StorageType>
553  void
554  SamplingInsensitiveDissimilarity(const unsigned int p1[2], const unsigned int p2[2],
555  const StorageType **ida1, const StorageType **ida2, double &result) const;
556 
557  protected:
558  template<class StorageType>
559  KLT_TYPE
560  _Bilinear(const StorageType **ida, const KLT_TYPE x, const KLT_TYPE y) const;
561 
562  KLT_TYPE
563  _Bilinear(const unsigned char **ida, const KLT_TYPE x, const KLT_TYPE y) const;
564 
565  /** bilinear interpolation, also moved gradient image data from [0:255]
566  to [-127, 128] before interpolating */
567  // template <class StorageType>
568  // KLT_TYPE _BilinearGrad(StorageType **ida, KLT_TYPE x, KLT_TYPE y);
569 
570  template<class StorageType>
571  void
572  _BilinearRegion(const StorageType **ida, const KLT_TYPE x,
573  const KLT_TYPE y, const unsigned hws, Matrix<KLT_TYPE>& res) const;
574 
575  void
576  _BilinearRegion(const unsigned char **ida, const KLT_TYPE x,
577  const KLT_TYPE y, const unsigned hws, Matrix<KLT_TYPE>& res) const;
578 
579  /** bilinear interpolation, also moved gradient image data from [0:255]
580  to [-127, 128] before interpolating */
581  // template <class StorageType>
582  // void _BilinearRegionGrad(StorageType **ida, KLT_TYPE x, KLT_TYPE y,
583  // unsigned hws, Matrix<KLT_TYPE>& res);
584 
585  inline
586  KLT_TYPE
587  _ParabolaApprox(const KLT_TYPE x1, const KLT_TYPE x2, const KLT_TYPE x3) const
588  {
589  KLT_TYPE a, p;
590  a = (x1 - 2.0 * x2 + x3);
591  p = (x1 - x3);
592  return (a != 0.0) ? (KLT_TYPE) p / (2.0 * a) : (KLT_TYPE) 0.0;
593  }
594 
595  inline
596  KLT_TYPE
597  _ParabolaApprox(const KLT_TYPE x1, const KLT_TYPE x2, const KLT_TYPE x3,
598  const KLT_TYPE x4, const KLT_TYPE x5) const
599  {
600  Matrix<double> A(5, 3, "4 -2 1 1 -1 1 0 0 1 1 1 1 4 2 1");
601  Vector<double> b(5);
602  b[0] = x1;
603  b[1] = x2;
604  b[2] = x3;
605  b[3] = x4;
606  b[4] = x5;
607  Vector<double> res = Lapack_LLS_QR_linear_solve(A, b);
608  return (res[0] != 0.0) ? (KLT_TYPE) res[1] / (2.0 * res[0])
609  : (KLT_TYPE) 0.0;
610  }
611 
612  // internal vars for KLT
613  Matrix2x2<KLT_TYPE> _G, _Ginv;
614  Vector2<KLT_TYPE> _e, _d, _disp;
615  Matrix<KLT_TYPE> _gx, _gy, _bl1, _bl2, _gx1, _gy1, _gx2, _gy2, _gsx, _gsy;
616  Matrix<KLT_TYPE> _bl1_5, _bl2_5, _gx1_5, _gy1_5, _gx2_5, _gy2_5, _gsx_5,
617  _gsy_5;
618  unsigned _klthws, _winsize;
619  // internals vars for interpolated NCC
620  Matrix<KLT_TYPE> _iim1, _iim2;
621  unsigned _ncchws, _nccwinsize;
622 
623  };
624 } // namespace
625 #endif // __RegionMatcher_hh__
Lapack_LLS_QR_linear_solve
BIAS::Vector< double > Lapack_LLS_QR_linear_solve(const BIAS::Matrix< double > &A, const BIAS::Vector< double > &b, int &res)
linear least squares solves |Ax-b|=min res=0 on success, anything else rudimentary tested (woelk) ...
Definition: Lapack.cpp:311
BIAS::RegionMatcher::_Ginv
Matrix2x2< KLT_TYPE > _Ginv
Definition: RegionMatcher.hh:613
BIAS::RegionMatcher::_gy2
Matrix< KLT_TYPE > _gy2
Definition: RegionMatcher.hh:615
BIAS::Image
class BIASImageBase_EXPORT Image
Definition: ImageBase.hh:91
BIAS::RegionMatcher::_gy2_5
Matrix< KLT_TYPE > _gy2_5
Definition: RegionMatcher.hh:616
BIAS::Image
The image template class for specific storage types.
Definition: Image.hh:78
BIAS::RegionMatcher
Basic functions for CornerMatcher.
Definition: RegionMatcher.hh:151
BIAS::RegionMatcher::_iim2
Matrix< KLT_TYPE > _iim2
Definition: RegionMatcher.hh:620
BIAS::RegionMatcher::_e
Vector2< KLT_TYPE > _e
Definition: RegionMatcher.hh:614
BIAS::RegionMatcher::_ParabolaApprox
KLT_TYPE _ParabolaApprox(const KLT_TYPE x1, const KLT_TYPE x2, const KLT_TYPE x3) const
bilinear interpolation, also moved gradient image data from [0:255] to [-127, 128] before interpolati...
Definition: RegionMatcher.hh:587
BIAS::RegionMatcher::_ParabolaApprox
KLT_TYPE _ParabolaApprox(const KLT_TYPE x1, const KLT_TYPE x2, const KLT_TYPE x3, const KLT_TYPE x4, const KLT_TYPE x5) const
Definition: RegionMatcher.hh:597
Generated on Tue Feb 6 2018 02:01:40 for Basic Image AlgorithmS Library by   doxygen 1.8.5