knnshow.cpp
Go to the documentation of this file.
00001 /*
00002 
00003 Copyright (c) 2010--2011, Stephane Magnenat, ASL, ETHZ, Switzerland
00004 You can contact the author at <stephane at magnenat dot net>
00005 
00006 All rights reserved.
00007 
00008 Redistribution and use in source and binary forms, with or without
00009 modification, are permitted provided that the following conditions are met:
00010     * Redistributions of source code must retain the above copyright
00011       notice, this list of conditions and the following disclaimer.
00012     * Redistributions in binary form must reproduce the above copyright
00013       notice, this list of conditions and the following disclaimer in the
00014       documentation and/or other materials provided with the distribution.
00015     * Neither the name of the <organization> nor the
00016       names of its contributors may be used to endorse or promote products
00017       derived from this software without specific prior written permission.
00018 
00019 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
00020 ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
00021 WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
00022 DISCLAIMED. IN NO EVENT SHALL ETH-ASL BE LIABLE FOR ANY
00023 DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
00024 (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
00025 LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
00026 ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
00027 (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00028 SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00029 
00030 */
00031 
00032 #include "nabo/nabo.h"
00033 #include <iostream>
00034 #include <fstream>
00035 #include <stdexcept>
00036 
00037 using namespace std;
00038 using namespace Nabo;
00039 
00040 template<typename T>
00041 typename NearestNeighbourSearch<T>::Matrix load(const char *fileName)
00042 {
00043         typedef typename NearestNeighbourSearch<T>::Matrix Matrix;
00044         
00045         ifstream ifs(fileName);
00046         if (!ifs.good())
00047                 throw runtime_error(string("Cannot open file ") + fileName);
00048         
00049         vector<T> data;
00050         int dim(0);
00051         bool firstLine(true);
00052         
00053         while (!ifs.eof())
00054         {
00055                 char line[1024];
00056                 ifs.getline(line, sizeof(line));
00057                 line[sizeof(line)-1] = 0;
00058                 
00059                 char *token = strtok(line, " \t,;");
00060                 while (token)
00061                 {
00062                         if (firstLine)
00063                                 ++dim;
00064                         data.push_back(atof(token));
00065                         //cout << atof(token) << " ";
00066                         token = strtok(NULL, " \t,;"); // FIXME: non reentrant, use strtok_r
00067                 }
00068                 firstLine = false;
00069         }
00070         
00071         return Matrix::Map(&data[0], dim, data.size() / dim);
00072 }
00073 
00074 template<typename T>
00075 void dumpCoordinateForSVG(const typename NearestNeighbourSearch<T>::Vector coord, const float zoom = 1, const float ptSize = 1, const char* style = "stroke=\"black\" fill=\"red\"")
00076 {
00077         if (coord.size() == 2)
00078                 cout
00079                         << "<circle cx=\"" << zoom*coord(0)
00080                         << "\" cy=\"" << zoom*coord(1)
00081                         << "\" r=\"" << ptSize
00082                         << "\" stroke-width=\"" << 0.2 * ptSize
00083                         << "\" " << style << "/>" << endl;
00084         else
00085                 assert(false);
00086 }
00087 
00088 int main(int argc, char* argv[])
00089 {
00090         typedef Nabo::NearestNeighbourSearch<float>::Matrix Matrix;
00091         typedef Nabo::NearestNeighbourSearch<float>::Vector Vector;
00092         typedef Nabo::NearestNeighbourSearch<float>::Index Index;
00093         typedef Nabo::NearestNeighbourSearch<float>::Indexes Indexes;
00094         typedef Nabo::BruteForceSearch<float> BFSF;
00095         typedef Nabo::KDTree<float> KDTF;
00096         
00097         if (argc != 2)
00098         {
00099                 cerr << "Usage " << argv[0] << " DATA" << endl;
00100                 return 1;
00101         }
00102         
00103         Matrix d(load<float>(argv[1]));
00104         BFSF bfs(d);
00105         KDTF kdt(d);
00106         const Index K(10);
00107         
00108         // uncomment to compare KDTree with brute force search
00109         if (K >= d.size())
00110                 return 2;
00111         const int itCount(100);
00112         for (int i = 0; i < itCount; ++i)
00113         {
00114                 //Vector q(bfs.minBound.size());
00115                 //for (int i = 0; i < q.size(); ++i)
00116                 //      q(i) = bfs.minBound(i) + float(rand()) * (bfs.maxBound(i) - bfs.minBound(i)) / float(RAND_MAX);
00117                 Vector q(d.col(rand() % d.cols()));
00118                 q.cwise() += 0.01;
00119                 //cerr << "bfs:\n";
00120                 Indexes indexes_bf(bfs.knn(q, K, false));
00121                 //cerr << "kdt:\n";
00122                 Indexes indexes_kdtree(kdt.knn(q, K, false));
00123                 //cout << indexes_bf.size() << " " << indexes_kdtree.size() << " " << K << endl;
00124                 if (ndexes_bf.size() != indexes_kdtree.size())
00125                         return 2;
00126                 assert(indexes_bf.size() == indexes_kdtree.size());
00127                 assert(indexes_bf.size() == K);
00128                 //cerr << "\nquery:\n" << q << "\n\n";
00129                 for (size_t j = 0; j < K; ++j)
00130                 {
00131                         Vector pbf(d.col(indexes_bf[j]));
00132                         Vector pkdtree(d.col(indexes_kdtree[j]));
00133                         //cerr << "index " << j << ": " << indexes_bf[j] << ", " << indexes_kdtree[j] << "\n";  
00134                         //cerr << "point " << j << ": " << "\nbf:\n" << pbf << "\nkdtree:\n" << pkdtree << "\n\n";
00135                         assert(indexes_bf[j] == indexes_kdtree[j]);
00136                         assert(dist2(pbf, pkdtree) < numeric_limits<float>::epsilon());
00137                 }
00138         }
00139         cerr << "stats kdtree: "
00140                 << kdt.getStatistics().totalVisitCount << " on "
00141                 << itCount * d.cols() << " ("
00142                 << double(100 * kdt.getStatistics().totalVisitCount) /  double(itCount * d.cols()) << " %"
00143                 << ")" << endl;
00144         
00145         /*
00146         // uncomment to randomly get a point and find its minimum
00147         cout << "<?xml version=\"1.0\" standalone=\"no\"?>" << endl;
00148         cout << "<!DOCTYPE svg PUBLIC \"-//W3C//DTD SVG 1.1//EN\" \"http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd\">" << endl;
00149         cout << "<svg width=\"100%\" height=\"100%\" version=\"1.1\" xmlns=\"http://www.w3.org/2000/svg\">" << endl;
00150         srand(time(0));
00151         for (int i = 0; i < d.cols(); ++i)
00152                 dumpCoordinateForSVG<float>(d.col(i), 100, 1);
00153         Vector q(bfs.minBound.size());
00154         for (int i = 0; i < q.size(); ++i)
00155                 q(i) = bfs.minBound(i) + float(rand()) * (bfs.maxBound(i) - bfs.minBound(i)) / float(RAND_MAX);
00156         Indexes indexes_bf(bfs.knn(q, K, false));
00157         for (size_t i = 0; i < K; ++i)
00158                 dumpCoordinateForSVG<float>(d.col(indexes_bf[i]), 100, 1,  "stroke=\"black\" fill=\"green\"");
00159         dumpCoordinateForSVG<float>(q, 100, 1,  "stroke=\"black\" fill=\"blue\"");
00160         cout << "</svg>" << endl;
00161         */
00162         
00163         //cout << "Average KDTree visit count: " << double(totKDTreeVisitCount) * 100. / double(itCount * d.cols()) << " %" << endl;
00164         
00165         
00166         return 0;
00167 }


libnabo
Author(s): St├ęphane Magnenat
autogenerated on Thu Sep 10 2015 10:54:55