kd-example.cpp

Go to the documentation of this file.

/* kd-example.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;

unsigned int progress=0;

#define TIME_START() timeval TIME_DEBUG_t; gettimeofday(&TIME_DEBUG_t,NULL);
#define TIME_END() float TIME_DELTA; {timeval TIME_DEBUG_tend; gettimeofday(&TIME_DEBUG_tend,NULL); \
TIME_DELTA=(TIME_DEBUG_tend.tv_sec - TIME_DEBUG_t.tv_sec) + ((float)TIME_DEBUG_tend.tv_usec-(float)TIME_DEBUG_t.tv_usec)/1.0e6;}


#ifndef M_PI
#define M_PI 3.141596254
#endif

bool callback()
{
    static unsigned int lastProgress=0;
    if(lastProgress < progress)
    {
        cerr << ".";
        lastProgress=progress+1;
    }

    return true;
}

int main()
{

    vector<IonHit> ions;

    unsigned int NUM_IONS=50000;
    const float SCALE=100;
    //initialise random number generator using 
    // current system time
    srand(time(NULL));

    
    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;

    //set progress callback
    kdExact.setCallback(callback);
    kdExact.setProgressPointer(&progress);
    //set data - note that the ions will be taken by the tree by default.
    // - this can be overridden. See docs for K3DTreeExact::resetPts
    kdExact.resetPts(ions);

    cerr << "Building tree with " << NUM_IONS/(float)1000000 << " M ions..."  << endl << " |";

    {
    TIME_START();
    //build tree
    kdExact.build();

    TIME_END();
    cerr << "| done, in " << TIME_DELTA << " seconds" << endl;
    }

    BoundCube bc;
    kdExact.getBoundCube(bc);

    //1NN search
    {
    TIME_START();
    
    cerr << "Finding all 1NNs...";

    //In parallel (if enabled), find all first NNs. 
    // if setTag = false, this does not modify the tree, so is safe to call in parallel
    #pragma omp parallel for
    for(size_t ui=0;ui<kdExact.size();ui++)
    {
        //find the nearest pt, but do not tag it
        kdExact.findNearestUntagged(kdExact.getPtRef(ui),bc,false);
    }
    TIME_END();
    cerr << "done, in " << TIME_DELTA << " seconds" << endl;


    }   

    //radius search, slow
    const float SEARCH_RAD=4;
    const float SIDE_RATIO= SEARCH_RAD/SCALE;
    const float REL_VOLUME =4.0/3.0*M_PI*SIDE_RATIO*SIDE_RATIO*SIDE_RATIO*100.0;
    cerr << "In a cube containing " << kdExact.size() << "points, of side length " << SCALE <<  endl;
    cerr << " finding pts within " << SEARCH_RAD << "distance. (" << REL_VOLUME << "\% of points)" <<endl;
    {
    

    cerr << "Slow query" << endl;
    double nPts=0;
    TIME_START();
    for(size_t ui=0;ui<kdExact.size();ui++)
    {
        vector<size_t> result;  
        kdExact.findUntaggedInRadius(kdExact.getPtRef(ui),bc,SEARCH_RAD,result);
        nPts+=result.size();

        result.clear();
    }
    TIME_END();


    cerr << " Average :" << nPts/(double)kdExact.size() << " pts around each other pt. Took "; 
    cerr << TIME_DELTA << " seconds" << endl;
    }

    //radius search, fast
    {
    cerr << "Faster (bulk) query" << endl;
    double nPts=0;
    TIME_START();
    for(size_t ui=0;ui<kdExact.size();ui++)
    {
        vector<size_t> result;  
        kdExact.ptsInSphere(kdExact.getPtRef(ui),SEARCH_RAD,result);
        nPts+=result.size();

        result.clear();
    }
    TIME_END();


    cerr << " Average :" << nPts/(double)kdExact.size() << " pts around each other pt. Took "; 
    cerr << TIME_DELTA << " seconds" << endl;
    }
}