voxels.h

Go to the documentation of this file.

 /*
 * common/voxels.h - Voxelised data manipulation class
 * Copyright (C) 2018, D 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/>.
 */

#ifndef ATOMPROBE_VOXELS_H
#define ATOMPROBE_VOXELS_H

#include "atomprobe/helper/aptAssert.h"

#include "atomprobe/primitives/boundcube.h"
#include "atomprobe/primitives/point3D.h"
#include "atomprobe/primitives/ionHit.h"

#include <stack>
#include <numeric>
#include <atomic>
#include <algorithm>
#include <fstream>

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

#ifndef XOR
    #define XOR(a,b) ((!(a)) ^ (!(b)))
#endif

namespace AtomProbe{

enum{   
    BOUND_CLIP=1,
    BOUND_HOLD,
    BOUND_DERIV_HOLD,
    BOUND_MIRROR,
    BOUND_ZERO,
    BOUND_ENUM_END
};

//Interpolation mode for slice 
enum
{
    VOX_INTERP_NONE,
    VOX_INTERP_LINEAR,
    VOX_INTERP_ENUM_END
};

enum{
    ADJ_ALL,
    ADJ_PLUS
};

//Error codes
enum{
    VOXELS_BAD_FILE_READ=1,
    VOXELS_BAD_FILE_OPEN,
    VOXELS_BAD_FILE_SIZE,
    VOXELS_OUT_OF_MEMORY
};



enum {
    CLIP_NONE=0,
       CLIP_LOWER_SOUTH_WEST,
       CLIP_UPPER_NORTH_EAST
};


enum {
    VOXEL_ABORT_ERR,
    VOXEL_MEMORY_ERR,
    VOXEL_BOUNDS_INVALID_ERR
};

#ifdef DEBUG
    bool runVoxelTests();
#endif

//Voxel abortion flag. Set to true to initiate abort from supporting algorithms




template<class T> class Voxels
{
    private:
        size_t binCount[3];

        std::vector<T> voxels;
        
        //Dataset is assumed to sit in a rectilinear volume from minBound
        //to maxBound
        AtomProbe::Point3D minBound, maxBound;

        size_t shape3(size_t x,size_t y, size_t z) const;   
    public:

        Voxels();
        ~Voxels();

        void swap(Voxels<T> &v);

        void copy(Voxels<T> &newCopy) const;

        void setPoint(const AtomProbe::Point3D &pt, const T &val);

        T getPointData(const AtomProbe::Point3D &pt) const;


        void getIndex(size_t &x, size_t &y,
                size_t &z, const AtomProbe::Point3D &p) const;
        
        // including upper borders
        void getIndexWithUpper(size_t &x, size_t &y,
                size_t &z, const AtomProbe::Point3D &p) const;
        
        AtomProbe::Point3D getPoint(size_t x, 
                size_t y, size_t z) const;
        inline T getData(size_t x, size_t y, size_t z) const;
        inline T getData(size_t i) const { return voxels[i];}

        void setEntry(size_t n, const T &val) { voxels[n] = val;};
        //const T &getDataRef(size_t x, size_t y, size_t z) const;
        void setData(size_t x, size_t y, size_t z, const T &val);
        void setData(size_t n, const T &val);

        //Obtain an interpolated entry. The interpolated values are obtained by padding
        void getInterpolatedData(const AtomProbe::Point3D &pt, T &v) const;


        //get an interpolated slice from a section of the data
        void getInterpSlice(size_t normal, float offset, T *p, 
            size_t interpMode=VOX_INTERP_NONE) const;

        //Get a specific slice, from an integral offset in the data, no interp
        void getSlice(size_t normal, size_t offset, T *p) const;
        void getSize(size_t &x, size_t &y, size_t &z) const;


        size_t resize(size_t newX, size_t newY, size_t newZ, const AtomProbe::Point3D &newMinBound=AtomProbe::Point3D(0.0f,0.0f,0.0f), const AtomProbe::Point3D &newMaxBound=AtomProbe::Point3D(1.0f,1.0f,1.0f));
    
        size_t resize(const Voxels<T> &v);


        size_t resizeKeepData(size_t newX, size_t newY, size_t newZ, 
                    unsigned int direction=CLIP_LOWER_SOUTH_WEST, const AtomProbe::Point3D &newMinBound=AtomProbe::Point3D(0.0f,0.0f,0.0f), const AtomProbe::Point3D &newMaxBound=AtomProbe::Point3D(1.0f,1.0f,1.0f), const T &fill=T(0),bool doFill=false);




        size_t deprecatedGetEdgeUniqueIndex(size_t x,
            size_t y, size_t z, unsigned int edge) const;


        size_t getCellUniqueEdgeIndex(size_t x,
            size_t y, size_t z, unsigned int edge) const;
    
        // returns the axis value so you know edge vector too.
        // NOte that the value to pass as the edge index is (getEdgeIndex>>2)<<2 to
        // make the ownership of the voxel correct 
        void getEdgeEnds(size_t edgeIndex,AtomProbe::Point3D &a, AtomProbe::Point3D &b) const;


        void getEdgeCell(size_t edgeUniqId, size_t &x,size_t &y, size_t &z, size_t &axis ) const;

        //TODO: there is duplicate code between this and getEdgeEnds. Refactor.
        void getEdgeEndApproxVals(size_t edgeUniqId, T  &a, T  &b ) const;



        T getSum(const T &initialVal=T(0.0)) const;

        size_t count(const T &minIntensity) const;

        void fill(const T &val);
        AtomProbe::Point3D getMinBounds() const;
        AtomProbe::Point3D getMaxBounds() const;
        //Obtain the bounds for a specified axis
        void getAxisBounds(size_t axis, float &minV, float &maxV) const;
        AtomProbe::Point3D getPitch() const;
        void setBounds(const AtomProbe::Point3D &pMin, const AtomProbe::Point3D &pMax);
        void setBounds(const BoundCube &b);
        void getBounds(AtomProbe::Point3D &pMin, AtomProbe::Point3D &pMax) const { pMin=minBound;pMax=maxBound;}
        void getBounds(BoundCube &bc) const { bc.setBounds(minBound,maxBound);}

        size_t init(size_t nX,size_t nY,size_t nZ, const BoundCube &bound);
        size_t init(size_t nX,size_t nY, size_t nZ);

        size_t loadFile(const char *cpFilename, size_t nX,
                        size_t nY, size_t nZ, bool silent=false);

        size_t writeFile(const char *cpFilename) const;

        /* Data outside (either above if keepUpper is true, or below threshold if keepUpper is false)
        * will be set to the new value
        */
        void threshold(const T &thresh, bool keepUpper,const T & newVal); 
        
        /* Data below the threshold will be true if keepUpper is false, and vice-versa
        */
        void thresholdToBoolMask(const T &thresh, bool keepUpper,Voxels<bool> &result) const; 

        void applyMask(const Voxels<bool> &mask, const T &newVal, bool invert=false);


        void thresholdForPosition(std::vector<AtomProbe::Point3D> &p, const T &thresh, bool lowerEq=false) const;





        static size_t sizeofType() { return sizeof(T);}; 


        /* On thresholding (val > thresh), set the value to "onThresh". 
         * Otherwise set to "offthresh"
         */
        void binarise(Voxels<T> &result,const T &thresh, const T &onThresh, 
                const T &offThresh) const;
    

        /*Remove the data from the voxel set
         */
        void clear() { voxels.clear();};

        T min() const;

        T max() const;  
        
        void minMax(T &min, T &max) const;  



        int countPoints( const std::vector<AtomProbe::Point3D> &points, bool noWrap=true, bool doErase=true);
        

        int countIons( const std::vector<AtomProbe::IonHit> &ions, bool noWrap=true, bool doErase=true);

        T trapezIntegral() const;   
        // this is done by dividing each voxel by its volume 
        void calculateDensity();

        float getBinVolume() const;

        void operator/=(const Voxels<T> &v);

        void operator/=(const T &v);
        
        bool operator==(const Voxels<T> &v) const;

        size_t size() const { return voxels.size();}
};



template<class T, class U>
void castVoxels(const Voxels<T> &src, Voxels<U> &dest)
{
    //TODO:  Remove me!
    src=dest;
}

template<class T, class U>
void sumVoxels(const Voxels<T> &src, U &counter)
{
    size_t nx,ny,nz;
    counter=0;
    src.getSize(nx,ny,nz);

    size_t nMax=src.size();
    for(size_t ui=0; ui<nMax; ui++)
    {
        counter+=src.getData(ui);
    }

}

//====
template<class T>
Voxels<T>::Voxels() : voxels(), minBound(AtomProbe::Point3D(0,0,0)), maxBound(AtomProbe::Point3D(1,1,1))
{
}

template<class T>
Voxels<T>::~Voxels()
{
}


template<class T>
size_t Voxels<T>::shape3(size_t x, size_t y, size_t z) const
{
    return x + y*binCount[0] + z*binCount[0]*binCount[1];
}

template<class T>
void Voxels<T>::copy(Voxels<T> &newVox) const
{
    newVox.binCount[0]=binCount[0];
    newVox.binCount[1]=binCount[1];
    newVox.binCount[2]=binCount[2];

    newVox.voxels=voxels;
    newVox.minBound=minBound;
    newVox.maxBound=maxBound;


}

template<class T>
void Voxels<T>::setPoint(const AtomProbe::Point3D &point,const T &val)
{
    ASSERT(voxels.size());
    size_t pos[3];
    for(size_t ui=0;ui<3;ui++)
        pos[ui] = (size_t)round(point[ui]*(float)binCount[ui]);

    voxels[shape3(pos[0],pos[1],pos[2])] = val;
}


template<class T>
void Voxels<T>::setData(size_t x, size_t y, 
            size_t z, const T &val)
{
    ASSERT(voxels.size());

    ASSERT( x < binCount[0] && y < binCount[1] && z < binCount[2]);
    voxels[shape3(x,y,z)]=val;
}

template<class T>
inline void Voxels<T>::setData(size_t n, const T &val)
{
    ASSERT(voxels.size());
    ASSERT(n<voxels.size());

    voxels[n] =val; 
}

template<class T>
T Voxels<T>::getPointData(const AtomProbe::Point3D &point) const
{
    ASSERT(voxels.size());
    size_t pos[3];
    AtomProbe::Point3D offsetFrac;
    offsetFrac=point-minBound;
    for(size_t ui=0;ui<3;ui++)
    {
        offsetFrac[ui]/=(maxBound[ui]-minBound[ui]);
        pos[ui] = (size_t)round(offsetFrac[ui]*(float)binCount[ui]);
    }   

    return voxels[shape3(pos[0],pos[1],pos[2])];
}

template<class T>
AtomProbe::Point3D Voxels<T>::getPoint(size_t x, size_t y, size_t z) const
{
    //ASSERT(x < binCount[0] && y<binCount[1] && z<binCount[2]);

    return AtomProbe::Point3D((float)x/(float)binCount[0]*(maxBound[0]-minBound[0]) + minBound[0],
            (float)y/(float)binCount[1]*(maxBound[1]-minBound[1]) + minBound[1],
            (float)z/(float)binCount[2]*(maxBound[2]-minBound[2]) + minBound[2]);
}

template<class T>
AtomProbe::Point3D Voxels<T>::getPitch() const
{
    return AtomProbe::Point3D((float)1.0/(float)binCount[0]*(maxBound[0]-minBound[0]),
            (float)1.0/(float)binCount[1]*(maxBound[1]-minBound[1]),
            (float)1.0/(float)binCount[2]*(maxBound[2]-minBound[2]));
}

template<class T>
void Voxels<T>::getSize(size_t &x, size_t &y, size_t &z) const
{
    ASSERT(voxels.size());
    x=binCount[0];
    y=binCount[1];
    z=binCount[2];
}

template<class T>
size_t Voxels<T>::deprecatedGetEdgeUniqueIndex(size_t x,size_t y, size_t z, unsigned int index) const
{
    //This provides a reversible mapping of x,y,z 
    //X aligned edges are first
    //Y second
    //Z third
    

    //Consider each parallel set of edges (eg all the X aligned edges)
    //to be the dual grid of the actual grid. From this you can visualise the
    //cell centres moving -1/2 -/12 units in the direction normal to the edge direction
    //to produce the centres of the edge. An additional vertex needs to be created at
    //the end of each dimension not equal to the alignement dim.


    //              ->ASCII ART TIME<-
    // In each individual cube, the offsets look like this:
        //      ------------7-----------
        //      \              |\ .
        //      |\             | \ .
        //      | 10               |  11
        //      |  \               |   \ .
        //      |   \              |    \ .
        //      |    \ --------6-------------|
        //      |     |                |     |
        //        2 |                3     |     |
        //      |     |                |     |
        //      |     |                |     |
        //      |     |                |     |
        //      |     0                |     1
        //      \-----|----5-----------      |
        //       \    |                 \    |
        //        8   |                  9   |
        //         \  |                   \  |
        //          \ |---------4------------
    //
    //      ^x
    //  z|\ |
    //    \ |
    //     \-->y
    //

    switch(index)
    {
        //X axis aligned
        //--
        case 0:
            break;  
        case 1:
            y++; // one across in Y
            break;  
        case 2:
            z++;//One across in Z
            break;  
        case 3:
            y++;
            z++;
            break;  
        //--

        //Y Axis aligned
        //--
        case 4:
            break;  
        case 5:
            z++;
            break;  
        case 6:
            x++;
            break;  
        case 7:
            z++;
            x++;
            break;  
        //--
        
        //Z Axis aligned
        //--
        case 8:
            break;  
        case 9:
            y++;
            break;  
        case 10:
            x++;
            break;  
        case 11:
            x++;
            y++;
            break;  
        //--
    }


    size_t result = 12*(z + y*(binCount[2]+1) + x*(binCount[2]+1)*(binCount[1]+1)) +
    index;
    
    return result;

}

/*
template<class T>
size_t Voxels<T>::getEdgeUniqueIndex(size_t x,size_t y, size_t z, unsigned int index) const
{
    //This provides a reversible mapping of x,y,z 
    //X aligned edges are first
    //Y second
    //Z third
    

    //Consider each parallel set of edges (eg all the X aligned edges)
    //to be the dual grid of the actual grid. From this you can visualise the
    //cell centres moving -1/2 -/12 units in the direction normal to the edge direction
    //to produce the centres of the edge. An additional vertex needs to be created at
    //the end of each dimension not equal to the alignement dim.


    //              ->ASCII ART TIME<-
    // In each individual cube, the offsets look like this:
        //      ------------7-----------
        //      \              |\ .
        //      |\             | \ .
        //      | 10               |  11
        //      |  \               |   \ .
        //      |   \              |    \ .
        //      |    \ --------6-------------|
        //      |     |                |     |
        //        2 |                3     |     |
        //      |     |                |     |
        //      |     |                |     |
        //      |     |                |     |
        //      |     0                |     1
        //      \-----|----5-----------      |
        //       \    |                 \    |
        //        8   |                  9   |
        //         \  |                   \  |
        //          \ |---------4------------
    //
    //      ^x
    //  z|\ |
    //    \ |
    //     \-->y
    //

    switch(index)
    {
        //X axis aligned
        //--
        case 0:
            break;  
        case 1:
            y++; // one across in Y
            break;  
        case 2:
            z++;//One across in Z
            break;  
        case 3:
            y++;
            z++;
            break;  
        //--

        //Y Axis aligned
        //--
        case 4:
            break;  
        case 5:
            z++;
            break;  
        case 6:
            x++;
            break;  
        case 7:
            z++;
            x++;
            break;  
        //--
        
        //Z Axis aligned
        //--
        case 8:
            break;  
        case 9:
            y++;
            break;  
        case 10:
            x++;
            break;  
        case 11:
            x++;
            y++;
            break;  
        //--
    }

    unsigned int axis = index/4;
    size_t result = 3*(z + y*(binCount[2]+1) + x*(binCount[2]+1)*(binCount[1]+1)) + axis;
    
    return result;

}
*/
template<class T>
size_t Voxels<T>::getCellUniqueEdgeIndex(size_t x,size_t y, size_t z, unsigned int index) const
{
    //This provides a reversible mapping of x,y,z 
    //X aligned edges are first
    //Y second
    //Z third
    

    //Consider each parallel set of edges (eg all the X aligned edges)
    //to be the dual grid of the actual grid. From this you can visualise the
    //cell centres moving -1/2 -/12 units in the direction normal to the edge direction
    //to produce the centres of the edge. An additional vertex needs to be created at
    //the end of each dimension not equal to the alignement dim.


    //              ->ASCII ART TIME<-
    // In each individual cube, the offsets look like this:
        //      ------------7-----------
        //      \              |\ .
        //      |\             | \ .
        //      | 10               |  11
        //      |  \               |   \ .
        //      |   \              |    \ .
        //      |    \ --------6-------------|
        //      |     |                |     |
        //        2 |                3     |     |
        //      |     |                |     |
        //      |     |                |     |
        //      |     |                |     |
        //      |     0                |     1
        //      \-----|----5-----------      |
        //       \    |                 \    |
        //        8   |                  9   |
        //         \  |                   \  |
        //          \ |---------4------------
    //
    //      ^x
    //  z|\ |
    //    \ |
    //     \-->y
    //

    size_t cellIdx = 12*(z + y*(binCount[2]+1) + x*(binCount[2]+1)*(binCount[1]+1)) ;

    cellIdx+=index;
    return cellIdx;

}

template<class T>
void Voxels<T>::getEdgeCell(size_t edgeUniqId, size_t &x,size_t &y, size_t &z, size_t &axis ) const
{
    //Invert the mapping generated by the edgeUniqId function
    //to retrieve the XYZ and axis values
    //--
    axis=(edgeUniqId%12)/4;

    size_t tmp = edgeUniqId/12;
    x = tmp/((binCount[2]+1)*(binCount[1]+1));
    tmp-=x*((binCount[2]+1)*(binCount[1]+1));

    y=tmp/(binCount[2]+1);
    tmp-=y*(binCount[2]+1);

    z=tmp;
    //--

    ASSERT(x< binCount[0]+1 && y<binCount[1]+1 && z<binCount[2]+1);
}
template<class T>
void Voxels<T>::getEdgeEnds(size_t edgeUniqId, AtomProbe::Point3D &a, AtomProbe::Point3D &b ) const
{
    size_t x,y,z;
    size_t axis;
    getEdgeCell(edgeUniqId,x,y,z,axis);

    AtomProbe::Point3D delta=getPitch();
    AtomProbe::Point3D cellCentre=getPoint(x,y,z);

    //Generate ends of edge, as seen in ascii diagram in uniqueID
    switch(axis)
    {
        case 0:
            //|| x
            a=cellCentre;
            b=cellCentre + AtomProbe::Point3D(delta[0],0,0);
            break;

        case 1:
            //|| y
            a=cellCentre;
            b=cellCentre + AtomProbe::Point3D(0,delta[1],0);
            break;
        case 2:
            //|| z
            a=cellCentre; 
            b=cellCentre + AtomProbe::Point3D(0,0,delta[2]);
            break;
        default:
            ASSERT(false);
    }


#ifdef DEBUG
    BoundCube bc;
    bc.setBounds(getMinBounds(), getMaxBounds());
    bc.expand(sqrtf(std::numeric_limits<float>::epsilon()));
    ASSERT(bc.containsPt(a) && bc.containsPt(b));
#endif
}

template<class T>
void Voxels<T>::getEdgeEndApproxVals(size_t edgeUniqId, T  &a, T  &b ) const
{

    //TODO: Speed me up? Could not use
    //  other routines to do this access.
    //  I think some redundant calculations are done
    AtomProbe::Point3D ptA,ptB;
    getEdgeEnds(edgeUniqId,ptA,ptB);
    size_t x,y,z;
    getIndex(x,y,z,ptA);    
    a = getPointData(ptA);
    b=getPointData(ptB);
}

template<class T>
void Voxels<T>::getAxisBounds(size_t axis, float &minV, float &maxV ) const
{
    maxV=maxBound[axis];
    minV=minBound[axis];
}

template<class T>
size_t Voxels<T>::resize(size_t x, size_t y, size_t z, const AtomProbe::Point3D &newMinBound, const AtomProbe::Point3D &newMaxBound) 
{
    binCount[0] = x;
    binCount[1] = y;
    binCount[2] = z;


    minBound=newMinBound;
    maxBound=newMaxBound;

    try
    {
        voxels.resize(binCount[0]*binCount[1]*binCount[2]);
    }
    catch(std::exception &e)
    {
        return 1;
    }

    return 0;
}

template<class T>
size_t Voxels<T>::resize(const Voxels<T> &oth)
{
    return resize(oth.binCount[0],oth.binCount[1],oth.binCount[2],oth.minBound,oth.maxBound);
}

template<class T>
size_t Voxels<T>::resizeKeepData(size_t newX, size_t newY, size_t newZ, 
            unsigned int direction, const AtomProbe::Point3D &newMinBound, const AtomProbe::Point3D &newMaxBound, const T &fill,bool doFill)
{

    ASSERT(direction==CLIP_LOWER_SOUTH_WEST);

    Voxels<T> v;
    
    if(v.resize(newX,newY,newZ))
        return 1;

    switch(direction)
    {
        case CLIP_LOWER_SOUTH_WEST:
        {
            minBound=newMinBound;
            maxBound=newMaxBound;

            if(doFill)
            {
                size_t itStop[3];
                itStop[0]=std::min(newX,binCount[0]);
                itStop[1]=std::min(newY,binCount[1]);
                itStop[2]=std::min(newZ,binCount[2]);

                size_t itMax[3];
                itMax[0]=std::max(newX,binCount[0]);
                itMax[1]=std::max(newY,binCount[1]);
                itMax[2]=std::max(newZ,binCount[2]);
                //Duplicate into new value, if currently inside bounding box
                //This logic will be a bit slow, owing to repeated "if"ing, but
                //it is easy to understand. Other logics would have many more
                //corner cases
#pragma omp parallel for
                for(size_t ui=0;ui<itMax[0];ui++)
                {

                    for(size_t uj=0;uj<itMax[1];uj++)
                    {
                        for(size_t uk=0;uk<itMax[2];uk++)
                        {
                            if(itStop[0]< binCount[0] && 
                                itStop[1]<binCount[1] && itStop[2] < binCount[2])
                                v.setData(ui,uj,uk,getData(ui,uj,uk));
                            else
                                v.setData(ui,uj,uk,fill);
                        }
                    }

                }
            }
            else
            {
                //Duplicate into new value, if currently inside bounding box
#pragma omp parallel for
                for(size_t ui=0;ui<newX;ui++)
                {

                    for(size_t uj=0;uj<newY;uj++)
                    {
                        for(size_t uk=0;uk<newZ;uk++)
                            v.setData(ui,uj,uk,getData(ui,uj,uk));
                    }

                }

            }



            break;
        }

        default:
            //Not implemented
            ASSERT(false);
    }

    swap(v);
    return 0;
}

template<class T>
AtomProbe::Point3D Voxels<T>::getMinBounds() const
{
    ASSERT(voxels.size());
    return minBound;
}
                                         
template<class T>
AtomProbe::Point3D Voxels<T>::getMaxBounds() const
{
    ASSERT(voxels.size());
    return maxBound;
}
                                         
template<class T>
void Voxels<T>::setBounds(const AtomProbe::Point3D &pMin, const AtomProbe::Point3D &pMax)
{
    ASSERT(voxels.size());
#ifdef DEBUG
    for(auto ui=0;ui<3;ui++)
    {
        ASSERT(pMin[ui] <=pMax[ui]);
    }
#endif

    minBound=pMin;
    maxBound=pMax;
}

