boundcube.cpp

Go to the documentation of this file.

/* boundcube.cpp : bounding cube class
 * 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/boundcube.h"
#include "atomprobe/primitives/ionHit.h"

#include "atomprobe/helper/aptAssert.h"

#include <limits>

namespace AtomProbe{

BoundCube::BoundCube(const Point3D &p1, const Point3D &p2)
{
    for(unsigned int ui=0;ui<3;ui++)
    {
#ifdef DEBUG
        valid[ui][0] = valid[ui][1] = true;
#endif
        bounds[ui][0] =std::min(p1[ui],p2[ui]);
        bounds[ui][1] =std::max(p1[ui],p2[ui]);
    }
}

bool BoundCube::operator==(const BoundCube &other) const
{
    for(unsigned int ui=0; ui<3; ui++)
    {
        if(bounds[ui][0]!=other.bounds[ui][0] ||
        bounds[ui][1]!=other.bounds[ui][1])
            return false;
    }
    return true;
}
    
void BoundCube::setBounds(float xMin,float yMin,float zMin,
                   float xMax,float yMax,float zMax) 
{
        bounds[0][0]=xMin; bounds[1][0]=yMin; bounds[2][0]=zMin;
        bounds[0][1]=xMax; bounds[1][1]=yMax; bounds[2][1]=zMax;
#ifdef DEBUG
        valid[0][0]=true; valid[1][0]=true; valid[2][0]=true;
        valid[0][1]=true; valid[1][1]=true; valid[2][1]=true;
#endif
}

void BoundCube::setBounds(const BoundCube &b)
{
    for(unsigned int ui=0;ui<3;ui++)
    {
        bounds[ui][0] = b.bounds[ui][0];
        bounds[ui][1] = b.bounds[ui][1];
#ifdef DEBUG
        valid[ui][0] = b.valid[ui][0];
        valid[ui][1] = b.valid[ui][1];
#endif
    }
}

    


void BoundCube::setBound(unsigned int bound, unsigned int minMax, float value)
{
    ASSERT(bound <3 && minMax < 2);
    bounds[bound][minMax]=value;
#ifdef DEBUG
    valid[bound][minMax]=true;
#endif
}


float BoundCube::getBound(unsigned int bound, unsigned int minMax) const
{
    ASSERT(bound <3 && minMax < 2);
    ASSERT(valid[bound][minMax]==true);
    return bounds[bound][minMax];
}


void BoundCube::setInverseLimits()
{
    bounds[0][0] = std::numeric_limits<float>::max();
    bounds[1][0] = std::numeric_limits<float>::max();
    bounds[2][0] = std::numeric_limits<float>::max();
    
    bounds[0][1] = -std::numeric_limits<float>::max();
    bounds[1][1] = -std::numeric_limits<float>::max();
    bounds[2][1] = -std::numeric_limits<float>::max();

#ifdef DEBUG
    valid[0][0] =false;
    valid[1][0] =false;
    valid[2][0] =false;
    
    valid[0][1] =false;
    valid[1][1] =false;
    valid[2][1] =false;
#endif
}

void BoundCube::setLimits()
{
    bounds[0][0] = -std::numeric_limits<float>::max();
    bounds[1][0] = -std::numeric_limits<float>::max();
    bounds[2][0] = -std::numeric_limits<float>::max();
    
    bounds[0][1] = std::numeric_limits<float>::max();
    bounds[1][1] = std::numeric_limits<float>::max();
    bounds[2][1] = std::numeric_limits<float>::max();
    

#ifdef DEBUG
    valid[0][0] =false;
    valid[1][0] =false;
    valid[2][0] =false;
    
    valid[0][1] =false;
    valid[1][1] =false;
    valid[2][1] =false;
#endif
}
bool BoundCube::isFlat() const
{
    //Test the limits being inverted or equated
    for(unsigned int ui=0;ui<3; ui++)
    {
        if(fabs(bounds[ui][0] - bounds[ui][1]) < std::numeric_limits<float>::epsilon())
            return true;
    }   

    return false;
}

void BoundCube::expand(const BoundCube &b) 
{
    //Check both lower and upper limit.
    //Moving to other cubes value as needed

#ifdef DEBUG
    if(!b.isValid())
        return;
#endif
    //If self not valid, ensure that it will be after this run
    //if(!isValid())
    //  setInverseLimits();

    for(unsigned int ui=0; ui<3; ui++)
    {
        if(b.bounds[ui][0] < bounds[ui][0])
        {
               bounds[ui][0] = b.bounds[ui][0]; 
#ifdef DEBUG
               valid[ui][0] = true;
#endif
        }

        if(b.bounds[ui][1] > bounds[ui][1])
        {
               bounds[ui][1] = b.bounds[ui][1]; 
               #ifdef DEBUG
               valid[ui][1] = true;
        #endif
        }
    }
}

void BoundCube::expand(const Point3D &p) 
{
    //If self not valid, ensure that it will be after this run
    //ASSERT(isValid())
    for(unsigned int ui=0; ui<3; ui++)
    {
        //Check lower bound is lower to new pt
        if(bounds[ui][0] > p[ui])
               bounds[ui][0] = p[ui];

        //Check upper bound is upper to new pt
        if(bounds[ui][1] < p[ui])
               bounds[ui][1] = p[ui];
    }
}

void BoundCube::expand(float f) 
{
    //If self not valid, ensure that it will be after this run
    for(unsigned int ui=0; ui<3; ui++)
    {
        //Check lower bound is lower to new pt
        bounds[ui][0]-=f;

        //Check upper bound is upper to new pt
        bounds[ui][1]+=f;
    }
}

void BoundCube::setBounds(const Point3D *p, unsigned int n)
{
    bounds[0][0] = std::numeric_limits<float>::max();
    bounds[1][0] = std::numeric_limits<float>::max();
    bounds[2][0] = std::numeric_limits<float>::max();
    
    bounds[0][1] = -std::numeric_limits<float>::max();
    bounds[1][1] = -std::numeric_limits<float>::max();
    bounds[2][1] = -std::numeric_limits<float>::max();
    
    for(unsigned int ui=0;ui<n; ui++)
    {
        for(unsigned int uj=0;uj<3;uj++)
        {
            if(bounds[uj][0] > p[ui][uj])
            {
                   bounds[uj][0] = p[ui][uj];
#ifdef DEBUG
                    valid[uj][0] = true;
#endif
            }   


            if(bounds[uj][1] < p[ui][uj])
            {
                   bounds[uj][1] = p[ui][uj];       
#ifdef DEBUG
                    valid[uj][1] = true;         
#endif 
            }   
        }   
    }
}

void BoundCube::setBounds( const Point3D &p1, const Point3D &p2)
{
    for(unsigned int ui=0; ui<3; ui++)
    {
        bounds[ui][0]=std::min(p1[ui],p2[ui]);
        bounds[ui][1]=std::max(p1[ui],p2[ui]);
#ifdef DEBUG
        valid[ui][0]= true;
        valid[ui][1]= true;
#endif
    }

}


void BoundCube::setBounds(const Point3D &p1, float radius)
{
    for(unsigned int ui=0;ui<3;ui++)
    {
        bounds[ui][0]=p1[ui]-radius;
        bounds[ui][1]=p1[ui]+radius;
#ifdef DEBUG
        valid[ui][0]=true;  
        valid[ui][1]=true;  
#endif
    }
}

void BoundCube::setBounds(const std::vector<Point3D> &ptArray)
{
    setBounds(&(ptArray[0]),ptArray.size());
}

void BoundCube::setBounds(const std::vector<IonHit> &ptArray)
{
    setInverseLimits();
#pragma omp parallel for
    for(size_t ui=0;ui<ptArray.size();ui++)
    {
        for(unsigned int uj=0;uj<3;uj++)
        {
            bounds[uj][0]=std::min(bounds[uj][0],ptArray[ui].getPosRef()[uj]);
            bounds[uj][1]=std::max(bounds[uj][1],ptArray[ui].getPosRef()[uj]);
        }
    }

#ifdef DEBUG    
    for(unsigned int ui=0;ui<3;ui++)
    {
        if(bounds[ui][0] < bounds[ui][1])
            valid[ui][0]=valid[ui][1]=true;
    }
#endif
}


void BoundCube::getBounds(Point3D &low, Point3D &high) const 
{
    for(unsigned int ui=0; ui<3; ui++) 
    {
        ASSERT(valid[ui][0] && valid[ui][1]);
        low.setValue(ui,bounds[ui][0]);
        high.setValue(ui,bounds[ui][1]);
    }
}

Point3D BoundCube::min() const
{
    Point3D minPt;
    for(unsigned int ui=0; ui<3; ui++) 
    {
        ASSERT(valid[ui][0]);
        minPt.setValue(ui,bounds[ui][0]);
    }

    return minPt;
}

Point3D BoundCube::max() const
{
    Point3D maxPt;
    for(unsigned int ui=0; ui<3; ui++) 
    {
        ASSERT(valid[ui][1]);
        maxPt.setValue(ui,bounds[ui][1]);
    }

    return maxPt;
}
float BoundCube::getLargestDim() const
{
    float f;
    f=getSize(0);
    f=std::max(getSize(1),f);   
    return std::max(getSize(2),f);  
}

float BoundCube::getSize(unsigned int dim) const
{
    ASSERT(dim < 3);
#ifdef DEBUG
    for(unsigned int ui=0;ui<3; ui++)
    {
        ASSERT(valid[0][1] && valid [0][0]);
    }
#endif
    return fabs(bounds[dim][1] - bounds[dim][0]);
}

float BoundCube::getVolume() const
{
    return getSize(2)*getSize(1)*getSize(0);
}

bool BoundCube::contains(const BoundCube &b) const
{
    Point3D low,high;
    b.getBounds(low,high);
    return containsPt(low) && containsPt(high); 
}

bool BoundCube::containsPt(const Point3D &p) const
{
    for(unsigned int ui=0; ui<3; ui++) 
    {
        ASSERT(valid[ui][0] && valid[ui][1]);
        if( p[ui] < bounds[ui][0] || p[ui] > bounds[ui][1])
            return false;
    }

    return true;
}

bool BoundCube::containedInSphere(const Point3D &origin, float radius) const
{
    float sqrRad= radius*radius;
    //Check each coordinate for the box, (hi, low)*(x,y,z)
    for(size_t ui=0;ui<8;ui++)
    {
        Point3D p;
        p[0] = bounds[0][ui&1];
        p[1] = bounds[1][(ui&2) >> 1];
        p[2] = bounds[2][(ui&4)>> 2];

        if(p.sqrDist(origin) >=sqrRad)
            return false;
    }



    return true;
}

BoundCube BoundCube::makeUnion(const BoundCube &bC) const
{
    BoundCube res;
    for(unsigned int dim=0;dim<3;dim++)
    {
        float a,b;
        //Expand the lower bound down
        a=bounds[dim][0]; b=bC.bounds[dim][0];
        res.setBound(dim,0,std::min(a,b));
        //Expand the upper bound up
        a=bounds[dim][1]; b=bC.bounds[dim][1];
        res.setBound(dim,1,std::max(a,b));
    }

    return res;
}

BoundCube BoundCube::makeIntersection(const BoundCube &bC) const
{
    BoundCube res;
    for(unsigned int dim=0;dim<3;dim++)
    {
        float a,b;
        //Expand the lower bound down
        a=bounds[dim][0]; b=bC.bounds[dim][0];
        res.setBound(dim,0,std::max(a,b));
        //Expand the upper bound up
        a=bounds[dim][1]; b=bC.bounds[dim][1];
        res.setBound(dim,1,std::min(a,b));
    }

    return res;
}

//checks intersection with sphere [centre,centre+radius)
bool BoundCube::intersects(const Point3D &pt, float sqrRad) const
{
    Point3D nearPt;
    
    //Find the closest point on the cube  to the sphere
    unsigned int ui=2;
    const float *val=pt.getValueArr()+ui;
    do
    {
        if(*val <= bounds[ui][0])
        {
            nearPt.setValue(ui,bounds[ui][0]);
            --val;
            continue;
        }
        
        if(*val >=bounds[ui][1])
        {
            nearPt.setValue(ui,bounds[ui][1]);
            --val;
            continue;
        }
    
        nearPt.setValue(ui,*val);
        --val;
    }while(ui--);

    //now test the distance from nrPt to pt
    return (nearPt.sqrDist(pt) < sqrRad);
}

bool BoundCube::intersects(const BoundCube &b) const
{
    Point3D p1,p2;

    //Check to see if the other cube contains either of our bounds
    // If this is the case, we intersect
    getBounds(p1,p2);
    if(b.containsPt(p1) || b.containsPt(p2)) 
        return true;

    //Check to see if the other cube's bounds are contined in this one
    // this also implies intersection
    b.getBounds(p1,p2);
    if(containsPt(p1) || containsPt(p2)) 
        return true;

    //No containment.   
    return false;
}


Point3D BoundCube::getCentroid() const
{
#ifdef DEBUG
    for(unsigned int ui=0;ui<3; ui++)
    {
        ASSERT(valid[ui][1] && valid [ui][0]);
    }
#endif
    return Point3D(bounds[0][1] + bounds[0][0],
            bounds[1][1] + bounds[1][0],
            bounds[2][1] + bounds[2][0])/2.0f;
}

float BoundCube::getMaxDistanceToBox(const Point3D &queryPt) const
{
#ifdef DEBUG
    for(unsigned int ui=0;ui<3; ui++)
    {
        ASSERT(valid[ui][1] && valid [ui][0]);
    }
#endif

    float maxDistSqr=0.0f;

    //Could probably be a bit more compact.
    

    //Set lower and upper corners on the bounding rectangle
    Point3D p[2];
    p[0] = Point3D(bounds[0][0],bounds[1][0],bounds[2][0]);
    p[1] = Point3D(bounds[0][1],bounds[1][1],bounds[2][1]);

    //Count binary-wise selecting upper and lower limits, to enumerate all 8 verticies.
    for(unsigned int ui=0;ui<8; ui++)
    {
        maxDistSqr=std::max(maxDistSqr,
            queryPt.sqrDist(Point3D(p[ui&1][0],p[(ui&2) >> 1][1],p[(ui&4) >> 2][2])));
    }

    return sqrtf(maxDistSqr);
}


unsigned int BoundCube::segmentTriple(unsigned int dim, float slice) const
{
    ASSERT(isValid(dim));
    ASSERT(dim < 3);

    //check lower
    if( slice < bounds[dim][0])
        return 0;
    
    //check upper
    if( slice >=bounds[dim][1])
        return 2;

    return 1;

}

BoundCube BoundCube::operator=(const BoundCube &b)
{
    for(unsigned int ui=0;ui<3; ui++)
    {
        for(unsigned int uj=0;uj<2; uj++)
        {
#ifdef DEBUG
            valid[ui][uj] = b.valid[ui][uj];
#endif
            bounds[ui][uj] = b.bounds[ui][uj];

        }
    }

    return *this;
}

std::ostream &operator<<(std::ostream &stream, const BoundCube& b)
{
    stream << "Bounds :Low (";
    stream << b.bounds[0][0] << ",";
    stream << b.bounds[1][0] << ",";
    stream << b.bounds[2][0] << ") , High (";
    
    stream << b.bounds[0][1] << ",";
    stream << b.bounds[1][1] << ",";
    stream << b.bounds[2][1] << ")" << std::endl;

#ifdef DEBUG    
    
    stream << "Bounds Valid: Low (";
    stream << b.valid[0][0] << ",";
    stream << b.valid[1][0] << ",";
    stream << b.valid[2][0] << ") , High (";
    
    stream << b.valid[0][1] << ",";
    stream << b.valid[1][1] << ",";
    stream << b.valid[2][1] << ")" << std::endl;
#endif
    return stream;
}

#ifdef DEBUG
bool BoundCube::isValid() const
{
    for(unsigned int ui=0;ui<3; ui++)
    {
        if(!valid[ui][0] || !valid[ui][1])
            return false;
    }

    return true;
}

bool BoundCube::isValid(unsigned int dim) const
{
    return (valid[dim][0] && valid[dim][1]);
}

void BoundCube::setInvalid()
{
    valid[0][0]=false; valid[1][0]=false; valid[2][0]=false;
    valid[0][1]=false; valid[1][1]=false; valid[2][1]=false;
}
#endif
}