KDTest.cpp

Go to the documentation of this file.

/* KDTest.cpp: Sample program for spatial searching (eg Nearest Neighbour)
 * Copyright (C) 2014  Daniel Haley
 * 
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 * 
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 * 
 * You should have received a copy of the GNU General Public License
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
 */
#include <iostream>
#include <cstdlib>
#include <sys/time.h>

#include "atomprobe/atomprobe.h"

using namespace std;
using namespace AtomProbe;
    
const unsigned int NUM_IONS=1000;
const float SCALE=100;
const float SEARCH_RAD = 5;

unsigned int progress=0;

bool callback()
{
    return true;
}

void kdExactFuzz()
{
    //Generate a random cube
    vector<IonHit> ions;

    //initialise random number generator using 
    // current system time
    
    ions.resize(NUM_IONS);
    for(unsigned int ui=0;ui<NUM_IONS;ui++)
    {

        float xyzm[4]; 
        xyzm[0] = SCALE*((float)rand()/RAND_MAX-0.5);
        xyzm[1] = SCALE*((float)rand()/RAND_MAX-0.5);
        xyzm[2] = SCALE*((float)rand()/RAND_MAX-0.5);
        xyzm[3] = SCALE*((float)rand()/RAND_MAX-0.5);
        ions[ui].setHit(xyzm);
    }

    K3DTreeExact kdExact;
    kdExact.setCallback(callback);
    kdExact.setProgressPointer(&progress);
    kdExact.resetPts(ions);

    cerr << ".";
    kdExact.build();
    
    BoundCube bc;
    kdExact.getBoundCube(bc);

    for(unsigned int ui=1;ui<kdExact.size();ui++)
    {
        {   
        //Find the nearest using the repeated untagged method
        vector<size_t> resultA, resultB;    
        kdExact.findUntaggedInRadius(kdExact.getPtRef(ui),bc,SEARCH_RAD,resultA);
        kdExact.ptsInSphere(kdExact.getPtRef(ui),SEARCH_RAD,resultB);

        std::sort(resultA.begin(),resultA.end());
        std::sort(resultB.begin(),resultB.end());


        if(resultA.size() != resultB.size() || 
        !(std::equal(resultA.begin(),resultA.end(),resultB.begin())))
        {
            cerr << "FAIL! Mismatched" << resultA.size () << " vs " << resultB.size() << "pts, sourcept : " << ui <<  endl; 
        
            cerr << "A:" << endl;
            for(unsigned int uj=0;uj<resultA.size();uj++)
                cerr << "\t" << kdExact.getPtRef(resultA[uj]) << endl;
            cerr << " :" << endl;
            for(unsigned int uj=0;uj<resultB.size();uj++)
                cerr << "\t" << kdExact.getPtRef(resultB[uj]) << endl;

            exit(1);
        }       
        }
    }
}

bool testInexactKDTree()
{
    //Tiny example where we should get the root node as the answer
    vector<Point3D> pts;

    pts.push_back(Point3D(0,0,0));
    pts.push_back(Point3D(0,0.5,0.5));
    pts.push_back(Point3D(100,0,0));

    BoundCube bc(pts);

    K3DTreeApprox k3dTree;
    k3dTree.build(pts);

    //0,0.5,0.5  (aka pts[1]) should be closest to 0,0,0
    const Point3D *p;
    p = k3dTree.findNearest(pts[0],bc,0.00001);

    if( p->sqrDist(pts[1]) > 0.01)
    {
        cerr << "K3DTreeApprox findNearest check failed" << endl;
        return false;
    }

    return true;
}

int main()
{

    timeval timeNow;
    gettimeofday(&timeNow,NULL);

    if(!testInexactKDTree())
    {
        cerr << "TEST FAILED" << endl;
        return false;
    }


    cerr << "Comparing two NN search methods to ensure they give the same answers, over many queries" << endl;
    
    unsigned int randSeed=timeNow.tv_usec; 
    srand(randSeed);
    cerr << "Random seed :" << randSeed << endl;
    cerr << "SCALE :" << SCALE << " Num Ions per Build : " << NUM_IONS << endl;
    cerr << "Search radius :" << SEARCH_RAD << endl;

    const unsigned int ITERATIONS=100;
    cerr << "Attempting  " << ITERATIONS << " random trees, finding ions in < search radius, using every ion as centre." << endl;
    cerr << "Will repeat for  :" << ITERATIONS << " passes" << endl;
    for(unsigned int ui=1;ui<ITERATIONS; ui++)
    {
        kdExactFuzz();

        if( ui && !(ui%50))
            cerr << " Pass" << ui<< endl;
    }
    
    timeval timeEnd;
    gettimeofday(&timeEnd,NULL);
    float timeDelta;
    timeDelta = (timeEnd.tv_sec - timeNow.tv_sec)*1e3 + (timeEnd.tv_usec - timeNow.tv_usec)/1000;
    cerr << "Complete (" << timeDelta <<" msec), all searches gave the same results" << endl;   
}