Basic Image AlgorithmS Library  2.8.0
 All Classes Namespaces Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
SimilarityMeasure.cpp
1 /* This file is part of the BIAS library (Basic ImageAlgorithmS).
2 
3  Copyright (C) 2003-2009 (see file CONTACTS 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 #include "RegionMatcher.hh"
22 
23 #include <float.h>
24 
25 namespace BIAS {
26 using namespace std;
27 
28 template<class StorageType>
29 void RegionMatcher::SAD(const unsigned int p1[2], const unsigned int p2[2],
30  const StorageType **ida1, const StorageType **ida2,
31  const unsigned int halfwinsize, double& result) const {
32  int nhalfwinsize = -int(halfwinsize);
33  result = 0.0;
34  register StorageType st1, st2;
35  for (register int c = nhalfwinsize; c <= (int) halfwinsize; c++) {
36  for (register int r = nhalfwinsize; r <= (int) halfwinsize; r++) {
37  st1 = ida1[p1[1] + r][p1[0] + c];
38  st2 = ida2[p2[1] + r][p2[0] + c];
39  result += (double) ((st1 >= st2) ? st1 - st2 : st2 - st1);
40  }
41  }
42 }
43 
44 template<> BIASMatcher2D_EXPORT
45 void RegionMatcher::SAD<unsigned char>(const unsigned int p1[2],
46  const unsigned int p2[2], const unsigned char **ida1,
47  const unsigned char **ida2, const unsigned int halfwinsize,
48  double& result) const {
49  int nhalfwinsize = -int(halfwinsize);
50  register long unsigned int res = 0;
51  register unsigned char c1, c2;
52  for (register int c = nhalfwinsize; c <= (int) halfwinsize; c++) {
53  for (register int r = nhalfwinsize; r <= (int) halfwinsize; r++) {
54  c1 = ida1[p1[1] + r][p1[0] + c];
55  c2 = ida2[p2[1] + r][p2[0] + c];
56  res += (c1 >= c2) ? c1 - c2 : c2 - c1;
57  }
58  }
59  result = (double) res;
60 }
61 
62 template<class StorageType>
63 void RegionMatcher::SAD(const unsigned int p1[2], const unsigned int p2[2],
64  const StorageType **ida1, const StorageType **ida2,
65  const unsigned halfwinsize, const unsigned channel_count,
66  double& result) const {
67  int nhalfwinsize = -int(halfwinsize);
68  result = 0.0;
69  register StorageType st1, st2;
70  register unsigned channel;
71  for (register int c = nhalfwinsize; c <= (int) halfwinsize; c++) {
72  for (register int r = nhalfwinsize; r <= (int) halfwinsize; r++) {
73  for (channel = 0; channel < channel_count; channel++) {
74  st1 = ida1[p1[1] + r][(p1[0] + c) * channel_count + channel];
75  st2 = ida2[p2[1] + r][(p2[0] + c) * channel_count + channel];
76  result += (double) ((st1 >= st2) ? st1 - st2 : st2 - st1);
77  }
78  }
79  }
80 }
81 
82 template<> BIASMatcher2D_EXPORT
83 void RegionMatcher::SAD<unsigned char>(const unsigned int p1[2],
84  const unsigned int p2[2], const unsigned char **ida1,
85  const unsigned char **ida2, const unsigned halfwinsize,
86  const unsigned channel_count, double& result) const {
87  int nhalfwinsize = -int(halfwinsize);
88  result = 0.0;
89  register unsigned char st1, st2;
90  register unsigned channel;
91  for (register int c = nhalfwinsize; c <= (int) halfwinsize; c++) {
92  for (register int r = nhalfwinsize; r <= (int) halfwinsize; r++) {
93  for (channel = 0; channel < channel_count; channel++) {
94  st1 = ida1[p1[1] + r][(p1[0] + c) * channel_count + channel];
95  st2 = ida2[p2[1] + r][(p2[0] + c) * channel_count + channel];
96  result += (double) ((st1 >= st2) ? st1 - st2 : st2 - st1);
97  }
98  }
99  }
100 }
101 
102 template<class StorageType>
103 void RegionMatcher::SADN(const unsigned int p1[2], const unsigned int p2[2],
104  const StorageType **ida1, const StorageType **ida2,
105  const unsigned int halfwinsize, double& result) const {
106  int nhalfwinsize = -int(halfwinsize);
107  result = 0.0;
108  register StorageType st1, st2;
109  for (register int c = nhalfwinsize; c <= (int) halfwinsize; c++) {
110  for (register int r = nhalfwinsize; r <= (int) halfwinsize; r++) {
111  st1 = ida1[p1[1] + r][p1[0] + c];
112  st2 = ida2[p2[1] + r][p2[0] + c];
113  result += (st1 >= st2) ? st1 - st2 : st2 - st1;
114  }
115  }
116  // STORAGETYPE_MAX should be the biggest poistive number that can be
117  // expressed with this storage type
118 #define STORAGETYPE_MAX 256.0
119 
120  unsigned tmp = (halfwinsize << 1) + 1;
121  result /= double(tmp * tmp * STORAGETYPE_MAX);
122  result = 1.0 - result;
123 }
124 
125 template<> BIASMatcher2D_EXPORT
126 void RegionMatcher::SADN<unsigned char>(const unsigned int p1[2],
127  const unsigned int p2[2], const unsigned char **ida1,
128  const unsigned char **ida2, const unsigned int halfwinsize,
129  double& result) const {
130  int nhalfwinsize = -int(halfwinsize);
131  register long unsigned int res = 0;
132  register unsigned char c1, c2;
133  for (register int c = nhalfwinsize; c <= (int) halfwinsize; c++) {
134  for (register int r = nhalfwinsize; r <= (int) halfwinsize; r++) {
135  c1 = ida1[p1[1] + r][p1[0] + c];
136  c2 = ida2[p2[1] + r][p2[0] + c];
137  res += (c1 >= c2) ? c1 - c2 : c2 - c1;
138  }
139  }
140  // result = 1.0 - res /(UCHAR_MAX*windowarea) = 1 - res / tmp
141  // tmp = sqrt(UCHAR_MAX)*(2*halfwinsize+1) ca. 16 * (2 * halfwinsize + 1)
142  // = 32 * halfwinsize + 16
143  unsigned tmp = (halfwinsize << 5) + 16;
144  result = 1.0 - (double) res / double(tmp * tmp);
145 }
146 
147 template<class StorageType>
148 void RegionMatcher::SSD(const unsigned int p1[2], const unsigned int p2[2],
149  const StorageType **ida1, const StorageType **ida2,
150  const unsigned int halfwinsize, double& result) const {
151  int nhalfwinsize = -int(halfwinsize);
152  result = 0.0;
153  double sd;
154  for (register int c = nhalfwinsize; c <= (int) halfwinsize; c++) {
155  for (register int r = nhalfwinsize; r <= (int) halfwinsize; r++) {
156  //DONT USE POW(.,2.0)!!!!
157  // result += pow(
158  // (double) (ida1[p1[1] + r][p1[0] + c]) -(double) (ida2[p2[1]
159  // + r][p2[0] + c]), 2.0);
160  sd = double(ida1[p1[1] + r][p1[0] + c])
161  - double(ida2[p2[1] + r][p2[0] + c]);
162  sd *= sd;
163  result += sd;
164  }
165  }
166 }
167 
168 template<> BIASMatcher2D_EXPORT
169 void RegionMatcher::SSD<unsigned char>(const unsigned int p1[2],
170  const unsigned int p2[2], const unsigned char **ida1,
171  const unsigned char **ida2, const unsigned int halfwinsize,
172  double& result) const {
173  int nhalfwinsize = -int(halfwinsize);
174  register long unsigned int res = 0;
175  register short diff;
176  for (register int c = nhalfwinsize; c <= (int) halfwinsize; c++) {
177  for (register int r = nhalfwinsize; r <= (int) halfwinsize; r++) {
178  diff = (short) (ida1[p1[1] + r][p1[0] + c]) -(short) (ida2[p2[1]
179  + r][p2[0] + c]);
180  res += (unsigned int) (diff * diff);
181  }
182  }
183  result = (double) res;
184 }
185 
186 template<class StorageType>
187 void RegionMatcher::NCC(const unsigned int p1[2], const unsigned int p2[2],
188  const StorageType **ida1, const StorageType **ida2,
189  const unsigned int halfwinsize, double& result) const {
190  int nhalfwinsize = -int(halfwinsize);
191  double im1 = 0, im2 = 0, S1 = 0, S2 = 0;
192  double S11 = 0, S12 = 0, S22 = 0, N = 0;
193  BIASASSERT(ida1!=NULL);
194  BIASASSERT(ida2!=NULL);
195  for (register int r = nhalfwinsize; r <= (int) halfwinsize; r++)
196  for (register int c = nhalfwinsize; c <= (int) halfwinsize; c++) {
197  im1 = (double) (ida1[p1[1] + r][p1[0] + c]);
198  im2 = (double) (ida2[p2[1] + r][p2[0] + c]);
199  N++;
200  S1 += im1;
201  S2 += im2;
202  S11 += im1 * im1;
203  S12 += im1 * im2;
204  S22 += im2 * im2;
205  }
206  result = (N * S12 - S1 * S2) / sqrt((N * S11 - S1 * S1) * (N * S22 - S2
207  * S2));
208  if (result != result)
209  result = 0.0;
210 }
211 
212 // faster specialisation because of integer arithmetic
213 // does work without overflow with halfwindowsizes <= 127
214 // if sizeof(unsigned char)=1 and sizeof(unsigned long int)=4
215 template<> BIASMatcher2D_EXPORT void RegionMatcher::NCC<unsigned char>(
216  const unsigned int p1[2], const unsigned int p2[2],
217  const unsigned char **ida1, const unsigned char **ida2,
218  const unsigned int halfwinsize, double& result) const {
219  int nhalfwinsize = -int(halfwinsize);
220  register int im1 = 0, im2 = 0;
221  register unsigned long int S11 = 0, S12 = 0, S22 = 0, N = 0, S1 = 0, S2 = 0;
222  for (register int r = nhalfwinsize; r <= (int) halfwinsize; r++)
223  for (register int c = nhalfwinsize; c <= (int) halfwinsize; c++) {
224  im1 = (int) (ida1[p1[1] + r][p1[0] + c]);
225  im2 = (int) (ida2[p2[1] + r][p2[0] + c]);
226  // if (im1!=im2)
227  // cerr << "im1 ("<<p1[0]+c<<", "<<p1[1]+r<<") : "<<im1
228  // << "\tim2 ("<<p2[0]+c<<", "<<p2[1]+r<<") : "<<im2<<endl;
229  N += 1;
230  S1 += im1;
231  S2 += im2;
232  S11 += im1 * im1;
233  S12 += im1 * im2;
234  S22 += im2 * im2;
235  }
236  // cerr << " N " << N << " S1 " << S1 << " S11 " << S11 << " S2 " << S2
237  // << " S22 " << S22 << " S12 " << S12 << endl;
238  result = ((double) N * (double) S12 - (double) S1 * (double) S2) / sqrt(
239  ((double) N * (double) S11 - (double) S1 * (double) S1)
240  * ((double) N * (double) S22 - (double) S2 * (double) S2));
241  if (result != result)
242  result = 0.0;
243 
244  // cout<<result<<" at ("<<p1[0]<<", "<<p1[1]<<") & ("
245  // <<p2[0]<<", "<<p2[1]<<")"<<endl;
246 }
247 
248 template<class StorageType>
249 void RegionMatcher::NCC(const double p1[2], const double p2[2],
250  const StorageType **ida1, const StorageType **ida2, const unsigned hws,
251  double& result) {
252  double im1 = 0, im2 = 0, S1 = 0, S2 = 0;
253  double S11 = 0, S12 = 0, S22 = 0, N = 0;
254  if (_ncchws != hws) {
255  _ncchws = hws;
256  _nccwinsize = 1 + (hws << 1);
257  _iim1.newsize(_nccwinsize, _nccwinsize);
258  _iim2.newsize(_nccwinsize, _nccwinsize);
259  }
260  _BilinearRegion(ida1, p1[0], p1[1], _ncchws, _iim1);
261  _BilinearRegion(ida2, p2[0], p2[1], _ncchws, _iim2);
262  for (register unsigned r = 0; r <= hws; r++)
263  for (register unsigned c = 0; c <= hws; c++) {
264  im1 = (double) _iim1[r][c];
265  im2 = (double) _iim2[r][c];
266  N++;
267  S1 += im1;
268  S2 += im2;
269  S11 += im1 * im1;
270  S12 += im1 * im2;
271  S22 += im2 * im2;
272  }
273  result = (N * S12 - S1 * S2) / sqrt((N * S11 - S1 * S1) * (N * S22 - S2
274  * S2));
275  if (result != result)
276  result = 0.0;
277 }
278 
279 template<class StorageType>
280 void RegionMatcher::NCC(const unsigned p1[2], const unsigned p2[2],
281  const StorageType **ida1, const StorageType **roi1,
282  const StorageType **ida2, const StorageType **roi2, const unsigned hww,
283  const unsigned hwh, double& result) const {
284  int nhww = -int(hww);
285  int nhwh = -int(hwh);
286  double im1 = 0, im2 = 0, S1 = 0, S2 = 0;
287  double S11 = 0, S12 = 0, S22 = 0, N = 0;
288  for (register int r = nhwh; r < (int) hwh; r++) {
289  for (register int c = nhww; c < (int) hww; c++) {
290  if (roi1[p1[1] + r][p1[0] + c] != 0 && roi2[p2[1] + r][p2[0] + c]
291  != 0) {
292  im1 = (double) (ida1[p1[1] + r][p1[0] + c]);
293  im2 = (double) (ida2[p2[1] + r][p2[0] + c]);
294  N++;
295  S1 += im1;
296  S2 += im2;
297  S11 += im1 * im1;
298  S12 += im1 * im2;
299  S22 += im2 * im2;
300  }
301  }
302  }
303  result = (N * S12 - S1 * S2) / sqrt((N * S11 - S1 * S1) * (N * S22 - S2
304  * S2));
305  if (result != result)
306  result = 0.0;
307 }
308 
309 template<class StorageType>
310 void RegionMatcher::ColorCCV(const unsigned int p1[2],
311  const unsigned int p2[2], const StorageType **ida1,
312  const StorageType **ida2, const unsigned int halfwinsize,
313  double &result) const {
314  double Sum1 = 0, Sh1 = 0, Ss1 = 0;
315  double Sum2 = 0, Sh2 = 0, Ss2 = 0;
316  double S12 = 0, N = 0;
317  StorageType *H1, *S1, *H2, *S2;
318  double I1, I2;
319  double angle, h1, s1, h2, s2;
320  for (register int r = -halfwinsize; r <= (int) halfwinsize; r++)
321  for (register int c = -halfwinsize; c <= (int) halfwinsize; c++) {
322  H1 = (ida1[(p1[1] + r)] + ((p1[0] + c) * 3));
323  S1 = H1 + 1;
324  I1 = *(S1 + 1);
325  H2 = (ida2[(p2[1] + r)] + ((p2[0] + c) * 3));
326  S2 = H2 + 1;
327  I2 = *(S2 + 1);
328  BIASDOUT(D_REGION_MATCHER_COLORCCV,
329  "HSI1:"<<(int)(*H1)<<" "<<(int)(*S1)<<" "<<(int)I1<<
330  " HSI2:"<<(int)(*H2)<<" "<<(int)(*S2)<<" "<<(int)I2);
331  // convert H and S to euclidean coords, H=0 is red
332  angle = M_PI * (float) (*H1) / 127.5;
333  h1 = cos(angle) * (*S1);
334  s1 = sin(angle) * (*S1);
335  angle = M_PI * (float) (*H2) / 127.5;
336  h2 = cos(angle) * (*S2);
337  s2 = sin(angle) * (*S2);
338  BIASDOUT(D_REGION_MATCHER_COLORCCV,
339  "[h1,s1]:["<<(int)h1<<","<<(int)s1<<
340  "] [h2,s2]:["<<(int)h2<<","<<(int)s2);
341  angle = h1 * h2 + s1 * s2;
342  if (angle < 0)
343  angle = 0;
344  S12 += angle + I1 * I2;
345  Sh1 += h1;
346  Ss1 += s1;
347  Sum1 += I1;
348  Sh2 += h2;
349  Ss2 += s2;
350  Sum2 += I2;
351  N++;
352  }
353  double S1S2 = Sh1 * Sh2 + Ss1 * Ss2;
354  if (S1S2 < 0)
355  S1S2 = 0;
356  S1S2 += Sum1 * Sum2;
357  BIASDOUT(D_REGION_MATCHER_COLORCCV,
358  "S12:"<<S12<<" Sh1,Ss1:"<<Sh1<<","<<Ss1<<" Sum1:"<<Sum1<<
359  " Sh2,Ss2:"<<Sh2<<","<<Ss2<<" Sum2:"<<Sum2);
360  result = (S12 - S1S2 / N);
361 }
362 
363 template<class StorageType>
364 void RegionMatcher::ColorNCC(const unsigned int p1[2],
365  const unsigned int p2[2], const StorageType **ida1,
366  const StorageType **ida2, const unsigned int halfwinsize,
367  double &result) const {
368  int nhalfwinsize = -int(halfwinsize);
369 #define CNCC_EPSILON 1E-08
370  double Sum1 = 0, Sh1 = 0, Ss1 = 0;
371  double Sum2 = 0, Sh2 = 0, Ss2 = 0;
372  double S11 = 0, S12 = 0, S22 = 0, N = 0;
373  const StorageType *H1, *S1, *H2, *S2;
374  double I1, I2;
375  double angle, h1, s1, h2, s2;
376  for (register int r = nhalfwinsize; r <= (int) halfwinsize; r++)
377  for (register int c = nhalfwinsize; c <= (int) halfwinsize; c++) {
378  H1 = (ida1[(p1[1] + r)] + ((p1[0] + c) * 3));
379  S1 = H1 + 1;
380  I1 = *(S1 + 1);
381  H2 = (ida2[(p2[1] + r)] + ((p2[0] + c) * 3));
382  S2 = H2 + 1;
383  I2 = *(S2 + 1);
384  BIASDOUT(D_REGION_MATCHER_COLORNCC,
385  "HSI1:"<<(int)(*H1)<<" "<<(int)(*S1)<<" "<<(int)I1<<
386  " HSI2:"<<(int)(*H2)<<" "<<(int)(*S2)<<" "<<(int)I2);
387  // convert H and S to euclidean coords, H=0 is red
388  angle = M_PI * (float) (*H1) / 127.5;
389  h1 = cos(angle) * (float) (*S1);
390  s1 = sin(angle) * (float) (*S1);
391  angle = M_PI * (float) (*H2) / 127.5;
392  h2 = cos(angle) * (float) (*S2);
393  s2 = sin(angle) * (float) (*S2);
394  BIASDOUT(D_REGION_MATCHER_COLORNCC,
395  "[h1,s1]:["<<(int)h1<<","<<(int)s1<<
396  "] [h2,s2]:["<<(int)h2<<","<<(int)s2);
397  angle = h1 * h2 + s1 * s2; // scalar product
398  // angle = fabs(double(*H1)-double(*H2));
399  // if (angle>180)
400  // angle=255-angle;
401  // angle = cos(angle/255.0) * double(*S1)*double(*S2); // scalar product
402  if (angle < 0)
403  angle = 0;
404  S12 += angle + I1 * I2;
405  Sh1 += h1;
406  Ss1 += s1;
407  Sum1 += I1;
408  Sh2 += h2;
409  Ss2 += s2;
410  Sum2 += I2;
411  S11 += h1 * h1 + s1 * s1 + I1 * I1;
412  S22 += h2 * h2 + s2 * s2 + I2 * I2;
413  N++;
414  }
415 
416  double S1S2 = Sh1 * Sh2 + Ss1 * Ss2;
417  if (S1S2 < 0)
418  S1S2 = 0;
419  S1S2 = Sum1 * Sum2;
420  double S1S1 = Sum1 * Sum1;
421  double S2S2 = Sum2 * Sum2;
422  BIASDOUT(D_REGION_MATCHER_COLORNCC,
423  "N*S12:"<<N*S12<<" S1S2:"<<S1S2<<endl<<
424  " Sh1,Ss1:"<<Sh1<<","<<Ss1<<" Sum1:"<<Sum1<<
425  " Sh2,Ss2:"<<Sh2<<","<<Ss2<<" Sum2:"<<Sum2);
426  BIASDOUT(D_REGION_MATCHER_COLORNCC,"-----------------------------");
427  BIASDOUT(D_REGION_MATCHER_COLORNCC,
428  "N*S11:"<<N*S11<<" S1S1:"<< S1S1<<" N*S22:"<<N*S22<<" S2S2:"<<S2S2);
429  double n1 = (N * S11 - S1S1);
430  double n2 = (N * S22 - S2S2);
431  BIASDOUT(D_REGION_MATCHER_COLORNCC,"n1:" <<n1<< " n2:"<<n2);
432  if ((fabs(n1) <= CNCC_EPSILON) || (fabs(n2) <= CNCC_EPSILON))
433  result = 0.0;
434  else
435  result = (N * S12 - S1S2) / sqrt(n1 * n2);
436 }
437 
438 template<class StorageType>
439 void RegionMatcher::CNCC(const unsigned int p1[2], const unsigned int p2[2],
440  const StorageType **ida1, const StorageType **ida2,
441  const unsigned int halfwinsize, double &result) const {
442 #define CNCC_EPSILON 1E-08
443  int nhalfwinsize = -int(halfwinsize);
444  register float Sum1 = 0, Sh1 = 0, Ss1 = 0;
445  register float Sum2 = 0, Sh2 = 0, Ss2 = 0;
446  register float S11 = 0, S12 = 0, S22 = 0;
447  const StorageType *H1, *H2;
448  register float h1, h2, s1, s2, L1, L2;
449  register unsigned int N = 0;
450  register float angle;
451  for (register int r = nhalfwinsize; r <= (int) halfwinsize; r++)
452  for (register int c = nhalfwinsize; c <= (int) halfwinsize; c++) {
453  H1 = (ida1[(p1[1] + r)] + ((p1[0] + c) * 3));
454  s1 = 2.0f * ((float) *(H1 + 1) - 127.0f);
455  L1 = (float) *(H1 + 2);
456  h1 = 2.0f * ((float) *H1 - 127.0f);
457  H2 = (ida2[(p2[1] + r)] + ((p2[0] + c) * 3));
458  s2 = 2.0f * ((float) *(H2 + 1) - 127.0f);
459  L2 = (float) *(H2 + 2);
460  h2 = 2.0f * ((float) *H2 - 127.0f);
461  BIASDOUT(D_REGION_MATCHER_COLORNCC,
462  "HSL1:"<<(int)(h1)<<" "<<(int)(s1)<<" "<<(int)L1<<
463  " HSL2:"<<(int)(h2)<<" "<<(int)(s2)<<" "<<(int)L2);
464  // convert H and S to euclidean coords, H=0 is red
465  // angle = M_PI * (float)(*H1) / 127.5;
466  // h1 = cos(angle)* (float)(*S1);
467  // s1 = sin(angle)* (float)(*S1);
468  // angle = M_PI * (float)(*H2) / 127.5;
469  // h2 = cos(angle)* (float)(*S2);
470  // s2 = sin(angle)* (float)(*S2);
471  BIASDOUT(D_REGION_MATCHER_COLORNCC,
472  "[h1,s1]:["<<(int)h1<<","<<(int)s1<<
473  "] [h2,s2]:["<<(int)h2<<","<<(int)s2);
474 
475  // angle = h1*h2 + s1*s2;
476  angle = float((int) (h1) * (int) (h2) + (int) (s1) * (int) (s2));
477  if (angle < 0)
478  angle = 0;
479  S12 += angle + L1 * L2;
480  Sh1 += h1;
481  Ss1 += s1;
482  Sum1 += L1;
483  Sh2 += h2;
484  Ss2 += s2;
485  Sum2 += L2;
486  S11 += h1 * h1 + s1 * s1 + L1 * L1;
487  S22 += h2 * h2 + s2 * s2 + L2 * L2;
488  N++;
489  }
490 
491  float S1S2 = Sum1 * Sum2;
492  float S1S1 = Sum1 * Sum1;
493  float S2S2 = Sum2 * Sum2;
494  BIASDOUT(D_REGION_MATCHER_COLORNCC,
495  "N*S12:"<<N*S12<<" S1S2:"<<S1S2<<endl<<
496  " Sh1,Ss1:"<<Sh1<<","<<Ss1<<" Sum1:"<<Sum1<<
497  " Sh2,Ss2:"<<Sh2<<","<<Ss2<<" Sum2:"<<Sum2);
498  BIASDOUT(D_REGION_MATCHER_COLORNCC,"-----------------------------");
499  BIASDOUT(D_REGION_MATCHER_COLORNCC,
500  "N*S11:"<<N*S11<<" S1S1:"<< S1S1<<" N*S22:"<<N*S22<<" S2S2:"<<S2S2);
501  float n1 = (N * S11 - S1S1);
502  float n2 = (N * S22 - S2S2);
503  BIASDOUT(D_REGION_MATCHER_COLORNCC,"n1:" <<n1<< " n2:"<<n2);
504  if ((fabs(n1) <= CNCC_EPSILON) || (fabs(n2) <= CNCC_EPSILON))
505  result = 0.0;
506  else
507  result = (N * S12 - S1S2) / sqrt(n1 * n2);
508 }
509 
510 template<class StorageType>
511 void RegionMatcher::HammingDistance(const unsigned int p1[2],
512  const unsigned int p2[2], const StorageType **ida1,
513  const StorageType **ida2, const unsigned int halfwinsize,
514  const unsigned int censussize, double& result) const {
515  BIASERR("Incoming images should have storage type unsigned char");
516  BIASABORT;
517 }
518 template<> BIASMatcher2D_EXPORT
519 void RegionMatcher::HammingDistance<unsigned char>(const unsigned int p1[2],
520  const unsigned int p2[2], const unsigned char **ida1,
521  const unsigned char **ida2, const unsigned int halfwinsize,
522  const unsigned int censussize, double& result) const {
523  int nhalfwinsize = -(int) halfwinsize;
524  result = 0;
525  register unsigned char cr = 0;
526  for (register int c = nhalfwinsize; c <= (int) halfwinsize; c++) {
527  for (register int r = nhalfwinsize; r <= (int) halfwinsize; r++) {
528  for (register unsigned int channel = 0; channel < censussize; channel++) {
529  hamming_distance(ida1[p1[1]+r][(p1[0]+c)*censussize+channel], ida2[p2[1]+r][(p2[0]+c)*censussize+channel], cr);
530  result += cr;
531  }
532  }
533  }
534 }
535 
536 template<> BIASMatcher2D_EXPORT
537 void RegionMatcher::SADSamplingInsensitive<unsigned char>(
538  const unsigned int p1[2], const unsigned int p2[2],
539  const unsigned char **ida1, const unsigned char **ida2,
540  const unsigned int halfwinsize, double& result) const {
541  int nhalfwinsize = -int(halfwinsize);
542  register double res = 0;
543  unsigned int pw1[2], pw2[2];
544  double value;
545  for (register int c = nhalfwinsize; c <= (int) halfwinsize; c++) {
546  for (register int r = nhalfwinsize; r <= (int) halfwinsize; r++) {
547  pw1[1] = p1[1] + r;
548  pw1[0] = p1[0] + c;
549  pw2[1] = p2[1] + r;
550  pw2[0] = p2[0] + c;
551  SamplingInsensitiveDissimilarity(pw1, pw2, ida1, ida2, value);
552  res += value;
553  }
554  }
555  result = (double) res;
556 }
557 
558 template<class StorageType>
560  const unsigned int p2[2], const StorageType **ida1,
561  const StorageType **ida2, double &result) const {
562  float iminus, iplus, r1, r2;
563  register unsigned char c1, c2, c3;//, c4;
564  c1 = ida1[p1[1]][p1[0]];
565  c2 = ida2[p2[1]][p2[0]];
566  c3 = ida2[p2[1]][p2[0] - 1];
567  //c4 = ida2[p2[1]][p2[0] + 1];
568  iminus = 0.5f * (c2 + c3);
569  iplus = 0.5f * (c2 + c3);
570  r1 = c1 - iplus;
571  r2 = iminus - c1;
572  if (r1 < r2)
573  r1 = r2;
574  float dl = (r1 >= 0) ? r1 : 0;
575 
576  c1 = ida2[p2[1]][p2[0]];
577  c2 = ida1[p1[1]][p1[0]];
578  c3 = ida1[p1[1]][p1[0] - 1];
579  //c4 = ida1[p1[1]][p1[0] + 1];
580  iminus = 0.5f * (c2 + c3);
581  iplus = 0.5f * (c2 + c3);
582  r1 = c1 - iplus;
583  r2 = iminus - c1;
584  if (r1 < r2)
585  r1 = r2;
586  float dr = (r1 >= 0) ? r1 : 0;
587 
588  result = (dl < dr) ? dl : dr;
589 }
590 
591 } // end namespace
void CNCC(const unsigned int p1[2], const unsigned int p2[2], const StorageType **ida1, const StorageType **ida2, const unsigned int halfwinsize, double &result) const
ColorNCC for interleaved hsL images.
void SADN(const unsigned int p1[2], const unsigned int p2[2], const StorageType **ida1, const StorageType **ida2, const unsigned int halfwinsize, double &result) const
'normalized' version of SAD, returns 1.0 incase of no differences 0.0 in case of max difference speci...
void ColorNCC(const unsigned int p1[2], const unsigned int p2[2], const StorageType **ida1, const StorageType **ida2, const unsigned int halfwinsize, double &result) const
ColorNCC for interleaved HSV images.
void ColorCCV(const unsigned int p1[2], const unsigned int p2[2], const StorageType **ida1, const StorageType **ida2, const unsigned int halfwinsize, double &result) const
computes the Cross-CoVariance ,which is equal to NCC but does no variance(contrast) normailzation ...
void SSD(const unsigned int p1[2], const unsigned int p2[2], const StorageType **ida1, const StorageType **ida2, const unsigned int halfwinsize, double &result) const
fast sum of square differences (SSD) sim.
void SAD(const unsigned int p1[2], const unsigned int p2[2], const StorageType **ida1, const StorageType **ida2, const unsigned int halfwinsize, double &result) const
fast sum of absolut differences (SAD) sim.
void SamplingInsensitiveDissimilarity(const unsigned int p1[2], const unsigned int p2[2], const StorageType **ida1, const StorageType **ida2, double &result) const
calculates sampling insensitive dissimilarity measure (birchfield, tomasi)
void NCC(const unsigned int p1[2], const unsigned int p2[2], const StorageType **ida1, const StorageType **ida2, const unsigned int halfwinsize, double &result) const
fast ncc between p1 and p2 without boundary checking Faster specialisation for unsigned char exists b...
void HammingDistance(const unsigned int p1[2], const unsigned int p2[2], const StorageType **ida1, const StorageType **ida2, const unsigned int halfwinsize, const unsigned int censussize, double &result) const
fast hamming distance between p1 and p2 without boundary checking.