00001 /*********************************************************************** 00002 * Software License Agreement (BSD License) 00003 * 00004 * Copyright 2008-2009 Marius Muja (mariusm@cs.ubc.ca). All rights reserved. 00005 * Copyright 2008-2009 David G. Lowe (lowe@cs.ubc.ca). All rights reserved. 00006 * 00007 * THE BSD LICENSE 00008 * 00009 * Redistribution and use in source and binary forms, with or without 00010 * modification, are permitted provided that the following conditions 00011 * are met: 00012 * 00013 * 1. Redistributions of source code must retain the above copyright 00014 * notice, this list of conditions and the following disclaimer. 00015 * 2. Redistributions in binary form must reproduce the above copyright 00016 * notice, this list of conditions and the following disclaimer in the 00017 * documentation and/or other materials provided with the distribution. 00018 * 00019 * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR 00020 * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES 00021 * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. 00022 * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT, 00023 * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT 00024 * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, 00025 * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY 00026 * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT 00027 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF 00028 * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 00029 *************************************************************************/ 00030 00031 #ifndef _OPENCV_DIST_H_ 00032 #define _OPENCV_DIST_H_ 00033 00034 #include <cmath> 00035 00036 #include "opencv2/flann/general.h" 00037 00038 namespace cvflann 00039 { 00040 00046 #define flann_dist custom_dist 00047 //#define flann_dist euclidean_dist 00048 00049 00059 template <typename Iterator1, typename Iterator2> 00060 double euclidean_dist(Iterator1 first1, Iterator1 last1, Iterator2 first2, double acc = 0) 00061 { 00062 double distsq = acc; 00063 double diff0, diff1, diff2, diff3; 00064 Iterator1 lastgroup = last1 - 3; 00065 00066 /* Process 4 items with each loop for efficiency. */ 00067 while (first1 < lastgroup) { 00068 diff0 = first1[0] - first2[0]; 00069 diff1 = first1[1] - first2[1]; 00070 diff2 = first1[2] - first2[2]; 00071 diff3 = first1[3] - first2[3]; 00072 distsq += diff0 * diff0 + diff1 * diff1 + diff2 * diff2 + diff3 * diff3; 00073 first1 += 4; 00074 first2 += 4; 00075 } 00076 /* Process last 0-3 pixels. Not needed for standard vector lengths. */ 00077 while (first1 < last1) { 00078 diff0 = *first1++ - *first2++; 00079 distsq += diff0 * diff0; 00080 } 00081 return distsq; 00082 } 00083 00084 CV_EXPORTS double euclidean_dist(const unsigned char* first1, const unsigned char* last1, unsigned char* first2, double acc); 00085 00086 00093 template <typename Iterator1, typename Iterator2> 00094 double manhattan_dist(Iterator1 first1, Iterator1 last1, Iterator2 first2, double acc = 0) 00095 { 00096 double distsq = acc; 00097 double diff0, diff1, diff2, diff3; 00098 Iterator1 lastgroup = last1 - 3; 00099 00100 /* Process 4 items with each loop for efficiency. */ 00101 while (first1 < lastgroup) { 00102 diff0 = fabs(first1[0] - first2[0]); 00103 diff1 = fabs(first1[1] - first2[1]); 00104 diff2 = fabs(first1[2] - first2[2]); 00105 diff3 = fabs(first1[3] - first2[3]); 00106 distsq += diff0 + diff1 + diff2 + diff3; 00107 first1 += 4; 00108 first2 += 4; 00109 } 00110 /* Process last 0-3 pixels. Not needed for standard vector lengths. */ 00111 while (first1 < last1) { 00112 diff0 = fabs(*first1++ - *first2++); 00113 distsq += diff0; 00114 } 00115 return distsq; 00116 } 00117 00118 00119 CV_EXPORTS int flann_minkowski_order(); 00129 template <typename Iterator1, typename Iterator2> 00130 double minkowski_dist(Iterator1 first1, Iterator1 last1, Iterator2 first2, double acc = 0) 00131 { 00132 double distsq = acc; 00133 double diff0, diff1, diff2, diff3; 00134 Iterator1 lastgroup = last1 - 3; 00135 00136 int p = flann_minkowski_order(); 00137 00138 /* Process 4 items with each loop for efficiency. */ 00139 while (first1 < lastgroup) { 00140 diff0 = fabs(first1[0] - first2[0]); 00141 diff1 = fabs(first1[1] - first2[1]); 00142 diff2 = fabs(first1[2] - first2[2]); 00143 diff3 = fabs(first1[3] - first2[3]); 00144 distsq += pow(diff0,p) + pow(diff1,p) + pow(diff2,p) + pow(diff3,p); 00145 first1 += 4; 00146 first2 += 4; 00147 } 00148 /* Process last 0-3 pixels. Not needed for standard vector lengths. */ 00149 while (first1 < last1) { 00150 diff0 = fabs(*first1++ - *first2++); 00151 distsq += pow(diff0,p); 00152 } 00153 return distsq; 00154 } 00155 00156 00157 // L_infinity distance (NOT A VALID KD-TREE DISTANCE - NOT DIMENSIONWISE ADDITIVE) 00158 template <typename Iterator1, typename Iterator2> 00159 double max_dist(Iterator1 first1, Iterator1 last1, Iterator2 first2, double acc = 0) 00160 { 00161 double dist = acc; 00162 Iterator1 lastgroup = last1 - 3; 00163 double diff0, diff1, diff2, diff3; 00164 00165 /* Process 4 items with each loop for efficiency. */ 00166 while (first1 < lastgroup) { 00167 diff0 = fabs(first1[0] - first2[0]); 00168 diff1 = fabs(first1[1] - first2[1]); 00169 diff2 = fabs(first1[2] - first2[2]); 00170 diff3 = fabs(first1[3] - first2[3]); 00171 if (diff0 > dist) dist = diff0; 00172 if (diff1 > dist) dist = diff1; 00173 if (diff2 > dist) dist = diff2; 00174 if (diff3 > dist) dist = diff3; 00175 first1 += 4; 00176 first2 += 4; 00177 } 00178 /* Process last 0-3 pixels. Not needed for standard vector lengths. */ 00179 while (first1 < last1) { 00180 diff0 = fabs(*first1++ - *first2++); 00181 dist = (diff0 > dist) ? diff0 : dist; 00182 } 00183 return dist; 00184 } 00185 00186 00187 template <typename Iterator1, typename Iterator2> 00188 double hist_intersection_kernel(Iterator1 first1, Iterator1 last1, Iterator2 first2) 00189 { 00190 double kernel = 0; 00191 Iterator1 lastgroup = last1 - 3; 00192 double min0, min1, min2, min3; 00193 00194 /* Process 4 items with each loop for efficiency. */ 00195 while (first1 < lastgroup) { 00196 min0 = first1[0] < first2[0] ? first1[0] : first2[0]; 00197 min1 = first1[1] < first2[1] ? first1[1] : first2[1]; 00198 min2 = first1[2] < first2[2] ? first1[2] : first2[2]; 00199 min3 = first1[3] < first2[3] ? first1[3] : first2[3]; 00200 kernel += min0 + min1 + min2 + min3; 00201 first1 += 4; 00202 first2 += 4; 00203 } 00204 /* Process last 0-3 pixels. Not needed for standard vector lengths. */ 00205 while (first1 < last1) { 00206 min0 = first1[0] < first2[0] ? first1[0] : first2[0]; 00207 kernel += min0; 00208 first1++; 00209 first2++; 00210 } 00211 return kernel; 00212 } 00213 00214 template <typename Iterator1, typename Iterator2> 00215 double hist_intersection_dist_sq(Iterator1 first1, Iterator1 last1, Iterator2 first2, double acc = 0) 00216 { 00217 double dist_sq = acc - 2 * hist_intersection_kernel(first1, last1, first2); 00218 while (first1 < last1) { 00219 dist_sq += *first1 + *first2; 00220 first1++; 00221 first2++; 00222 } 00223 return dist_sq; 00224 } 00225 00226 00227 // Hellinger distance 00228 template <typename Iterator1, typename Iterator2> 00229 double hellinger_dist(Iterator1 first1, Iterator1 last1, Iterator2 first2, double acc = 0) 00230 { 00231 double distsq = acc; 00232 double diff0, diff1, diff2, diff3; 00233 Iterator1 lastgroup = last1 - 3; 00234 00235 /* Process 4 items with each loop for efficiency. */ 00236 while (first1 < lastgroup) { 00237 diff0 = sqrt(first1[0]) - sqrt(first2[0]); 00238 diff1 = sqrt(first1[1]) - sqrt(first2[1]); 00239 diff2 = sqrt(first1[2]) - sqrt(first2[2]); 00240 diff3 = sqrt(first1[3]) - sqrt(first2[3]); 00241 distsq += diff0 * diff0 + diff1 * diff1 + diff2 * diff2 + diff3 * diff3; 00242 first1 += 4; 00243 first2 += 4; 00244 } 00245 /* Process last 0-3 pixels. Not needed for standard vector lengths. */ 00246 while (first1 < last1) { 00247 diff0 = sqrt(*first1++) - sqrt(*first2++); 00248 distsq += diff0 * diff0; 00249 } 00250 return distsq; 00251 } 00252 00253 00254 // chi-dsquare distance 00255 template <typename Iterator1, typename Iterator2> 00256 double chi_square_dist(Iterator1 first1, Iterator1 last1, Iterator2 first2, double acc = 0) 00257 { 00258 double dist = acc; 00259 00260 while (first1 < last1) { 00261 double sum = *first1 + *first2; 00262 if (sum > 0) { 00263 double diff = *first1 - *first2; 00264 dist += diff * diff / sum; 00265 } 00266 first1++; 00267 first2++; 00268 } 00269 return dist; 00270 } 00271 00272 00273 // Kullback–Leibler divergence (NOT SYMMETRIC) 00274 template <typename Iterator1, typename Iterator2> 00275 double kl_divergence(Iterator1 first1, Iterator1 last1, Iterator2 first2, double acc = 0) 00276 { 00277 double div = acc; 00278 00279 while (first1 < last1) { 00280 if (*first2 != 0) { 00281 double ratio = *first1 / *first2; 00282 if (ratio > 0) { 00283 div += *first1 * log(ratio); 00284 } 00285 } 00286 first1++; 00287 first2++; 00288 } 00289 return div; 00290 } 00291 00292 00293 00294 00295 CV_EXPORTS flann_distance_t flann_distance_type(); 00303 template <typename Iterator1, typename Iterator2> 00304 double custom_dist(Iterator1 first1, Iterator1 last1, Iterator2 first2, double acc = 0) 00305 { 00306 switch (flann_distance_type()) { 00307 case FLANN_DIST_EUCLIDEAN: 00308 return euclidean_dist(first1, last1, first2, acc); 00309 case FLANN_DIST_MANHATTAN: 00310 return manhattan_dist(first1, last1, first2, acc); 00311 case FLANN_DIST_MINKOWSKI: 00312 return minkowski_dist(first1, last1, first2, acc); 00313 case FLANN_DIST_MAX: 00314 return max_dist(first1, last1, first2, acc); 00315 case FLANN_DIST_HIST_INTERSECT: 00316 return hist_intersection_dist_sq(first1, last1, first2, acc); 00317 case FLANN_DIST_HELLINGER: 00318 return hellinger_dist(first1, last1, first2, acc); 00319 case FLANN_DIST_CS: 00320 return chi_square_dist(first1, last1, first2, acc); 00321 case FLANN_DIST_KL: 00322 return kl_divergence(first1, last1, first2, acc); 00323 default: 00324 return euclidean_dist(first1, last1, first2, acc); 00325 } 00326 } 00327 00328 /* 00329 * This is a "zero iterator". It basically behaves like a zero filled 00330 * array to all algorithms that use arrays as iterators (STL style). 00331 * It's useful when there's a need to compute the distance between feature 00332 * and origin it and allows for better compiler optimisation than using a 00333 * zero-filled array. 00334 */ 00335 template <typename T> 00336 struct ZeroIterator { 00337 00338 T operator*() { 00339 return 0; 00340 } 00341 00342 T operator[](int) { 00343 return 0; 00344 } 00345 00346 ZeroIterator<T>& operator ++(int) { 00347 return *this; 00348 } 00349 00350 ZeroIterator<T>& operator+=(int) { 00351 return *this; 00352 } 00353 00354 }; 00355 00356 CV_EXPORTS ZeroIterator<float>& zero(); 00357 00358 } // namespace cvflann 00359 00360 #endif //_OPENCV_DIST_H_