proxigram.cpp

Go to the documentation of this file.

#include "atomprobe/atomprobe.h"


using namespace std;
using namespace AtomProbe;

int main( int argc, char *argv[])
{
    Voxels<float> v;

    const unsigned int VOXEL_COUNT=30;
    v.resize(VOXEL_COUNT,VOXEL_COUNT,VOXEL_COUNT);

    const float VOX_SIZE=2;
    v.setBounds(Point3D(-VOX_SIZE,-VOX_SIZE,-VOX_SIZE),
        Point3D(VOX_SIZE,VOX_SIZE,VOX_SIZE) );
    v.fill(0);

    //Consider two spheres at pCentreA, and pCentreB
    // of given radius SPHERE_RAD
    const float SPHERE_RAD=1.0f;
    const float SPHERE_RAD_SQR=SPHERE_RAD*SPHERE_RAD;
    const Point3D pCentreA(-0.5f,-0.5f,-0.5f); 
    const Point3D pCentreB(0.5f,0.5f,0.5f); 


    cerr << "Creating synthetic voxel field...";    
    for(auto ui=0;ui<VOXEL_COUNT;ui++)
    {
        for(auto uj=0;uj<VOXEL_COUNT;uj++)
        {
            for(auto uk=0;uk<VOXEL_COUNT;uk++)
            {
                Point3D p;
                p=v.getPoint(ui,uj,uk);

                if(p.sqrDist(pCentreA) < SPHERE_RAD_SQR
                    || p.sqrDist(pCentreB) < SPHERE_RAD_SQR)
                {
                    v.setData(ui,uj,uk,1);
                }
            }
        }
            
    }
    cerr << "Done" << endl;

    cerr << "Extracting isosurface (triangle soup)..." ;
    vector<TriangleWithVertexNorm> surfTris;
    marchingCubes(v,0.5,surfTris);
    cerr << "Done (" << surfTris.size() << " triangles)" << endl;

    //FIXME: This is not needed strictly speaking
    // but it has got advantages if we used a different algorithm
    //Convert to mesh object
    cerr << "Constructing mesh (deduplicating soup)...";
    Mesh m;
    m.nodes.resize(surfTris.size()*3);
    m.triangles.resize(surfTris.size());

    for(auto ui=0;ui<surfTris.size();ui++) 
    {
        for(auto uj=0;uj<3;uj++)
        {
            m.nodes[3*ui+ uj]=surfTris[ui].p[uj];
            m.triangles[ui].p[uj]=3*ui+uj;
        }
    }
    //Remove surf triangles, as we will now work with mesh
    surfTris.clear(); 

    //Remove duplciated vertices, and fuse into a single mesh object
    m.mergeDuplicateVertices(
        2.0f*sqrt(std::numeric_limits<float>::epsilon()));
    cerr  << "Done" << endl;

    //TODO: Optimise : split mesh into disconnected parts

    cerr << "Generating synthetic points...";   
    const unsigned int NUM_POINTS=1000;
    vector<Point3D> randomPoints;
    randomPoints.resize(NUM_POINTS);
    BoundCube bc;
    v.getBounds(bc);    
    gsl_rng *r = randGen.getRng();
    for(auto ui=0;ui<NUM_POINTS;ui++)
    {
        float x,y,z;
        x=2.0*VOX_SIZE*gsl_rng_uniform(r) -VOX_SIZE;
        y=2.0*VOX_SIZE*gsl_rng_uniform(r) -VOX_SIZE;
        z=2.0*VOX_SIZE*gsl_rng_uniform(r) -VOX_SIZE;
        randomPoints[ui]=Point3D(x,y,z);

    }
    cerr << "Done" << endl;
    
    cerr << "Running signed-distance computation...";
    //Run points-inside algorithm
    vector<bool> inside;

    //TODO: Thread me, as w/o threading, progress is not useful
    unsigned int progress;
    m.pointsInside(randomPoints,inside,progress);

    vector<float> distances;
    distances.resize(randomPoints.size());

#pragma omp parallel for
    for(auto ui=0;ui<randomPoints.size();ui++)
    {

        m.getNearestTri(randomPoints[ui],
            distances[ui]); 
        if(inside[ui])
            distances[ui]=-fabs(distances[ui]);
    }
    
    const char *OUTPUT_FILE="proxigram.txt";
    ofstream fOut(OUTPUT_FILE);

    if(!fOut)
    {
        cerr << "failed to open output file" << endl;
        return 1;
    }

    vector<float> hist;
    const unsigned int NBINS=50;
    float dMin,dMax;
    dMin=*std::min_element(distances.begin(),distances.end());
    dMax=*std::max_element(distances.begin(),distances.end());
    float delta = (dMax-dMin)/NBINS;
    linearHistogram(distances,dMin,dMax,delta,hist);

    for(auto ui=0;ui<hist.size();ui++)
        fOut << ui*delta+dMin << "\t" << hist[ui] << endl;
    
    cerr << "Done. Output written to :" << OUTPUT_FILE << endl;
}