Basic Image AlgorithmS Library 2.8.0

SimilarityMeasure.cpp

00001 /* This file is part of the BIAS library (Basic ImageAlgorithmS).
00002 
00003  Copyright (C) 2003-2009    (see file CONTACTS for details)
00004  Multimediale Systeme der Informationsverarbeitung
00005  Institut fuer Informatik
00006  Christian-Albrechts-Universitaet Kiel
00007 
00008  BIAS is free software; you can redistribute it and/or modify
00009  it under the terms of the GNU Lesser General Public License as published by
00010  the Free Software Foundation; either version 2.1 of the License, or
00011  (at your option) any later version.
00012 
00013  BIAS is distributed in the hope that it will be useful,
00014  but WITHOUT ANY WARRANTY; without even the implied warranty of
00015  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00016  GNU Lesser General Public License for more details.
00017 
00018  You should have received a copy of the GNU Lesser General Public License
00019  along with BIAS; if not, write to the Free Software
00020  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA*/
00021 #include "RegionMatcher.hh"
00022 
00023 #include <float.h>
00024 
00025 namespace BIAS {
00026 using namespace std;
00027 
00028 template<class StorageType>
00029 void RegionMatcher::SAD(const unsigned int p1[2], const unsigned int p2[2],
00030         const StorageType **ida1, const StorageType **ida2,
00031         const unsigned int halfwinsize, double& result) const {
00032     int nhalfwinsize = -int(halfwinsize);
00033     result = 0.0;
00034     register StorageType st1, st2;
00035     for (register int c = nhalfwinsize; c <= (int) halfwinsize; c++) {
00036         for (register int r = nhalfwinsize; r <= (int) halfwinsize; r++) {
00037             st1 = ida1[p1[1] + r][p1[0] + c];
00038             st2 = ida2[p2[1] + r][p2[0] + c];
00039             result += (double) ((st1 >= st2) ? st1 - st2 : st2 - st1);
00040         }
00041     }
00042 }
00043 
00044 template<> BIASMatcher2D_EXPORT
00045 void RegionMatcher::SAD<unsigned char>(const unsigned int p1[2],
00046         const unsigned int p2[2], const unsigned char **ida1,
00047         const unsigned char **ida2, const unsigned int halfwinsize,
00048         double& result) const {
00049     int nhalfwinsize = -int(halfwinsize);
00050     register long unsigned int res = 0;
00051     register unsigned char c1, c2;
00052     for (register int c = nhalfwinsize; c <= (int) halfwinsize; c++) {
00053         for (register int r = nhalfwinsize; r <= (int) halfwinsize; r++) {
00054             c1 = ida1[p1[1] + r][p1[0] + c];
00055             c2 = ida2[p2[1] + r][p2[0] + c];
00056             res += (c1 >= c2) ? c1 - c2 : c2 - c1;
00057         }
00058     }
00059     result = (double) res;
00060 }
00061 
00062 template<class StorageType>
00063 void RegionMatcher::SAD(const unsigned int p1[2], const unsigned int p2[2],
00064         const StorageType **ida1, const StorageType **ida2,
00065         const unsigned halfwinsize, const unsigned channel_count,
00066         double& result) const {
00067     int nhalfwinsize = -int(halfwinsize);
00068     result = 0.0;
00069     register StorageType st1, st2;
00070     register unsigned channel;
00071     for (register int c = nhalfwinsize; c <= (int) halfwinsize; c++) {
00072         for (register int r = nhalfwinsize; r <= (int) halfwinsize; r++) {
00073             for (channel = 0; channel < channel_count; channel++) {
00074                 st1 = ida1[p1[1] + r][(p1[0] + c) * channel_count + channel];
00075                 st2 = ida2[p2[1] + r][(p2[0] + c) * channel_count + channel];
00076                 result += (double) ((st1 >= st2) ? st1 - st2 : st2 - st1);
00077             }
00078         }
00079     }
00080 }
00081 
00082 template<> BIASMatcher2D_EXPORT
00083 void RegionMatcher::SAD<unsigned char>(const unsigned int p1[2],
00084         const unsigned int p2[2], const unsigned char **ida1,
00085         const unsigned char **ida2, const unsigned halfwinsize,
00086         const unsigned channel_count, double& result) const {
00087     int nhalfwinsize = -int(halfwinsize);
00088     result = 0.0;
00089     register unsigned char st1, st2;
00090     register unsigned channel;
00091     for (register int c = nhalfwinsize; c <= (int) halfwinsize; c++) {
00092         for (register int r = nhalfwinsize; r <= (int) halfwinsize; r++) {
00093             for (channel = 0; channel < channel_count; channel++) {
00094                 st1 = ida1[p1[1] + r][(p1[0] + c) * channel_count + channel];
00095                 st2 = ida2[p2[1] + r][(p2[0] + c) * channel_count + channel];
00096                 result += (double) ((st1 >= st2) ? st1 - st2 : st2 - st1);
00097             }
00098         }
00099     }
00100 }
00101 
00102 template<class StorageType>
00103 void RegionMatcher::SADN(const unsigned int p1[2], const unsigned int p2[2],
00104         const StorageType **ida1, const StorageType **ida2,
00105         const unsigned int halfwinsize, double& result) const {
00106     int nhalfwinsize = -int(halfwinsize);
00107     result = 0.0;
00108     register StorageType st1, st2;
00109     for (register int c = nhalfwinsize; c <= (int) halfwinsize; c++) {
00110         for (register int r = nhalfwinsize; r <= (int) halfwinsize; r++) {
00111             st1 = ida1[p1[1] + r][p1[0] + c];
00112             st2 = ida2[p2[1] + r][p2[0] + c];
00113             result += (st1 >= st2) ? st1 - st2 : st2 - st1;
00114         }
00115     }
00116     // STORAGETYPE_MAX should be the biggest poistive number that can be
00117     // expressed with this storage type
00118 #define STORAGETYPE_MAX 256.0
00119 
00120     unsigned tmp = (halfwinsize << 1) + 1;
00121     result /= double(tmp * tmp * STORAGETYPE_MAX);
00122     result = 1.0 - result;
00123 }
00124 
00125 template<> BIASMatcher2D_EXPORT
00126 void RegionMatcher::SADN<unsigned char>(const unsigned int p1[2],
00127         const unsigned int p2[2], const unsigned char **ida1,
00128         const unsigned char **ida2, const unsigned int halfwinsize,
00129         double& result) const {
00130     int nhalfwinsize = -int(halfwinsize);
00131     register long unsigned int res = 0;
00132     register unsigned char c1, c2;
00133     for (register int c = nhalfwinsize; c <= (int) halfwinsize; c++) {
00134         for (register int r = nhalfwinsize; r <= (int) halfwinsize; r++) {
00135             c1 = ida1[p1[1] + r][p1[0] + c];
00136             c2 = ida2[p2[1] + r][p2[0] + c];
00137             res += (c1 >= c2) ? c1 - c2 : c2 - c1;
00138         }
00139     }
00140     // result = 1.0 - res /(UCHAR_MAX*windowarea) = 1 - res / tmp
00141     // tmp = sqrt(UCHAR_MAX)*(2*halfwinsize+1) ca. 16 * (2 * halfwinsize + 1)
00142     //     = 32 * halfwinsize + 16
00143     unsigned tmp = (halfwinsize << 5) + 16;
00144     result = 1.0 - (double) res / double(tmp * tmp);
00145 }
00146 
00147 template<class StorageType>
00148 void RegionMatcher::SSD(const unsigned int p1[2], const unsigned int p2[2],
00149         const StorageType **ida1, const StorageType **ida2,
00150         const unsigned int halfwinsize, double& result) const {
00151     int nhalfwinsize = -int(halfwinsize);
00152     result = 0.0;
00153         double sd;
00154     for (register int c = nhalfwinsize; c <= (int) halfwinsize; c++) {
00155         for (register int r = nhalfwinsize; r <= (int) halfwinsize; r++) {
00156                   //DONT USE POW(.,2.0)!!!!
00157                   // result += pow(
00158                   //        (double) (ida1[p1[1] + r][p1[0] + c]) -(double) (ida2[p2[1]
00159                   //                + r][p2[0] + c]), 2.0);
00160                   sd = double(ida1[p1[1] + r][p1[0] + c])
00161                        - double(ida2[p2[1] + r][p2[0] + c]);
00162                   sd *= sd;
00163                   result += sd;
00164         }
00165     }
00166 }
00167 
00168 template<> BIASMatcher2D_EXPORT
00169 void RegionMatcher::SSD<unsigned char>(const unsigned int p1[2],
00170         const unsigned int p2[2], const unsigned char **ida1,
00171         const unsigned char **ida2, const unsigned int halfwinsize,
00172         double& result) const {
00173     int nhalfwinsize = -int(halfwinsize);
00174     register long unsigned int res = 0;
00175     register short diff;
00176     for (register int c = nhalfwinsize; c <= (int) halfwinsize; c++) {
00177         for (register int r = nhalfwinsize; r <= (int) halfwinsize; r++) {
00178             diff = (short) (ida1[p1[1] + r][p1[0] + c]) -(short) (ida2[p2[1]
00179                     + r][p2[0] + c]);
00180             res += (unsigned int) (diff * diff);
00181         }
00182     }
00183     result = (double) res;
00184 }
00185 
00186 template<class StorageType>
00187 void RegionMatcher::NCC(const unsigned int p1[2], const unsigned int p2[2],
00188         const StorageType **ida1, const StorageType **ida2,
00189         const unsigned int halfwinsize, double& result) const {
00190     int nhalfwinsize = -int(halfwinsize);
00191     double im1 = 0, im2 = 0, S1 = 0, S2 = 0;
00192     double S11 = 0, S12 = 0, S22 = 0, N = 0;
00193     BIASASSERT(ida1!=NULL);
00194     BIASASSERT(ida2!=NULL);
00195     for (register int r = nhalfwinsize; r <= (int) halfwinsize; r++)
00196         for (register int c = nhalfwinsize; c <= (int) halfwinsize; c++) {
00197             im1 = (double) (ida1[p1[1] + r][p1[0] + c]);
00198             im2 = (double) (ida2[p2[1] + r][p2[0] + c]);
00199             N++;
00200             S1 += im1;
00201             S2 += im2;
00202             S11 += im1 * im1;
00203             S12 += im1 * im2;
00204             S22 += im2 * im2;
00205         }
00206     result = (N * S12 - S1 * S2) / sqrt((N * S11 - S1 * S1) * (N * S22 - S2
00207             * S2));
00208     if (result != result)
00209         result = 0.0;
00210 }
00211 
00212 // faster specialisation because of integer arithmetic
00213 // does work without overflow with halfwindowsizes <= 127
00214 // if sizeof(unsigned char)=1 and sizeof(unsigned long int)=4
00215 template<> BIASMatcher2D_EXPORT void RegionMatcher::NCC<unsigned char>(
00216         const unsigned int p1[2], const unsigned int p2[2],
00217         const unsigned char **ida1, const unsigned char **ida2,
00218         const unsigned int halfwinsize, double& result) const {
00219     int nhalfwinsize = -int(halfwinsize);
00220     register int im1 = 0, im2 = 0;
00221     register unsigned long int S11 = 0, S12 = 0, S22 = 0, N = 0, S1 = 0, S2 = 0;
00222     for (register int r = nhalfwinsize; r <= (int) halfwinsize; r++)
00223         for (register int c = nhalfwinsize; c <= (int) halfwinsize; c++) {
00224             im1 = (int) (ida1[p1[1] + r][p1[0] + c]);
00225             im2 = (int) (ida2[p2[1] + r][p2[0] + c]);
00226             //        if (im1!=im2)
00227             //          cerr << "im1 ("<<p1[0]+c<<", "<<p1[1]+r<<") : "<<im1
00228             //               << "\tim2 ("<<p2[0]+c<<", "<<p2[1]+r<<") : "<<im2<<endl;
00229             N += 1;
00230             S1 += im1;
00231             S2 += im2;
00232             S11 += im1 * im1;
00233             S12 += im1 * im2;
00234             S22 += im2 * im2;
00235         }
00236     //     cerr << " N " << N << "   S1 " << S1 << "   S11 " << S11 << "   S2 " << S2
00237     //          << "   S22 " << S22 << "   S12 " << S12 << endl;
00238     result = ((double) N * (double) S12 - (double) S1 * (double) S2) / sqrt(
00239             ((double) N * (double) S11 - (double) S1 * (double) S1)
00240                     * ((double) N * (double) S22 - (double) S2 * (double) S2));
00241     if (result != result)
00242         result = 0.0;
00243 
00244     //   cout<<result<<" at ("<<p1[0]<<", "<<p1[1]<<") & ("
00245     //       <<p2[0]<<", "<<p2[1]<<")"<<endl;
00246 }
00247 
00248 template<class StorageType>
00249 void RegionMatcher::NCC(const double p1[2], const double p2[2],
00250         const StorageType **ida1, const StorageType **ida2, const unsigned hws,
00251         double& result) {
00252     double im1 = 0, im2 = 0, S1 = 0, S2 = 0;
00253     double S11 = 0, S12 = 0, S22 = 0, N = 0;
00254     if (_ncchws != hws) {
00255         _ncchws = hws;
00256         _nccwinsize = 1 + (hws << 1);
00257         _iim1.newsize(_nccwinsize, _nccwinsize);
00258         _iim2.newsize(_nccwinsize, _nccwinsize);
00259     }
00260     _BilinearRegion(ida1, p1[0], p1[1], _ncchws, _iim1);
00261     _BilinearRegion(ida2, p2[0], p2[1], _ncchws, _iim2);
00262     for (register unsigned r = 0; r <= hws; r++)
00263         for (register unsigned c = 0; c <= hws; c++) {
00264             im1 = (double) _iim1[r][c];
00265             im2 = (double) _iim2[r][c];
00266             N++;
00267             S1 += im1;
00268             S2 += im2;
00269             S11 += im1 * im1;
00270             S12 += im1 * im2;
00271             S22 += im2 * im2;
00272         }
00273     result = (N * S12 - S1 * S2) / sqrt((N * S11 - S1 * S1) * (N * S22 - S2
00274             * S2));
00275     if (result != result)
00276         result = 0.0;
00277 }
00278 
00279 template<class StorageType>
00280 void RegionMatcher::NCC(const unsigned p1[2], const unsigned p2[2],
00281         const StorageType **ida1, const StorageType **roi1,
00282         const StorageType **ida2, const StorageType **roi2, const unsigned hww,
00283         const unsigned hwh, double& result) const {
00284     int nhww = -int(hww);
00285     int nhwh = -int(hwh);
00286     double im1 = 0, im2 = 0, S1 = 0, S2 = 0;
00287     double S11 = 0, S12 = 0, S22 = 0, N = 0;
00288     for (register int r = nhwh; r < (int) hwh; r++) {
00289         for (register int c = nhww; c < (int) hww; c++) {
00290             if (roi1[p1[1] + r][p1[0] + c] != 0 && roi2[p2[1] + r][p2[0] + c]
00291                     != 0) {
00292                 im1 = (double) (ida1[p1[1] + r][p1[0] + c]);
00293                 im2 = (double) (ida2[p2[1] + r][p2[0] + c]);
00294                 N++;
00295                 S1 += im1;
00296                 S2 += im2;
00297                 S11 += im1 * im1;
00298                 S12 += im1 * im2;
00299                 S22 += im2 * im2;
00300             }
00301         }
00302     }
00303     result = (N * S12 - S1 * S2) / sqrt((N * S11 - S1 * S1) * (N * S22 - S2
00304             * S2));
00305     if (result != result)
00306         result = 0.0;
00307 }
00308 
00309 template<class StorageType>
00310 void RegionMatcher::ColorCCV(const unsigned int p1[2],
00311         const unsigned int p2[2], const StorageType **ida1,
00312         const StorageType **ida2, const unsigned int halfwinsize,
00313         double &result) const {
00314     double Sum1 = 0, Sh1 = 0, Ss1 = 0;
00315     double Sum2 = 0, Sh2 = 0, Ss2 = 0;
00316     double S12 = 0, N = 0;
00317     StorageType *H1, *S1, *H2, *S2;
00318     double I1, I2;
00319     double angle, h1, s1, h2, s2;
00320     for (register int r = -halfwinsize; r <= (int) halfwinsize; r++)
00321         for (register int c = -halfwinsize; c <= (int) halfwinsize; c++) {
00322             H1 = (ida1[(p1[1] + r)] + ((p1[0] + c) * 3));
00323             S1 = H1 + 1;
00324             I1 = *(S1 + 1);
00325             H2 = (ida2[(p2[1] + r)] + ((p2[0] + c) * 3));
00326             S2 = H2 + 1;
00327             I2 = *(S2 + 1);
00328             BIASDOUT(D_REGION_MATCHER_COLORCCV,
00329                     "HSI1:"<<(int)(*H1)<<" "<<(int)(*S1)<<" "<<(int)I1<<
00330                     "  HSI2:"<<(int)(*H2)<<" "<<(int)(*S2)<<" "<<(int)I2);
00331             // convert H and S to euclidean coords, H=0 is red
00332             angle = M_PI * (float) (*H1) / 127.5;
00333             h1 = cos(angle) * (*S1);
00334             s1 = sin(angle) * (*S1);
00335             angle = M_PI * (float) (*H2) / 127.5;
00336             h2 = cos(angle) * (*S2);
00337             s2 = sin(angle) * (*S2);
00338             BIASDOUT(D_REGION_MATCHER_COLORCCV,
00339                     "[h1,s1]:["<<(int)h1<<","<<(int)s1<<
00340                     "] [h2,s2]:["<<(int)h2<<","<<(int)s2);
00341             angle = h1 * h2 + s1 * s2;
00342             if (angle < 0)
00343                 angle = 0;
00344             S12 += angle + I1 * I2;
00345             Sh1 += h1;
00346             Ss1 += s1;
00347             Sum1 += I1;
00348             Sh2 += h2;
00349             Ss2 += s2;
00350             Sum2 += I2;
00351             N++;
00352         }
00353     double S1S2 = Sh1 * Sh2 + Ss1 * Ss2;
00354     if (S1S2 < 0)
00355         S1S2 = 0;
00356     S1S2 += Sum1 * Sum2;
00357     BIASDOUT(D_REGION_MATCHER_COLORCCV,
00358             "S12:"<<S12<<" Sh1,Ss1:"<<Sh1<<","<<Ss1<<"  Sum1:"<<Sum1<<
00359             " Sh2,Ss2:"<<Sh2<<","<<Ss2<<"  Sum2:"<<Sum2);
00360     result = (S12 - S1S2 / N);
00361 }
00362 
00363 template<class StorageType>
00364 void RegionMatcher::ColorNCC(const unsigned int p1[2],
00365         const unsigned int p2[2], const StorageType **ida1,
00366         const StorageType **ida2, const unsigned int halfwinsize,
00367         double &result) const {
00368     int nhalfwinsize = -int(halfwinsize);
00369 #define CNCC_EPSILON 1E-08
00370     double Sum1 = 0, Sh1 = 0, Ss1 = 0;
00371     double Sum2 = 0, Sh2 = 0, Ss2 = 0;
00372     double S11 = 0, S12 = 0, S22 = 0, N = 0;
00373     const StorageType *H1, *S1, *H2, *S2;
00374     double I1, I2;
00375     double angle, h1, s1, h2, s2;
00376     for (register int r = nhalfwinsize; r <= (int) halfwinsize; r++)
00377         for (register int c = nhalfwinsize; c <= (int) halfwinsize; c++) {
00378             H1 = (ida1[(p1[1] + r)] + ((p1[0] + c) * 3));
00379             S1 = H1 + 1;
00380             I1 = *(S1 + 1);
00381             H2 = (ida2[(p2[1] + r)] + ((p2[0] + c) * 3));
00382             S2 = H2 + 1;
00383             I2 = *(S2 + 1);
00384             BIASDOUT(D_REGION_MATCHER_COLORNCC,
00385                     "HSI1:"<<(int)(*H1)<<" "<<(int)(*S1)<<" "<<(int)I1<<
00386                     "  HSI2:"<<(int)(*H2)<<" "<<(int)(*S2)<<" "<<(int)I2);
00387             // convert H and S to euclidean coords, H=0 is red
00388             angle = M_PI * (float) (*H1) / 127.5;
00389             h1 = cos(angle) * (float) (*S1);
00390             s1 = sin(angle) * (float) (*S1);
00391             angle = M_PI * (float) (*H2) / 127.5;
00392             h2 = cos(angle) * (float) (*S2);
00393             s2 = sin(angle) * (float) (*S2);
00394             BIASDOUT(D_REGION_MATCHER_COLORNCC,
00395                     "[h1,s1]:["<<(int)h1<<","<<(int)s1<<
00396                     "] [h2,s2]:["<<(int)h2<<","<<(int)s2);
00397             angle = h1 * h2 + s1 * s2; // scalar product
00398             //    angle = fabs(double(*H1)-double(*H2));
00399             //       if (angle>180)
00400             //         angle=255-angle;
00401             //       angle =  cos(angle/255.0) * double(*S1)*double(*S2); // scalar product
00402             if (angle < 0)
00403                 angle = 0;
00404             S12 += angle + I1 * I2;
00405             Sh1 += h1;
00406             Ss1 += s1;
00407             Sum1 += I1;
00408             Sh2 += h2;
00409             Ss2 += s2;
00410             Sum2 += I2;
00411             S11 += h1 * h1 + s1 * s1 + I1 * I1;
00412             S22 += h2 * h2 + s2 * s2 + I2 * I2;
00413             N++;
00414         }
00415 
00416     double S1S2 = Sh1 * Sh2 + Ss1 * Ss2;
00417     if (S1S2 < 0)
00418         S1S2 = 0;
00419     S1S2 = Sum1 * Sum2;
00420     double S1S1 = Sum1 * Sum1;
00421     double S2S2 = Sum2 * Sum2;
00422     BIASDOUT(D_REGION_MATCHER_COLORNCC,
00423             "N*S12:"<<N*S12<<"  S1S2:"<<S1S2<<endl<<
00424             " Sh1,Ss1:"<<Sh1<<","<<Ss1<<"  Sum1:"<<Sum1<<
00425             " Sh2,Ss2:"<<Sh2<<","<<Ss2<<"  Sum2:"<<Sum2);
00426     BIASDOUT(D_REGION_MATCHER_COLORNCC,"-----------------------------");
00427     BIASDOUT(D_REGION_MATCHER_COLORNCC,
00428             "N*S11:"<<N*S11<<" S1S1:"<< S1S1<<" N*S22:"<<N*S22<<" S2S2:"<<S2S2);
00429     double n1 = (N * S11 - S1S1);
00430     double n2 = (N * S22 - S2S2);
00431     BIASDOUT(D_REGION_MATCHER_COLORNCC,"n1:" <<n1<< " n2:"<<n2);
00432     if ((fabs(n1) <= CNCC_EPSILON) || (fabs(n2) <= CNCC_EPSILON))
00433         result = 0.0;
00434     else
00435         result = (N * S12 - S1S2) / sqrt(n1 * n2);
00436 }
00437 
00438 template<class StorageType>
00439 void RegionMatcher::CNCC(const unsigned int p1[2], const unsigned int p2[2],
00440         const StorageType **ida1, const StorageType **ida2,
00441         const unsigned int halfwinsize, double &result) const {
00442 #define CNCC_EPSILON 1E-08
00443     int nhalfwinsize = -int(halfwinsize);
00444     register float Sum1 = 0, Sh1 = 0, Ss1 = 0;
00445     register float Sum2 = 0, Sh2 = 0, Ss2 = 0;
00446     register float S11 = 0, S12 = 0, S22 = 0;
00447     const StorageType *H1, *H2;
00448     register float h1, h2, s1, s2, L1, L2;
00449     register unsigned int N = 0;
00450     register float angle;
00451     for (register int r = nhalfwinsize; r <= (int) halfwinsize; r++)
00452         for (register int c = nhalfwinsize; c <= (int) halfwinsize; c++) {
00453             H1 = (ida1[(p1[1] + r)] + ((p1[0] + c) * 3));
00454             s1 = 2.0f * ((float) *(H1 + 1) - 127.0f);
00455             L1 = (float) *(H1 + 2);
00456             h1 = 2.0f * ((float) *H1 - 127.0f);
00457             H2 = (ida2[(p2[1] + r)] + ((p2[0] + c) * 3));
00458             s2 = 2.0f * ((float) *(H2 + 1) - 127.0f);
00459             L2 = (float) *(H2 + 2);
00460             h2 = 2.0f * ((float) *H2 - 127.0f);
00461             BIASDOUT(D_REGION_MATCHER_COLORNCC,
00462                     "HSL1:"<<(int)(h1)<<" "<<(int)(s1)<<" "<<(int)L1<<
00463                     "  HSL2:"<<(int)(h2)<<" "<<(int)(s2)<<" "<<(int)L2);
00464             // convert H and S to euclidean coords, H=0 is red
00465             //   angle = M_PI * (float)(*H1) / 127.5;
00466             //       h1 = cos(angle)* (float)(*S1);
00467             //       s1 = sin(angle)* (float)(*S1);
00468             //       angle = M_PI * (float)(*H2) / 127.5;
00469             //       h2 = cos(angle)* (float)(*S2);
00470             //       s2 = sin(angle)* (float)(*S2);
00471             BIASDOUT(D_REGION_MATCHER_COLORNCC,
00472                     "[h1,s1]:["<<(int)h1<<","<<(int)s1<<
00473                     "] [h2,s2]:["<<(int)h2<<","<<(int)s2);
00474 
00475             //      angle =  h1*h2 + s1*s2;
00476             angle = float((int) (h1) * (int) (h2) + (int) (s1) * (int) (s2));
00477             if (angle < 0)
00478                 angle = 0;
00479             S12 += angle + L1 * L2;
00480             Sh1 += h1;
00481             Ss1 += s1;
00482             Sum1 += L1;
00483             Sh2 += h2;
00484             Ss2 += s2;
00485             Sum2 += L2;
00486             S11 += h1 * h1 + s1 * s1 + L1 * L1;
00487             S22 += h2 * h2 + s2 * s2 + L2 * L2;
00488             N++;
00489         }
00490 
00491     float S1S2 = Sum1 * Sum2;
00492     float S1S1 = Sum1 * Sum1;
00493     float S2S2 = Sum2 * Sum2;
00494     BIASDOUT(D_REGION_MATCHER_COLORNCC,
00495             "N*S12:"<<N*S12<<"  S1S2:"<<S1S2<<endl<<
00496             " Sh1,Ss1:"<<Sh1<<","<<Ss1<<"  Sum1:"<<Sum1<<
00497             " Sh2,Ss2:"<<Sh2<<","<<Ss2<<"  Sum2:"<<Sum2);
00498     BIASDOUT(D_REGION_MATCHER_COLORNCC,"-----------------------------");
00499     BIASDOUT(D_REGION_MATCHER_COLORNCC,
00500             "N*S11:"<<N*S11<<" S1S1:"<< S1S1<<" N*S22:"<<N*S22<<" S2S2:"<<S2S2);
00501     float n1 = (N * S11 - S1S1);
00502     float n2 = (N * S22 - S2S2);
00503     BIASDOUT(D_REGION_MATCHER_COLORNCC,"n1:" <<n1<< " n2:"<<n2);
00504     if ((fabs(n1) <= CNCC_EPSILON) || (fabs(n2) <= CNCC_EPSILON))
00505         result = 0.0;
00506     else
00507         result = (N * S12 - S1S2) / sqrt(n1 * n2);
00508 }
00509 
00510 template<class StorageType>
00511 void RegionMatcher::HammingDistance(const unsigned int p1[2],
00512         const unsigned int p2[2], const StorageType **ida1,
00513         const StorageType **ida2, const unsigned int halfwinsize,
00514         const unsigned int censussize, double& result) const {
00515     BIASERR("Incoming images should have storage type unsigned char");
00516     BIASABORT;
00517 }
00518 template<> BIASMatcher2D_EXPORT
00519 void RegionMatcher::HammingDistance<unsigned char>(const unsigned int p1[2],
00520         const unsigned int p2[2], const unsigned char **ida1,
00521         const unsigned char **ida2, const unsigned int halfwinsize,
00522         const unsigned int censussize, double& result) const {
00523     int nhalfwinsize = -(int) halfwinsize;
00524     result = 0;
00525     register unsigned char cr = 0;
00526     for (register int c = nhalfwinsize; c <= (int) halfwinsize; c++) {
00527         for (register int r = nhalfwinsize; r <= (int) halfwinsize; r++) {
00528             for (register unsigned int channel = 0; channel < censussize; channel++) {
00529                 hamming_distance(ida1[p1[1]+r][(p1[0]+c)*censussize+channel], ida2[p2[1]+r][(p2[0]+c)*censussize+channel], cr);
00530                 result += cr;
00531             }
00532         }
00533     }
00534 }
00535 
00536 template<> BIASMatcher2D_EXPORT
00537 void RegionMatcher::SADSamplingInsensitive<unsigned char>(
00538         const unsigned int p1[2], const unsigned int p2[2],
00539         const unsigned char **ida1, const unsigned char **ida2,
00540         const unsigned int halfwinsize, double& result) const {
00541     int nhalfwinsize = -int(halfwinsize);
00542     register double res = 0;
00543     unsigned int pw1[2], pw2[2];
00544     double value;
00545     for (register int c = nhalfwinsize; c <= (int) halfwinsize; c++) {
00546         for (register int r = nhalfwinsize; r <= (int) halfwinsize; r++) {
00547             pw1[1] = p1[1] + r;
00548             pw1[0] = p1[0] + c;
00549             pw2[1] = p2[1] + r;
00550             pw2[0] = p2[0] + c;
00551             SamplingInsensitiveDissimilarity(pw1, pw2, ida1, ida2, value);
00552             res += value;
00553         }
00554     }
00555     result = (double) res;
00556 }
00557 
00558 template<class StorageType>
00559 void RegionMatcher::SamplingInsensitiveDissimilarity(const unsigned int p1[2],
00560         const unsigned int p2[2], const StorageType **ida1,
00561         const StorageType **ida2, double &result) const {
00562     float iminus, iplus, r1, r2;
00563         register unsigned char c1, c2, c3;//, c4;
00564     c1 = ida1[p1[1]][p1[0]];
00565     c2 = ida2[p2[1]][p2[0]];
00566     c3 = ida2[p2[1]][p2[0] - 1];
00567         //c4 = ida2[p2[1]][p2[0] + 1];
00568     iminus = 0.5f * (c2 + c3);
00569     iplus = 0.5f * (c2 + c3);
00570     r1 = c1 - iplus;
00571     r2 = iminus - c1;
00572     if (r1 < r2)
00573         r1 = r2;
00574     float dl = (r1 >= 0) ? r1 : 0;
00575 
00576     c1 = ida2[p2[1]][p2[0]];
00577     c2 = ida1[p1[1]][p1[0]];
00578     c3 = ida1[p1[1]][p1[0] - 1];
00579         //c4 = ida1[p1[1]][p1[0] + 1];
00580     iminus = 0.5f * (c2 + c3);
00581     iplus = 0.5f * (c2 + c3);
00582     r1 = c1 - iplus;
00583     r2 = iminus - c1;
00584     if (r1 < r2)
00585         r1 = r2;
00586     float dr = (r1 >= 0) ? r1 : 0;
00587 
00588     result = (dl < dr) ? dl : dr;
00589 }
00590 
00591 } // end namespace
 All Classes Functions Variables Typedefs Enumerations Enumerator Friends