point3D.cpp

Go to the documentation of this file.

/* point3D.cpp : 3D point primitive for libatomprobe
 * 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 "atomprobe/primitives/point3D.h"
#include "atomprobe/helper/endianTest.h"

#include "atomprobe/helper/aptAssert.h"
#include "atomprobe/helper/constants.h"

#include <limits>
#include <cstdlib>
#include <iostream>
#include <algorithm>

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


//FIXME: This is not the right place for this
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_blas.h>

namespace AtomProbe{

bool Point3D::parse(const std::string &str)
{
    using std::string;
    using std::vector;

    //Needs to be at minimum #,#,#
    if(str.size()< 5)
        return false;

    string tmpStr;
    tmpStr=stripWhite(str);

    //are we using polar notation? This is indicated by < ... > 
    bool polarNotation = (tmpStr[0] == '<' && tmpStr[tmpStr.size()-1] == '>');

    //Two strings must be in sync
    std::string allowableStartChars, allowableEndChars;
    allowableStartChars="([{<'";
    allowableEndChars=")]}>'";

    //Find the start/end chars
    size_t startPos,endPos;
    startPos=allowableStartChars.find(tmpStr[0]);
    endPos=allowableEndChars.find(tmpStr[tmpStr.size()-1]);


    //Strip the start/end chars 
    if(startPos !=std::string::npos && endPos != std::string::npos)
        tmpStr=tmpStr.substr(1,tmpStr.size()-1);
    else if (startPos !=endPos)
        return false; //we had one start bracket, but not the other...

    //First try splitting with non-whitespace separators,
    // there should be exactly 3 components
    vector<string> components;
    const char *NONWHITE_SEPARATOR=",;|_";
    splitStrsRef(tmpStr.c_str(),NONWHITE_SEPARATOR,
            components);
    if(components.size()!=3)
        return false;
    components.clear();

    //Now try splitting with whitespace components, dropping empty
    // strings. As we have already checked the non-whitespace components
    // additional components must be whitespace only, which is fine.

    const char *ALLOWABLE_SEPARATORS=",; \t|_";

    splitStrsRef(tmpStr.c_str(),ALLOWABLE_SEPARATORS,
            components);
    for(size_t ui=0;ui<components.size();ui++)
        components[ui]=stripWhite(components[ui]);
    
    //Drop the blank bits from the field
    vector<string>::iterator rmIt;
    rmIt=std::remove(components.begin(),components.end(),string(""));
    components.erase(rmIt,components.end());

    if(components.size()!=3)
        return false;

    float p[3];
    for(size_t ui=0;ui<3;ui++)
    {
        if(stream_cast(p[ui],components[ui]))
            return false;
    }

    if(!polarNotation)
        setValueArr(p);
    else
    {
        //r,theta,phi, as users should input
        // <r,theta,phi>
        //Convert from degrees to radiians
        p[1]*=M_PI/180; 
        p[2]*=M_PI/180; 

        //set spherical co-ordinates
        setISOSpherical(p[0],p[1],p[2]);
    }


    
    return true;
}


float Point3D::operator[](unsigned int ui) const 
{
    ASSERT(ui < 3);
    return value[ui];
}

float &Point3D::operator[](unsigned int ui) 
{
    ASSERT(ui < 3);
    return value[ui];
}

void Point3D::copyValueArr(float *valArr) const
{
    ASSERT(valArr);
    //compiler should unroll this automatically
    for(unsigned int ui=0; ui<3; ui++)
    {
        *(valArr+ui) = *(value+ui);
    }
}

const Point3D &Point3D::operator-=(const Point3D &pt)
{
    for(unsigned int ui=0;ui<3; ui++)
        value[ui]-= pt.value[ui];
    
    return *this;
}

bool Point3D::operator==(const Point3D &pt) const
{
    return (value[0] == pt.value[0] && value[1] == pt.value[1] && value[2] == pt.value[2]);
}

const Point3D &Point3D::operator=(const Point3D &pt)
{
    value [0] = pt.value[0];
    value [1] = pt.value[1];
    value [2] = pt.value[2];
    return *this;
}

const Point3D &Point3D::operator+=(const Point3D &pt)
{
    for(unsigned int ui=0;ui<3; ui++)
        value[ui]+= pt.value[ui];
    
    return *this;
}

const Point3D Point3D::operator+(const Point3D &pt) const
{
    Point3D ptTmp;
    
    for(unsigned int ui=0;ui<3; ui++)
        ptTmp.value[ui] = value[ui]  + pt.value[ui];
    
    return ptTmp;
}

const Point3D Point3D::operator+(float f)  const
{
    Point3D tmp;
    for(unsigned int ui=0;ui<3; ui++)
        tmp.value[ui] = value[ui]  + f;
    
    return tmp;
}

const Point3D Point3D::operator-(const Point3D &pt) const
{
    Point3D ptTmp;
    
    for(unsigned int ui=0;ui<3; ui++)
        ptTmp.value[ui] = value[ui]  - pt.value[ui];
    
    return ptTmp;
}

const Point3D Point3D::operator-() const
{
    Point3D ptTmp;
    
    for(unsigned int ui=0;ui<3; ui++)
        ptTmp.value[ui] = -value[ui];
    
    return ptTmp;
}

Point3D Point3D::operator*=(const float scale)
{
    value[0] = value[0]*scale;
    value[1] = value[1]*scale;
    value[2] = value[2]*scale;
    
    return *this;
}

const Point3D Point3D::operator*(float scale) const
{
    Point3D tmpPt;

    tmpPt.value[0] = value[0]*scale;
    tmpPt.value[1] = value[1]*scale;
    tmpPt.value[2] = value[2]*scale;
    
    return tmpPt;
}

const Point3D Point3D::operator*(const Point3D &pt) const
{
    Point3D tmpPt;
    
    tmpPt.value[0] = value[0]*pt[0];
    tmpPt.value[1] = value[1]*pt[1];
    tmpPt.value[2] = value[2]*pt[2];
    
    return tmpPt;
}

const Point3D Point3D::operator/(float scale) const
{
    Point3D tmpPt;

    scale = 1.0f/scale;
    tmpPt.value[0] = value[0]*scale;
    tmpPt.value[1] = value[1]*scale;
    tmpPt.value[2] = value[2]*scale;

    return tmpPt;
}

const Point3D Point3D::operator/(const Point3D &pt) const
{
    Point3D tmpPt;
    
    tmpPt.value[0] = value[0]/pt[0];
    tmpPt.value[1] = value[1]/pt[1];
    tmpPt.value[2] = value[2]/pt[2];
    
    return tmpPt;
}


float Point3D::sqrDist(const Point3D &pt) const
{
    return (pt.value[0]-value[0])*(pt.value[0]-value[0])+
        (pt.value[1]-value[1])*(pt.value[1]-value[1])+
        (pt.value[2]-value[2])*(pt.value[2]-value[2]);
}
        
float Point3D::dotProd(const Point3D &pt) const
{
    //Return the inner product
    return value[0]*pt.value[0] + value[1]*pt.value[1] + value[2]*pt.value[2];
}

Point3D Point3D::crossProd(const Point3D &pt) const
{
    Point3D cross;

    cross.value[0] = (pt.value[2]*value[1] - pt.value[1]*value[2]);
    cross.value[1] = (value[2]*pt.value[0] - pt.value[2]*value[0]);
    cross.value[2] = (value[0]*pt.value[1] - value[1]*pt.value[0]);



    return cross;
}

void Point3D::extend(float distance)
{
    ASSERT(sqrMag() > 0.0f);
    
    Point3D p;
    p=*this;
    p.normalise();
    *this+=p*distance;
}

bool Point3D::insideBox(const Point3D &farPoint) const
{
    
    return (value[0] < farPoint.value[0] && value[0] >=0) &&
        (value[1] < farPoint.value[1] && value[1] >=0) &&
        (value[2] < farPoint.value[2] && value[2] >=0);
}

bool Point3D::insideBox(const Point3D &lowPt,const Point3D &highPt) const
{
    
    return (value[0] < highPt.value[0] && value[0] >=lowPt.value[0]) &&
        (value[1] < highPt.value[1] && value[1] >=lowPt.value[1]) &&
        (value[2] < highPt.value[2] && value[2] >=lowPt.value[2]);
}

//This is different to +=, because it generates no return value
void Point3D::add(const Point3D &obj)
{
    value[0] = obj.value[0] + value[0];
    value[1] = obj.value[1] + value[1];
    value[2] = obj.value[2] + value[2];
}

float Point3D::sqrMag() const
{
    return value[0]*value[0] + value[1]*value[1] + value[2]*value[2];
}

void Point3D::normalise()
{
    float mag = sqrtf(sqrMag());

    value[0]/=mag;
    value[1]/=mag;
    value[2]/=mag;

}

Point3D Point3D::normal() const
{
    Point3D ret;
    float mag = sqrtf(sqrMag());

    ret[0]=value[0]/mag;
    ret[1]=value[1]/mag;
    ret[2]=value[2]/mag;

    return ret;
}
void Point3D::negate() 
{
    value[0] = -value[0];
    value[1] = -value[1];
    value[2] = -value[2];
}

void Point3D::reciprocal()
{
    value[0] = 1.0/value[0];
    value[1] = 1.0/value[1];
    value[2] = 1.0/value[2];
}

bool Point3D::orthogonalise(const Point3D &pt)
{
    Point3D crossp;
    crossp=this->crossProd(pt);

    //They are co-linear, or near-enough to be not resolvable.
    if(crossp.sqrMag()  < sqrt(std::numeric_limits<float>::epsilon()))
        return false;
    crossp.normalise();

    crossp=crossp.crossProd(pt);
    *this=crossp.normal()*sqrt(this->sqrMag()); 

    return true;
}

float Point3D::angle(const Point3D &pt) const
{
    return acos(dotProd(pt)/(sqrtf(sqrMag()*pt.sqrMag())));
}

void Point3D::getCentroid(const std::vector<Point3D> &points,Point3D &centroid)
{
    centroid=Point3D(0,0,0);
    size_t nPoints=points.size();
#ifdef _OPENMP
    
    //Parallel version
    //--
    std::vector<Point3D> centroids(omp_get_max_threads(),Point3D(0,0,0));
#pragma omp parallel for 
    for(size_t ui=0;ui<nPoints;ui++)
        centroids[omp_get_thread_num()]+=points[ui];

    for(size_t ui=0;ui<centroids.size();ui++)
        centroid+=centroids[ui];
    //--

#else
    for(unsigned int ui=0;ui<nPoints;ui++)
        centroid+=points[ui];
#endif
    
    centroid*=1.0f/(float)nPoints;
}
std::ostream& operator<<(std::ostream &stream, const Point3D &pt)
{
       stream << "(" << pt.value[0] << "," << pt.value[1] 
                   << "," << pt.value[2] << ")";
       return stream;
}


void Point3D::setISOSpherical(float r,float theta, float phi)
{
    value[0] = r*sin(theta)*cos(phi);
    value[1] = r*sin(theta)*sin(phi);
    value[2] = r*cos(theta);
}


void Point3D::getISOSpherical(float &r,float &theta, float &phi) const
{
    r=sqrt(sqrMag());
    if(r < std::numeric_limits<float>::epsilon())
    {
        theta=phi=0;
        return;
    }
    theta =acos(value[2]/r);
    phi =atan2(value[1],value[0]);
}

void Point3D::transform3x3(const gsl_matrix* matrix)
{
    float newValue[3];
    for(unsigned int row=0;row<3;row++)
    {
        newValue[row]=0;
    
        for(unsigned int col=0;col<3;col++)
        {

            newValue[row]+= gsl_matrix_get(matrix,row,col)*value[col];
        }
    }

    value[0]=newValue[0];
    value[1]=newValue[1];
    value[2]=newValue[2];
}

#ifdef __LITTLE_ENDIAN__
void Point3D::switchEndian()
{
    floatSwapBytes(&value[0]);
    floatSwapBytes(&value[1]);
    floatSwapBytes(&value[2]);
}
#endif
#ifdef DEBUG

#include "helper/helpFuncs.h"

bool testPoint3D()
{

    Point3D p1;
    p1.parse("(1,0,0)");
    Point3D p2;
    p2.parse("<1,0,0>");

    TEST(p1 == Point3D(1,0,0),"Parse check");
    TEST(p2 == Point3D(0,0,1),"parse check (spherical)");

    TEST(p2.mag() == 1.0,"mag check");

    p1=Point3D(1,1,1);

    TEST(TOL_EQ(p1.dotProd(p1),(3.0)),"Dot product");

    p1=Point3D(1,0,0);
    p2=Point3D(0,1,0);

    TEST(TOL_EQ(
        (p1.crossProd(p2) - Point3D(0,0,1)).sqrMag(),0),"cross prod"); 
        
    return true;
}

#endif
}