git.maemo.org Git - opencv/blob - 3rdparty/flann/nn/index_testing.cpp

   1 /***********************************************************************
   2  * Software License Agreement (BSD License)
   3  *
   4  * Copyright 2008-2009  Marius Muja (mariusm@cs.ubc.ca). All rights reserved.
   5  * Copyright 2008-2009  David G. Lowe (lowe@cs.ubc.ca). All rights reserved.
   6  *
   7  * THE BSD LICENSE
   8  *
   9  * Redistribution and use in source and binary forms, with or without
  10  * modification, are permitted provided that the following conditions
  11  * are met:
  12  *
  13  * 1. Redistributions of source code must retain the above copyright
  14  *    notice, this list of conditions and the following disclaimer.
  15  * 2. Redistributions in binary form must reproduce the above copyright
  16  *    notice, this list of conditions and the following disclaimer in the
  17  *    documentation and/or other materials provided with the distribution.
  18  *
  19  * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
  20  * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
  21  * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
  22  * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
  23  * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
  24  * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
  25  * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
  26  * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
  27  * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
  28  * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
  29  *************************************************************************/
  30
  31 #include "index_testing.h"
  32 #include "result_set.h"
  33 #include "timer.h"
  34 #include "logger.h"
  35 #include "dist.h"
  36 #include "common.h"
  37
  38 #include <algorithm>
  39 #include <math.h>
  40 #include <string.h>
  41 #include <stdlib.h>
  42
  43 namespace flann
  44 {
  45
  46 const float SEARCH_EPS = 0.001f;
  47
  48 int countCorrectMatches(int* neighbors, int* groundTruth, int n)
  49 {
  50     int count = 0;
  51     for (int i=0;i<n;++i) {
  52         for (int k=0;k<n;++k) {
  53             if (neighbors[i]==groundTruth[k]) {
  54                 count++;
  55                 break;
  56             }
  57         }
  58     }
  59     return count;
  60 }
  61
  62
  63 float computeDistanceRaport(const Matrix<float>& inputData, float* target, int* neighbors, int* groundTruth, int veclen, int n)
  64 {
  65         float* target_end = target + veclen;
  66     float ret = 0;
  67     for (int i=0;i<n;++i) {
  68         float den = (float)flann_dist(target,target_end, inputData[groundTruth[i]]);
  69         float num = (float)flann_dist(target,target_end, inputData[neighbors[i]]);
  70
  71 //        printf("den=%g,num=%g\n",den,num);
  72
  73         if (den==0 && num==0) {
  74             ret += 1;
  75         } else {
  76             ret += num/den;
  77         }
  78     }
  79
  80     return ret;
  81 }
  82
  83 float search_with_ground_truth(NNIndex& index, const Matrix<float>& inputData, const Matrix<float>& testData, const Matrix<int>& matches, int nn, int checks, float& time, float& dist, int skipMatches)
  84 {
  85     if (matches.cols<nn) {
  86         logger.info("matches.cols=%d, nn=%d\n",matches.cols,nn);
  87
  88         throw FLANNException("Ground truth is not computed for as many neighbors as requested");
  89     }
  90
  91     KNNResultSet resultSet(nn+skipMatches);
  92     SearchParams searchParams(checks);
  93
  94     int correct = 0;
  95     float distR = 0;
  96     StartStopTimer t;
  97     int repeats = 0;
  98     while (t.value<0.2) {
  99         repeats++;
 100         t.start();
 101         correct = 0;
 102         distR = 0;
 103         for (int i = 0; i < testData.rows; i++) {
 104             float* target = testData[i];
 105             resultSet.init(target, testData.cols);
 106             index.findNeighbors(resultSet,target, searchParams);
 107             int* neighbors = resultSet.getNeighbors();
 108             neighbors = neighbors+skipMatches;
 109
 110             correct += countCorrectMatches(neighbors,matches[i], nn);
 111             distR += computeDistanceRaport(inputData, target,neighbors,matches[i], testData.cols, nn);
 112         }
 113         t.stop();
 114     }
 115     time = (float)(t.value/repeats);
 116
 117
 118     float precicion = (float)correct/(nn*testData.rows);
 119
 120     dist = distR/(testData.rows*nn);
 121
 122     logger.info("%8d %10.4g %10.5g %10.5g %10.5g\n",
 123             checks, precicion, time, 1000.0 * time / testData.rows, dist);
 124
 125     return precicion;
 126 }
 127
 128 void search_for_neighbors(NNIndex& index, const Matrix<float>& testset, Matrix<int>& result,  Matrix<float>& dists, const SearchParams& searchParams, int skip)
 129 {
 130     assert(testset.rows == result.rows);
 131
 132     int nn = result.cols;
 133     KNNResultSet resultSet(nn+skip);
 134
 135
 136     for (int i = 0; i < testset.rows; i++) {
 137         float* target = testset[i];
 138         resultSet.init(target, testset.cols);
 139
 140         index.findNeighbors(resultSet,target, searchParams);
 141
 142         int* neighbors = resultSet.getNeighbors();
 143         float* distances = resultSet.getDistances();
 144         memcpy(result[i], neighbors+skip, nn*sizeof(int));
 145         memcpy(dists[i], distances+skip, nn*sizeof(float));
 146     }
 147
 148 }
 149
 150 float test_index_checks(NNIndex& index, const Matrix<float>& inputData, const Matrix<float>& testData, const Matrix<int>& matches, int checks, float& precision, int nn, int skipMatches)
 151 {
 152     logger.info("  Nodes  Precision(%)   Time(s)   Time/vec(ms)  Mean dist\n");
 153     logger.info("---------------------------------------------------------\n");
 154
 155     float time = 0;
 156     float dist = 0;
 157     precision = search_with_ground_truth(index, inputData, testData, matches, nn, checks, time, dist, skipMatches);
 158
 159     return time;
 160 }
 161
 162
 163 float test_index_precision(NNIndex& index, const Matrix<float>& inputData, const Matrix<float>& testData, const Matrix<int>& matches,
 164              float precision, int& checks, int nn, int skipMatches)
 165 {
 166     logger.info("  Nodes  Precision(%)   Time(s)   Time/vec(ms)  Mean dist\n");
 167     logger.info("---------------------------------------------------------\n");
 168
 169     int c2 = 1;
 170     float p2;
 171     int c1 = 1;
 172     float p1;
 173     float time;
 174     float dist;
 175
 176     p2 = search_with_ground_truth(index, inputData, testData, matches, nn, c2, time, dist, skipMatches);
 177
 178     if (p2>precision) {
 179         logger.info("Got as close as I can\n");
 180         checks = c2;
 181         return time;
 182     }
 183
 184     while (p2<precision) {
 185         c1 = c2;
 186         p1 = p2;
 187         c2 *=2;
 188         p2 = search_with_ground_truth(index, inputData, testData, matches, nn, c2, time, dist, skipMatches);
 189     }
 190
 191     int cx;
 192     float realPrecision;
 193     if (fabs(p2-precision)>SEARCH_EPS) {
 194         logger.info("Start linear estimation\n");
 195         // after we got to values in the vecinity of the desired precision
 196         // use linear approximation get a better estimation
 197
 198         cx = (c1+c2)/2;
 199         realPrecision = search_with_ground_truth(index, inputData, testData, matches, nn, cx, time, dist, skipMatches);
 200         while (fabs(realPrecision-precision)>SEARCH_EPS) {
 201
 202             if (realPrecision<precision) {
 203                 c1 = cx;
 204             }
 205             else {
 206                 c2 = cx;
 207             }
 208             cx = (c1+c2)/2;
 209             if (cx==c1) {
 210                 logger.info("Got as close as I can\n");
 211                 break;
 212             }
 213             realPrecision = search_with_ground_truth(index, inputData, testData, matches, nn, cx, time, dist, skipMatches);
 214         }
 215
 216         c2 = cx;
 217         p2 = realPrecision;
 218
 219     } else {
 220         logger.info("No need for linear estimation\n");
 221         cx = c2;
 222         realPrecision = p2;
 223     }
 224
 225     checks = cx;
 226     return time;
 227 }
 228
 229
 230 float test_index_precisions(NNIndex& index, const Matrix<float>& inputData, const Matrix<float>& testData, const Matrix<int>& matches,
 231                     float* precisions, int precisions_length, int nn, int skipMatches, float maxTime)
 232 {
 233     // make sure precisions array is sorted
 234     sort(precisions, precisions+precisions_length);
 235
 236     int pindex = 0;
 237     float precision = precisions[pindex];
 238
 239     logger.info("  Nodes  Precision(%)   Time(s)   Time/vec(ms)  Mean dist");
 240     logger.info("---------------------------------------------------------");
 241
 242     int c2 = 1;
 243     float p2;
 244
 245     int c1 = 1;
 246     float p1;
 247
 248     float time;
 249     float dist;
 250
 251     p2 = search_with_ground_truth(index, inputData, testData, matches, nn, c2, time, dist, skipMatches);
 252
 253     // if precision for 1 run down the tree is already
 254     // better then some of the requested precisions, then
 255     // skip those
 256     while (precisions[pindex]<p2 && pindex<precisions_length) {
 257         pindex++;
 258     }
 259
 260     if (pindex==precisions_length) {
 261         logger.info("Got as close as I can\n");
 262         return time;
 263     }
 264
 265     for (int i=pindex;i<precisions_length;++i) {
 266
 267         precision = precisions[i];
 268         while (p2<precision) {
 269             c1 = c2;
 270             p1 = p2;
 271             c2 *=2;
 272             p2 = search_with_ground_truth(index, inputData, testData, matches, nn, c2, time, dist, skipMatches);
 273             if (maxTime> 0 && time > maxTime && p2<precision) return time;
 274         }
 275
 276         int cx;
 277         float realPrecision;
 278         if (fabs(p2-precision)>SEARCH_EPS) {
 279             logger.info("Start linear estimation\n");
 280             // after we got to values in the vecinity of the desired precision
 281             // use linear approximation get a better estimation
 282
 283             cx = (c1+c2)/2;
 284             realPrecision = search_with_ground_truth(index, inputData, testData, matches, nn, cx, time, dist, skipMatches);
 285             while (fabs(realPrecision-precision)>SEARCH_EPS) {
 286
 287                 if (realPrecision<precision) {
 288                     c1 = cx;
 289                 }
 290                 else {
 291                     c2 = cx;
 292                 }
 293                 cx = (c1+c2)/2;
 294                 if (cx==c1) {
 295                     logger.info("Got as close as I can\n");
 296                     break;
 297                 }
 298                 realPrecision = search_with_ground_truth(index, inputData, testData, matches, nn, cx, time, dist, skipMatches);
 299             }
 300
 301             c2 = cx;
 302             p2 = realPrecision;
 303
 304         } else {
 305             logger.info("No need for linear estimation\n");
 306             cx = c2;
 307             realPrecision = p2;
 308         }
 309
 310     }
 311     return time;
 312 }
 313
 314 }