axialdf.cpp

Go to the documentation of this file.

/* axialdf.cpp: Axial distribution functions
 * Copyright (C) 2020  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/axialdf.h"
#include "atomprobe/helper/aptAssert.h"

#include <cstring>

#ifdef _OPENMP
    #include <omp.h>
#endif

#ifdef DEBUG
#include <limits>
#endif

#include "atomprobe/helper/constants.h"

namespace AtomProbe
{

using std::vector;

unsigned int generate1DAxialDistHist(const vector<Point3D> &pointList, K3DTreeExact &tree,
        const Point3D &axisDir, float distMax, vector<unsigned int> &histogram) 
{
    ASSERT(fabs(axisDir.sqrMag() -1.0f) < sqrt(std::numeric_limits<float>::epsilon()));
    ASSERT(!histogram.empty());

    std::fill(histogram.begin(),histogram.end(),0);

    if(pointList.empty())
        return 0;

    BoundCube cube;
    cube.setBounds(pointList);

    float maxSqrDist = distMax*distMax;

#ifdef _OPENMP
    vector<vector<unsigned int> > threadHistograms;
    threadHistograms.resize(omp_get_num_threads());
    for(unsigned int ui=0;ui<threadHistograms.size();ui++)
    {
        threadHistograms[ui].resize(histogram.size(),0);
    }
#endif

    

    
    //Main r-max searching routine
#pragma omp parallel for
    for(unsigned int ui=0; ui<pointList.size(); ui++)
    {
        Point3D sourcePoint;
        //Go through each point and grab all points within the maximum distance
        //that we need
    
        vector<size_t> inSphereIdx;
            
        sourcePoint=pointList[ui];

        //Grab the nearest point. Will tag on completion
        tree.ptsInSphere(sourcePoint, distMax, inSphereIdx);

        //Loop through the points within the search radius and
        // update the tag radius    
        for(size_t uj=0;uj<inSphereIdx.size();uj++)
        {   
            float sqrDist;
            const Point3D &nearPt=tree.getPtRef(inSphereIdx[uj]);
    
            //Calculate the sq of the distance to the point
            sqrDist = nearPt.sqrDist(sourcePoint);
            
            //if sqrDist is = maxSqrdist then this will cause
            //the histogram indexing to trash alternate memory
            //- this is bad - prevent this please.
            // also skip self matching (zero distance)  
            if(sqrDist == 0.0f || sqrDist >= maxSqrDist)
                continue;
            
            //Compute the projection of
            // the point onto the axis of the
            // primary analysis direction
            float distance;
            distance=(nearPt-sourcePoint).dotProd(axisDir);
        
            //update the histogram with the new position.
            // centre of the distribution function lies at the analysis point,
            // and can be either negative or positive. 
            // Shift the zero to the center of the histogram
            int offset=(int)(((0.5f*distance)/distMax+0.5f)*(float)histogram.size());
            if(offset < (int)histogram.size() && offset >=0)
            {
#ifdef _OPENMP
                threadHistograms[omp_get_thread_num()][offset]++;
#else
                histogram[offset]++;
#endif
            }

        }
        
    }
#ifdef _OPENMP
    for(unsigned int ui=0;ui<threadHistograms.size();ui++)
    {
        for(unsigned int uj=0;uj<histogram.size();uj++)
            histogram[uj]+=threadHistograms[ui][uj];
    }
#endif

    return 0;
}

unsigned int generate1DAxialDistHistSweep(const vector<Point3D> &pointList, K3DTreeExact &tree,
        float distMax, float dTheta, float dPhi, ProgressBar &p,
            vector<vector<vector<unsigned int> > > &histogram)
{
    ASSERT(!histogram.empty());
    ASSERT(dTheta > 0.0f);
    ASSERT(dPhi > 0.0f);

    if(pointList.empty())
        return 0;

    BoundCube cube;
    cube.setBounds(pointList);

    float maxSqrDist = distMax*distMax;

    const size_t NUM_THETA=histogram.size();
    const size_t NUM_PHI=histogram[0].size();
    const size_t DIST_BINS=histogram[0][0].size();


    unsigned int *hist = new unsigned int[NUM_THETA*NUM_PHI*DIST_BINS];
    memset(hist,0,NUM_THETA*NUM_PHI*DIST_BINS*sizeof(unsigned int));

    Point3D **axisDir= new Point3D*[NUM_THETA];
    for(size_t ui=0;ui<NUM_THETA;ui++)
        axisDir[ui] = new Point3D[NUM_PHI];

    //Construct a ball of unit vectors, facing outwards
    for(size_t ui=0; ui<NUM_THETA; ui++)
    {
        for(size_t uj=0; uj<NUM_PHI; uj++)
        {
            float tmpPhi,tmpTheta;
            tmpTheta= dTheta*ui;
            tmpPhi = dPhi*uj;

            //theta is the angle from +z, to prevent duplication of the vector
            // and its negative, we restrict it to 0->90 degrees from +z
            ASSERT(tmpTheta>=0.0f && tmpTheta<=M_PI/2);
            //This should like somewhere between 0 and 360 degrees
            ASSERT(tmpPhi>=0&& tmpPhi<=2.1f*M_PI);
            axisDir[ui][uj].setISOSpherical(1.0,tmpTheta,tmpPhi);

            ASSERT(fabs(acos(axisDir[ui][uj][2]) - tmpTheta) < 0.001);

        }
    }

    
    //Main r-max searching routine

//construct worker threads
#pragma omp parallel for schedule(dynamic)
    for(unsigned int uSrcPt=0; uSrcPt<pointList.size(); uSrcPt++)
    {
        Point3D sourcePoint;
        sourcePoint=pointList[uSrcPt];

        //Go through each point and grab up to the maximum distance
        //that we need
        vector<size_t> inSphereIdx;
        tree.ptsInSphere(sourcePoint,distMax, inSphereIdx);



        //Loop through the points within the search radius and
        // update the tag radius
        for(size_t uPt=0; uPt<inSphereIdx.size(); uPt++)
        {
            float sqrDist;
            ASSERT(inSphereIdx[uPt] < tree.size());
            const Point3D &nearPt=tree.getPtRef(inSphereIdx[uPt]);

            //Calculate the sq of the distance to the point
            sqrDist = nearPt.sqrDist(sourcePoint);

            if(sqrDist == 0.0f)
                continue;

            //if sqrDist is = maxSqrdist then this will cause
            //the histogram indexing to trash alternate memory
            //- this is bad - prevent this please.
            if(sqrDist < maxSqrDist)
            {
                Point3D deltaPt;
                deltaPt=nearPt-sourcePoint;
                //Compute the SDM fit for all angular space
                for(size_t ui=0; ui<NUM_THETA; ui++)
                {
                    for(size_t uj=0; uj<NUM_PHI; uj++)
                    {

                        //Compute the projection of
                        // the point onto the axis of the
                        // primary analysis direction
                        float distance;
                        distance=deltaPt.dotProd(axisDir[ui][uj]);
                        int offset;

                        offset=(int)(((0.5f*distance)/distMax+0.5f)*(float)DIST_BINS);
                        if(offset >=0 && offset < (int)DIST_BINS)
                        {
                            #pragma omp atomic
                            hist[ui*DIST_BINS*NUM_PHI + uj*DIST_BINS+ offset]++;
                        }
                    }
                }
            }

        }
#ifdef _OPENMP
        if(!omp_get_thread_num())
#endif
            p.update((float)uSrcPt/(float)pointList.size()*100.0f);
    }

    for(unsigned int ui=0;ui<histogram.size();ui++)
    {
        for(unsigned int uj=0;uj<histogram[ui].size();uj++)
        {
            for(unsigned int uk=0;uk<histogram[ui][uj].size();uk++)
            {
                histogram[ui][uj][uk]=hist[ui*DIST_BINS*NUM_PHI + uj*DIST_BINS+ uk];
            }
        }
    }


    delete[] hist;

    for(size_t ui=0;ui<NUM_THETA;ui++)
        delete[] axisDir[ui]; 

    delete[] axisDir;
    return 0;
}
}