convexHull.cpp

Go to the documentation of this file.

/* 
 * convexHull.cpp - Wrapper for quickHull convex hull routines
 * Copyright (C) 2015  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 "atomprobe/algorithm/convexHull.h"
#include "atomprobe/helper/maths/misc.h"

//QHull library
extern "C"
{
    #include <libqhull/qhull_a.h>
}

#include <numeric>
#include <limits>

#include "atomprobe/helper/aptAssert.h"

namespace AtomProbe
{

using std::vector;

//Is the qhull library already initialised (true)
// or do we need to call qhull's init routines?
bool qhullInited=false;
//TODO: Work out where the payoff for this is
//grab size when doing convex hull calculations
const unsigned int HULL_GRAB_SIZE=4096;




void freeConvexHull()
{
    qh_freeqhull(true);
    qhullInited=false;
}

unsigned int doHull(unsigned int bufferSize, double *buffer, 
            vector<Point3D> &resHull, Point3D &midPoint, bool freeHullOnExit)
{
    if(qhullInited)
    {
        qh_freeqhull(true);
        qhullInited=false;
    }

    const int dim=3;
    //Now compute the new hull
    //Generate the convex hull
    //(result is stored in qh's globals :(  )
    //note that the input is "joggled" to 
    //ensure simplicial facet generation

    //Qhull >=2012 has a "feature" where it won't accept null arguments for the output
    // there is no clear way to shut it up.
    FILE *outSquelch=0;
#if defined(__linux__) || defined(__APPLE__) || defined(__BSD__)
    outSquelch=fopen("/dev/null","w");
#elif defined(__win32__) || defined(__win64__)
    outSquelch=fopen("NUL","w");
#endif

    if(!outSquelch)
    {
        //Give up, just let qhull output random statistics to stderr
        outSquelch=stderr;
    }

    qh_new_qhull(   dim,
            bufferSize,
            buffer,
            false,
            (char *)"qhull QJ", //Joggle the output, such that only simplical facets are generated
            outSquelch, //QHULL's interface is bizarre, no way to set null pointer in qhull 2012 - result is inf. loop in qhull_fprintf and error reporting func. 
            outSquelch);
    qhullInited=true;

    if(outSquelch !=stderr)
    {
        fclose(outSquelch);
    }

    unsigned int numPoints=0;
    //count points
    //--    
    //OKay, whilst this may look like invalid syntax,
    //qh is actually a macro from qhull
    //that creates qh. or qh-> as needed
    numPoints = qh num_vertices;    
    //--    
    midPoint=Point3D(0,0,0);    

    if(!numPoints)
        return 0;

    //store points in vector
    //--
    try
    {
        resHull.resize(numPoints);  
    }
    catch(std::bad_alloc)
    {
        return HULL_ERR_NO_MEM;
    }
    //--

    //Compute mean point
    //--
    vertexT *vertex;
    vertex= qh vertex_list;
    int curPt=0;
    while(vertex != qh vertex_tail)
    {
        resHull[curPt]=Point3D(vertex->point[0],
                vertex->point[1],
                vertex->point[2]);
        midPoint+=resHull[curPt];
        curPt++;
        vertex = vertex->next;
    }
    midPoint*=1.0f/(float)numPoints;
    //--
    if(freeHullOnExit)
    {
        qh_freeqhull(true);
        qhullInited=false;
    
    }

    return 0;
}

unsigned int computeConvexHull(const vector<IonHit> &data, unsigned int *progress,
                bool (*callback)(bool),std::vector<Point3D> &curHull, bool freeHull)
{

    size_t numPts;
    numPts=data.size();
    //Easy case of no data
    if(numPts < 4)
        return 0;

    double *buffer;
    double *tmp;
    //Use malloc so we can re-alloc
    buffer =(double*) malloc(HULL_GRAB_SIZE*3*sizeof(double));

    if(!buffer)
        return HULL_ERR_NO_MEM;

    size_t bufferOffset=0;

    //Do the convex hull in steps for two reasons
    // 1) qhull chokes on large data
    // 2) we need to run the callback every now and again, so we have to
    //   work in batches.
    Point3D midPoint;
    float maxSqrDist=-1;
    
    const size_t PROGRESS_REDUCE = 5000;

    size_t progressReduce=PROGRESS_REDUCE;
    size_t n=0;
    for(size_t uj=0; uj<data.size(); uj++)
    {
        //Do contained-in-sphere check
        if(curHull.empty() || midPoint.sqrDist(data[uj].getPos())>= maxSqrDist)
        {
            //Copy point data into hull buffer
            buffer[3*bufferOffset]=data[uj].getPos()[0];
            buffer[3*bufferOffset+1]=data[uj].getPos()[1];
            buffer[3*bufferOffset+2]=data[uj].getPos()[2];
            bufferOffset++;

            //If we have hit the hull grab size, perform a hull

            if(bufferOffset == HULL_GRAB_SIZE)
            {
                bufferOffset+=curHull.size();
                tmp=(double*)realloc(buffer,
                             3*bufferOffset*sizeof(double));
                if(!tmp)
                {
                    free(buffer);

                    return HULL_ERR_NO_MEM;
                }

                buffer=tmp;
                //Copy in the old hull
                for(size_t uk=0; uk<curHull.size(); uk++)
                {
                    buffer[3*(HULL_GRAB_SIZE+uk)]=curHull[uk][0];
                    buffer[3*(HULL_GRAB_SIZE+uk)+1]=curHull[uk][1];
                    buffer[3*(HULL_GRAB_SIZE+uk)+2]=curHull[uk][2];
                }

                unsigned int errCode=0;
                
                errCode=doHull(bufferOffset,buffer,curHull,midPoint,freeHull);
                if(errCode)
                {
                    free(buffer);
                    return errCode;
                }


                //Now compute the min sqr distance
                //to the vertex, so we can fast-reject
                maxSqrDist=std::numeric_limits<float>::max();
                for(size_t ui=0; ui<curHull.size(); ui++)
                    maxSqrDist=std::min(maxSqrDist,curHull[ui].sqrDist(midPoint));
                //reset buffer size
                bufferOffset=0;
            }

        }
        n++;

        //Update the progress information, and run callback periodically
        if(!progressReduce--)
        {
            if(!(*callback)(false))
            {
                free(buffer);
                return HULL_ERR_USER_ABORT;
            }

            *progress= (unsigned int)((float)(n)/((float)numPts)*100.0f);

            progressReduce=PROGRESS_REDUCE;
        }
    }

    //Need at least 4 objects to construct a sufficiently large buffer
    if(bufferOffset + curHull.size() > 4)
    {
        //Re-allocate the buffer to determine the last hull size
        tmp=(double*)realloc(buffer,
                             3*(bufferOffset+curHull.size())*sizeof(double));
        if(!tmp)
        {
            free(buffer);
            return HULL_ERR_NO_MEM;
        }
        buffer=tmp;

        #pragma omp parallel for
        for(size_t ui=0; ui<curHull.size(); ui++)
        {
            buffer[3*(bufferOffset+ui)]=curHull[ui][0];
            buffer[3*(bufferOffset+ui)+1]=curHull[ui][1];
            buffer[3*(bufferOffset+ui)+2]=curHull[ui][2];
        }

        unsigned int errCode=doHull(bufferOffset+curHull.size(),buffer,curHull,midPoint,freeHull);

        if(errCode)
        {
            free(buffer);
            //Free the last convex hull mem
            return errCode;
        }
    }


    free(buffer);
    return 0;
}


//FIXME: This needs to use the new qhull 2015 we have implemented
unsigned int GetReducedHullPts(const vector<IonHit> &points, float reductionDim,  
        unsigned int *progress, bool (*callback)(bool), vector<IonHit> &pointResult)
{
    //TODO: This could be made to use a fixed amount of ram, by
    //partitioning the input points, and then
    //computing multiple hulls.
    //Take the resultant hull points, then hull again. This would be
    //much more space efficient, and more easily parallellised
    //Alternately, compute a for a randoms K set of points, and reject
    //points that lie in the hull from further computation

    //Need at least 4 points to define a hull in 3D
    if(points.size() < 4) 
        return 1;

    unsigned int dummyProg;
    vector<Point3D> theHull;
    if(computeConvexHull(points,progress,callback,theHull,false))
        return 2;

    Point3D midPoint(0,0,0);
    for(size_t ui=0;ui<theHull.size();ui++)
        midPoint+=theHull[ui];
    midPoint *= 1.0f/(float)theHull.size();

    //Now we will find the mass/volume centroid of the hull
    //by constructing sets of pyramids.
    //We cannot use the simple midpoint, as point distribution
    //around the hull surface may be uneven
    
    //Here the faces of the hull are the bases of 
    //pyramids, and the midpoint is the apex
    Point3D hullCentroid(0.0f,0.0f,0.0f);
    float massPyramids=0.0f;
    //Run through the faced list
    facetT *curFac = qh facet_list;
    
    while(curFac != qh facet_tail)
    {
        vertexT *vertex;
        Point3D pyramidCentroid;
        unsigned int ui;
        Point3D ptArray[3];
        float vol;

        //find the pyramid volume + centroid
        pyramidCentroid = midPoint;
        ui=0;
        
        //This assertion fails, some more processing is needed to be done to break
        //the facet into something simplical
        ASSERT(curFac->simplicial);
        vertex  = (vertexT *)curFac->vertices->e[ui].p;
        while(vertex)
        {   //copy the vertex info into the pt array
            (ptArray[ui])[0]  = vertex->point[0];
            (ptArray[ui])[1]  = vertex->point[1];
            (ptArray[ui])[2]  = vertex->point[2];
        
            //aggregate pyramidal points    
            pyramidCentroid += ptArray[ui];
        
            //increment before updating vertex
            //to allow checking for NULL termination    
            ui++;
            vertex  = (vertexT *)curFac->vertices->e[ui].p;
        
        }
        
        //note that this counter has been post incremented. 
        ASSERT(ui ==3);
        vol = pyramidVol(ptArray,midPoint);
        
        ASSERT(vol>=0);
        
        //Find the midpoint of the pyramid, this will be the 
        //same as its centre of mass.
        pyramidCentroid*= 0.25f;
        hullCentroid = hullCentroid + (pyramidCentroid*vol);
        massPyramids+=vol;

        curFac=curFac->next;
    }

    hullCentroid *= 1.0f/massPyramids;

    size_t badFacetCount=0;
    float minDist=std::numeric_limits<float>::max();
    //find the smallest distance between the centroid and the
    //convex hull
        curFac=qh facet_list;
    while(curFac != qh facet_tail)
    {
        float temp;
        Point3D vertexPt[3];
        
        //The shortest distance from the plane to the point
        //is the dot product of the UNIT normal with 
        //A-B, where B is on plane, A is point in question
        for(unsigned int ui=0; ui<3; ui++)
        {
            vertexT *vertex;
            //grab vertex
            vertex  = ((vertexT *)curFac->vertices->e[ui].p);
            vertexPt[ui] = Point3D(vertex->point[0],vertex->point[1],vertex->point[2]);
        }

        //Some facets can be degenerate. Skip if this occurs
        //--
        bool isBad;
        isBad=triIsDegenerate(vertexPt[0],vertexPt[1],vertexPt[2]);

        if(isBad)
        {
            badFacetCount++;
            curFac=curFac->next;
            continue;
        }
        //--

    

        //Find the distance between hullcentroid and a given facet
        temp = distanceToFacet(vertexPt[0],vertexPt[1],vertexPt[2],
                    Point3D(curFac->normal[0],curFac->normal[1],curFac->normal[2]),hullCentroid);

        if(temp < minDist)
            minDist = temp;
    
        curFac=curFac->next;
    }

    if(badFacetCount)
    {
        if(badFacetCount == qh num_facets)
        {
            return HULL_SURFREDUCE_NEGATIVE_SCALE_FACT;
        }
    }

    //shrink the convex hull such that it lies at
    //least reductionDim from the original surface of
    //the convex hull
    float scaleFactor;
    scaleFactor = 1  - reductionDim/ minDist;

    if(scaleFactor < 0.0f)
        return HULL_SURFREDUCE_NEGATIVE_SCALE_FACT;
    
    
    //now scan through the input points and see if they
    //lie in the reduced convex hull
    vertexT *vertex = qh vertex_list;   

    unsigned int ui=0;
    while(vertex !=qh vertex_tail)
    {
        //Translate around hullCentroid before scaling, 
        //then undo translation after scale
        //Modify the vertex data such that it is scaled around the hullCentroid
        vertex->point[0] = (vertex->point[0] - hullCentroid[0])*scaleFactor + hullCentroid[0];
        vertex->point[1] = (vertex->point[1] - hullCentroid[1])*scaleFactor + hullCentroid[1];
        vertex->point[2] = (vertex->point[2] - hullCentroid[2])*scaleFactor + hullCentroid[2];
    
        vertex = vertex->next;
        ui++;
    }

    //if the dot product of the normal with the point vector of the
    //considered point P, to any vertex on all of the facets of the 
    //convex hull F1, F2, ... , Fn is negative,
    //then P does NOT lie inside the convex hull.
    pointResult.reserve(points.size()/2);
    curFac = qh facet_list;
    
    //minimum distance from centroid to convex hull
    for(size_t uj=points.size(); uj--;)
    {
        float fX,fY,fZ;
        double *ptArr,*normalArr;
        fX =points[uj][0];
        fY = points[uj][1];
        fZ = points[uj][2];
        
        //loop through the facets
        curFac = qh facet_list;
        while(curFac != qh facet_tail)
        {
            //Dont ask. It just grabs the first coords of the vertex
            //associated with this facet
            ptArr = ((vertexT *)curFac->vertices->e[0].p)->point;
            
            normalArr = curFac->normal;
            //if the dotproduct is negative, then the point vector from the 
            //point in question to the surface is in opposite to the outwards facing
            //normal, which means the point lies outside the hull   
            if (dotProduct( (float)ptArr[0]  - fX,
                    (float)ptArr[1] - fY,
                    (float)ptArr[2] - fZ,
                    normalArr[0], normalArr[1],
                    normalArr[2]) >= 0)
            {
                curFac=curFac->next;
                continue;
            }
            goto reduced_loop_next;
        }
        //we passed all tests, point is inside convex hull
        pointResult.push_back(points[uj]);
        
reduced_loop_next:
    ;
    }

    freeConvexHull();

    return 0;
}

#ifdef DEBUG

bool dummyCallback(bool b)
{
    return true;
}
#include <iostream>
#include <sys/time.h>
using std::cerr;
using std::endl;
bool testConvexHull()
{
        std::vector<IonHit> ions;

        unsigned int NUM_IONS=100000;
        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);
        }


    vector<Point3D> hullVertices;

    unsigned int dummy=0;
    //Now compute the convex hull
    if(computeConvexHull(ions,&dummy,dummyCallback,hullVertices,true))
        return false;

        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);
        }
    
    vector<IonHit> reducedHullPts;
    if(GetReducedHullPts(ions,SCALE*0.1,&dummy,dummyCallback, reducedHullPts))
        return false;

    if(reducedHullPts.empty())
        return false;

    
    return true;
}
#endif

}