    template<class T>
void Voxels<T>::setBounds(const BoundCube &bc)
{
    setBounds(bc.min(),bc.max());
}

template<class T>
size_t Voxels<T>::init(size_t nX, size_t nY,
                        size_t nZ, const BoundCube &bound)
{
    binCount[0]=nX;
    binCount[1]=nY;
    binCount[2]=nZ;

    bound.getBounds(minBound, maxBound);

    voxels.resize(nX*nY*nZ);

    return 0;
}

template<class T>
size_t Voxels<T>::init(size_t nX, size_t nY, size_t nZ)

{
    AtomProbe::Point3D pMin(0,0,0), pMax(nX,nY,nZ); 

    return init(nX,nY,nZ,AtomProbe::BoundCube(pMin,pMax));
}

template<class T>
size_t Voxels<T>::loadFile(const char *cpFilename, size_t nX, size_t nY, size_t nZ , bool silent)
{
    //Must be power of two (buffer size when loading files, in sizeof(T)s)
    const unsigned int ITEM_BUFFER_SIZE=65536;
    
    std::ifstream CFile(cpFilename,std::ios::binary);

    if(!CFile)
        return VOXELS_BAD_FILE_OPEN;
    
    CFile.seekg(0,std::ios::end);
    
    
    size_t fileSize = CFile.tellg();
    if(fileSize !=nX*nY*nZ*sizeof(T))
        return VOXELS_BAD_FILE_SIZE;

    resize(nX,nY,nZ,AtomProbe::Point3D(nX,nY,nZ));
    
    CFile.seekg(0,std::ios::beg);

    unsigned int curBufferSize=ITEM_BUFFER_SIZE*sizeof(T);
    unsigned char *buffer = new unsigned char[curBufferSize];

    //Shrink the buffer size by factors of 2
    //in the case of small datasets
    while(fileSize < curBufferSize)
        curBufferSize = curBufferSize >> 1;

    
    //Draw a progress bar
    if(!silent)
    {
        std::cerr << std::endl << "|";
        for(unsigned int ui=0; ui<100; ui++)
            std::cerr << ".";
        std::cerr << "| 100%" << std::endl << "|";
    }
        
    unsigned int lastFrac=0;
    unsigned int ui=0;
    unsigned int pts=0;
    do
    {
    
        //Still have data? Keep going   
        while((size_t)CFile.tellg() <= fileSize-curBufferSize)
        {
            //Update progress bar
            if(!silent && ((unsigned int)(((float)CFile.tellg()*100.0f)/(float)fileSize) > lastFrac))
            {
                std::cerr << ".";
                pts++;
                lastFrac=(unsigned int)(((float)CFile.tellg()*100.0f)/(float)fileSize) ;    
            }

            //Read a chunk from the file
            CFile.read((char *)buffer,curBufferSize);
    
            if(!CFile.good())
                return VOXELS_BAD_FILE_READ;

            //Place the chunk contents into ram
            for(size_t position=0; position<curBufferSize; position+=(sizeof(T)))
                voxels[ui++] = (*((T *)(buffer+position)));
        }
                
        //Halve the buffer size
        curBufferSize = curBufferSize >> 1 ;

    }while(curBufferSize> sizeof(T)); //This does a few extra loops. Not many 

    delete[] buffer;

    //Fill out the progress bar
    if(!silent)
    {
        while(pts++ <100)
            std::cerr << ".";
    
        std::cerr << "| done" << std::endl;
    }

    return 0;
}

template<class T>
size_t Voxels<T>::writeFile(const char *filename) const
{

    ASSERT(voxels.size());

    std::ofstream file(filename, std::ios::binary);

    if(!file)
        return 1;

    
    for(size_t ui=0; ui<voxels.size(); ui++)
    {
        T v;
        v=voxels[ui];
        file.write((char *)&v,sizeof(T));
        if(!file.good())
            return 2;
    }
    return 0;
}


template<class T>
T Voxels<T>::getSum(const T &initialValue) const
{
    ASSERT(voxels.size());
    T tmp(initialValue);
    size_t n=voxels.size();
#pragma omp parallel for reduction(+:tmp)
    for(size_t ui=0;ui<n;ui++)
        tmp+=voxels[ui];

    return tmp;
}

template<class T>
T Voxels<T>::trapezIntegral() const
{
    //Compute volume prefactor - volume of cube of each voxel
    //--
    float prefactor=1.0;
    for(size_t ui=0;ui<3;ui++)
    {
        prefactor*=(maxBound[ui]-
            minBound[ui])/
            (float)binCount[ui];
    }

    //--


    double accumulation(0.0);
    //Loop across dataset integrating along z direction
#pragma omp parallel for reduction(+:accumulation)
    for(size_t ui=0;ui<voxels.size(); ui++)
        accumulation+=voxels[ui];

    return prefactor*accumulation;
}


template<class T>
size_t Voxels<T>::count(const T &minIntensity) const
{
    size_t bins;
    bins=binCount[0]*binCount[1]*binCount[2];

    size_t sum=0;
#pragma omp parallel for reduction(+:sum)
    for(size_t ui=0;ui<bins; ui++)
    {
        if(voxels[ui]>=minIntensity)
            sum++;
    }

    return sum;
}

template<class T>
void Voxels<T>::swap(Voxels<T> &other)
{
    std::swap(binCount[0],other.binCount[0]);
    std::swap(binCount[1],other.binCount[1]);
    std::swap(binCount[2],other.binCount[2]);

    voxels.swap(other.voxels);
    
    std::swap(maxBound,other.maxBound);
    std::swap(minBound,other.minBound);
}

template<class T>
T Voxels<T>::getData(size_t x, size_t y, size_t z) const
{
    ASSERT(x < binCount[0] && y < binCount[1] && z < binCount[2]);
    return  voxels[shape3(x,y,z)]; 
}

template<class T>
void Voxels<T>::thresholdForPosition(std::vector<AtomProbe::Point3D> &p, const T &thresh, bool lowerEq) const
{
    p.clear();

    if(lowerEq)
    {
#pragma omp parallel for 
        for(size_t ui=0;ui<binCount[0]; ui++)
        {
            for(size_t uj=0;uj<binCount[1]; uj++)
            {
                for(size_t uk=0;uk<binCount[2]; uk++)
                {
                    if( getData(ui,uj,uk) < thresh)
                    {
#pragma omp critical
                        p.push_back(getPoint(ui,uj,uk));
                    }

                }
            }
        }
    }
    else
    {
#pragma omp parallel for 
        for(size_t ui=0;ui<binCount[0]; ui++)
        {
            for(size_t uj=0;uj<binCount[1]; uj++)
            {
                for(size_t uk=0;uk<binCount[2]; uk++)
                {
                    if( getData(ui,uj,uk) > thresh)
                    {
#pragma omp critical
                        p.push_back(getPoint(ui,uj,uk));
                    }

                }
            }
        }
    }
}

template<class T>
void Voxels<T>::threshold(const T &thresh, bool keepUpper, const T &newVal)
{

    if(keepUpper)
    {
#pragma omp parallel for
        for(size_t ui=0;ui<(size_t)binCount[0]; ui++)
        {
            for(size_t uj=0;uj<binCount[1]; uj++)
            {
                for(size_t uk=0;uk<binCount[2]; uk++)
                {
                    if( getData(ui,uj,uk) < thresh)
                        setData(ui,uj,uk,newVal);
                }
            }
        }
    }
    else
    {
#pragma omp parallel for
        for(size_t ui=0;ui<(size_t)binCount[0]; ui++)
        {
            for(size_t uj=0;uj<binCount[1]; uj++)
            {
                for(size_t uk=0;uk<binCount[2]; uk++)
                {
                    if( getData(ui,uj,uk) > thresh)
                        setData(ui,uj,uk,newVal);
                }
            }
        }
    }
}

template<class T>
void Voxels<T>::thresholdToBoolMask(const T &thresh, bool keepUpper, Voxels<bool> &result) const
{
    result.resize(binCount[0],binCount[1],binCount[2],
                minBound,maxBound);
    
#pragma omp parallel for
    for(size_t ui=0;ui<(size_t)binCount[0]; ui++)
    {
        for(size_t uj=0;uj<binCount[1]; uj++)
        {
            for(size_t uk=0;uk<binCount[2]; uk++)
            {
                bool res;
                res=(getData(ui,uj,uk) > thresh); 
                
                result.setData(ui,uj,uk,XOR(keepUpper,res));
            }
        }
    }

}

template<class T>
void Voxels<T>::applyMask(const Voxels<bool> &mask, const T&newVal,bool invert)
{
    ASSERT(mask.size() == size());
    ASSERT(mask.getMinBounds()== minBound);
    ASSERT(mask.getMaxBounds() == maxBound);
#ifdef DEBUG
    size_t v[3];
    mask.getSize(v[0],v[1],v[2]);
    for(auto i=0;i<3;i++)
    {
        ASSERT(v[i] == binCount[i]);
    }
#endif

    if(invert)
    {
        for(auto i=0; i<voxels.size(); i++)
        {
            if(!mask.getData(i))
                voxels[i]=newVal;
        }
    }
    else
    {
        for(auto i=0; i<voxels.size(); i++)
        {
            if(mask.getData(i))
                voxels[i]=newVal;
        }
    }
}

template<class T>
void Voxels<T>::binarise(Voxels<T> &result, const T &thresh, 
                const T &onThresh, const T &offThresh) const
{

    result.resize(binCount[0],binCount[1],
            binCount[2],minBound,maxBound);
#pragma omp parallel for
    for(size_t ui=0;ui<(size_t)binCount[0]; ui++)
    {
        for(size_t uj=0;uj<binCount[1]; uj++)
        {
            for(size_t uk=0;uk<binCount[2]; uk++)
            {
                if( getData(ui,uj,uk) < thresh)
                    result.setData(ui,uj,uk,offThresh);
                else
                {
                    result.setData(ui,uj,uk,onThresh);
                }
            }
        }
    }
}



template<class T>
T Voxels<T>::min() const
{
    ASSERT(voxels.size());
    return *std::min_element(voxels.begin(),voxels.end());
}

template<class T>
T Voxels<T>::max() const
{
    ASSERT(voxels.size());
    return *std::max_element(voxels.begin(),voxels.end());
}


template<class T>
void Voxels<T>::minMax(T &min, T&max) const
{
    ASSERT(voxels.size());
    //FIXME: we can probably do this in one traverse, maybe quicker?
    // needs to be benchmarked.
    min=*std::min_element(voxels.begin(),voxels.end());
    max=*std::max_element(voxels.begin(),voxels.end());
}

template<class T>
int Voxels<T>::countPoints( const std::vector<AtomProbe::Point3D> &points, bool noWrap, bool doErase)
{
    if(doErase)
    {
        fill(0);    
    }

    size_t x,y,z;

    for(size_t ui=0; ui<points.size(); ui++)
    {
        
        getIndex(x,y,z,points[ui]);

        //Ensure it lies within the dataset 
        if(x < binCount[0] && y < binCount[1] && z< binCount[2])
        {
            {
                T value;
                value=getData(x,y,z)+T(1);

                //Prevent wrap-around errors
                if (noWrap) {
                    if(value > getData(x,y,z))
                        setData(x,y,z,value);
                } else {
                    setData(x,y,z,value);
                }
            }
        }   
    }

    return 0;
}

template<class T>
int Voxels<T>::countIons( const std::vector<AtomProbe::IonHit> &points, bool noWrap, bool doErase)
{
    if(doErase)
    {
        fill(0);    
    }

    size_t x,y,z;

    if (noWrap) 
    {
        for(size_t ui=0; ui<points.size(); ui++)
        {

            getIndex(x,y,z,points[ui].getPosRef());

            //Ensure it lies within the dataset
            if(x < binCount[0] && y < binCount[1] && z< binCount[2])
            {
                {
                    T value;
                    value=getData(x,y,z)+T(1);

                    //Prevent wrap-around errors
                    if(value > getData(x,y,z))
                        setData(x,y,z,value);
                }
            }
        }
    } 
    else 
    {
        for(size_t ui=0; ui<points.size(); ui++)
        {

            getIndex(x,y,z,points[ui].getPosRef());

            //Ensure it lies within the dataset
            if(x < binCount[0] && y < binCount[1] && z< binCount[2])
            {
                T value;
                value=getData(x,y,z)+T(1);
                setData(x,y,z,value);
            }
        }
    }

    return 0;
}

template<class T>
void Voxels<T>::calculateDensity()
{
    AtomProbe::Point3D size = maxBound - minBound;
    // calculate the volume of a voxel
    double volume = 1.0;
    for (int i = 0; i < 3; i++)
        volume *= size[i] / binCount[i];
    
    // normalise the voxel value based on volume
#pragma omp parallel for
    for(size_t ui=0; ui<voxels.size(); ui++) 
        voxels[ui] /= volume;

}

template<class T>
float Voxels<T>::getBinVolume() const
{
    AtomProbe::Point3D size = maxBound - minBound;
    double volume = 1.0;
    for (int i = 0; i < 3; i++)
        volume *= size[i] / binCount[i];

    return volume;
}

template<class T>
void Voxels<T>::getIndex(size_t &x, size_t &y,
                    size_t &z, const AtomProbe::Point3D &p) const
{

    ASSERT(p[0] >=minBound[0] && p[1] >=minBound[1] && p[2] >=minBound[2] &&
           p[0] <=maxBound[0] && p[1] <=maxBound[1] && p[2] <=maxBound[2]);
    x=(size_t)((p[0]-minBound[0])/(maxBound[0]-minBound[0])*(float)binCount[0]);
    y=(size_t)((p[1]-minBound[1])/(maxBound[1]-minBound[1])*(float)binCount[1]);
    z=(size_t)((p[2]-minBound[2])/(maxBound[2]-minBound[2])*(float)binCount[2]);

    if(x == binCount[0])
        x--;
    if(y == binCount[1])
        y--;
    if(z == binCount[2])
        z--;
}

template<class T>
void Voxels<T>::getIndexWithUpper(size_t &x, size_t &y,
                    size_t &z, const AtomProbe::Point3D &p) const
{
    //Get the data index as per normal
    getIndex(x,y,z,p);

    //but, as a special case, if the index is the same as our bincount, check
    //to see if it is positioned on an edge
    if(x==binCount[0] &&  
        fabs(p[0] -maxBound[0]) < sqrtf(std::numeric_limits<float>::epsilon()))
        x--;
    if(y==binCount[1] &&  
        fabs(p[1] -maxBound[1]) < sqrtf(std::numeric_limits<float>::epsilon()))
        y--;
    if(z==binCount[2] &&  
        fabs(p[2] -maxBound[2]) < sqrtf(std::numeric_limits<float>::epsilon()))
        z--;

}

template<class T>
void Voxels<T>::fill(const T &v)
{
    std::fill(voxels.begin(),voxels.end(),v);
}

//Obtain a slice of the voxel data. Data output will be in column order
// p[posB*nA + posA]. Input slice must be sufficiently sized and allocated
// to hold the output data
template<class T>
void Voxels<T>::getSlice(size_t normalAxis, size_t offset, T *p) const
{
    ASSERT(normalAxis < 3);

    size_t dimA,dimB,nA;
    switch(normalAxis)
    {
        case 0:
        {
            dimA=1;
            dimB=2;
            nA=binCount[dimA];
            break;
        }
        case 1:
        {   
            dimA=0;
            dimB=2;
            nA=binCount[dimA];
            break;
        }
        case 2:
        {
            dimA=0;
            dimB=1;
            nA=binCount[dimA];
            break;
        }
        default:
            ASSERT(false); //WTF - how did you get here??
    }
        

    //We are within bounds, use normal access functions
    switch(normalAxis)
    {
        case 0:
        {
            for(size_t ui=0;ui<binCount[dimA];ui++)
            {
                for(size_t uj=0;uj<binCount[dimB];uj++)
                    p[uj*nA + ui] = getData(offset,ui,uj);
            }
            break;
        }
        case 1:
        {   
            for(size_t ui=0;ui<binCount[dimA];ui++)
            {
                for(size_t uj=0;uj<binCount[dimB];uj++)
                    p[uj*nA + ui] = getData(ui,offset,uj);
            }
            break;
        }
        case 2:
        {
            for(size_t ui=0;ui<binCount[dimA];ui++)
            {
                for(size_t uj=0;uj<binCount[dimB];uj++)
                    p[uj*nA + ui] = getData(ui,uj,offset);
            }
            break;
        }
        default:
            ASSERT(false);
    }
}

template<class T>
void Voxels<T>::getInterpSlice(size_t normal, float offset, 
        T *p, size_t interpMode) const
{
    ASSERT(offset <=1.0f && offset >=0.0f);

    //Obtain the appropriately interpolated slice
    switch(interpMode)
    {
        case VOX_INTERP_NONE:
        {
            size_t slicePos;
            slicePos=roundf(offset*binCount[normal]);
            slicePos=std::min(slicePos,binCount[normal]-1);
            getSlice(normal,slicePos,p);
            break;
        }
        case VOX_INTERP_LINEAR:
        {
            //Find the upper and lower bounds, then
            // limit them so we don't fall off the end of the dataset
            size_t sliceUpper,sliceLower;
            if(binCount[0] == 1)
                sliceUpper=sliceLower=0;
            else
            {
                sliceUpper=ceilf(offset*binCount[normal]);

                if(sliceUpper >=binCount[normal])
                    sliceUpper=binCount[normal]-1;
                else if(sliceUpper==0)
                    sliceUpper=1;
                
                sliceLower=sliceUpper-1;
            }

            {
                T *pLower;
                size_t numEntries=binCount[(normal+1)%3]*binCount[(normal+2)%3];
                
                pLower  = new T[numEntries];

                getSlice(normal,sliceLower,pLower);
                getSlice(normal,sliceUpper,p);

                //Get the decimal part of the float
                float integ;
                float delta=modff(offset*binCount[normal],&integ);
                for(size_t ui=0;ui<numEntries;ui++)
                    p[ui] = delta*(p[ui]-pLower[ui]) + pLower[ui];

                delete[] pLower;
            }
            break;
        }
        default:
            ASSERT(false);
            
    }

}

//FIXME: I think this has a slight shift as we are moving the data voxels
// definition of the voxel centre by 1/2 a pitch, I think
template<class T>
void Voxels<T>::getInterpolatedData(const AtomProbe::Point3D &p, T &v) const

{
#ifdef DEBUG
    BoundCube bc(minBound,maxBound);
    ASSERT(bc.containsPt(p));
#endif

    size_t index[3];
    getIndex(index[0],index[1],index[2],p);

    AtomProbe::Point3D pitch =getPitch();
    
    //Find the offset to the voxel that we are in.
    //fraction should be in range [0,1)
    AtomProbe::Point3D fraction = p - (minBound + AtomProbe::Point3D(index[0],index[1],index[2])*pitch);
    fraction =fraction/pitch;


    size_t iPlus[3];
    //0.5 corresponds to voxel centre.
    for(unsigned int ui=0;ui<3;ui++)
    {
        if(index[ui] == (binCount[ui]-1))
            iPlus[ui]=0;
        else
            iPlus[ui]=1;
    }   
    

    float c[2][2];
    //Tri-linear interpolation

    //interpolate data values at cube vertices that surround point. We are coming from below the point
    // so we are simply extending the field on the upper edge by duplicating values as needed

    float xf = 1-fraction[0];
    
    size_t xLow,xHigh;
    xLow = index[0];
    xHigh = index[0] + iPlus[0];
    c[0][0] = getData(xLow,index[1],index[2])*xf            + getData(xHigh,index[1],index[2])*fraction[0] ;
    c[0][1] = getData(xLow,index[1],index[2]+iPlus[2])*xf       + getData(xHigh,index[1],index[2]+iPlus[2])*fraction[0] ;
    c[1][0] = getData(xLow,index[1]+iPlus[1],index[2])*xf       + getData(xHigh,index[1]+iPlus[1],index[2])*fraction[0] ;
    c[1][1] = getData(xLow,index[1]+iPlus[1],index[2]+iPlus[2])*xf  + getData(xHigh,index[1]+iPlus[1],index[2]+iPlus[2])*fraction[0] ;


    float c0,c1;
    c0 = c[0][0]*(1-fraction[1]) + c[1][0]*fraction[1];
    c1 = c[0][1]*(1-fraction[1]) + c[1][1]*fraction[1];

    v= c0*(1-fraction[2])+c1*fraction[2];   
    
}   



template<class T>
void Voxels<T>::operator/=(const Voxels<T> &v)
{
    ASSERT(v.voxels.size() == voxels.size());

    //don't use the built-in /, this 
    // can generate inf values (which is correct)
    // that we want to avoid.
    for(size_t ui=0;ui<voxels.size();ui++)
    {
        if(v.voxels[ui] )
            voxels[ui]/=v.voxels[ui];
        else
        {
            ASSERT(!voxels[ui]);
        } 
    }
}   

    
template<class T>
void Voxels<T>::operator/=(const T &v)
{
    //don't use the built-in /, this 
    // can generate inf values (which is correct)
    // that we want to avoid.
    for(size_t ui=0;ui<voxels.size();ui++)
    {
        if( v > T(0) )
            voxels[ui]/=v;
        else
        {
            voxels[ui]=0;
        } 
    }
}   

template<class T>
bool Voxels<T>::operator==(const Voxels<T> &v) const
{
    for(size_t ui=0;ui<3;ui++)
    {
        
        if(v.binCount[ui] != binCount[ui])
            return false;
    }

    return v.voxels == voxels;
}

#ifdef DEBUG
bool runVoxelTests();
#endif

//===
}

#endif