mitkPointSet.h

/*=========================================================================

  Program:   3DMed
  Date:      $Date: 2014-02-25 18:30:00 +0800 $
  Version:   $Version: 4.6.0 $
  Copyright: MIPG, Institute of Automation, Chinese Academy of Sciences

=========================================================================*/


#ifndef __mitkPointSet_h
#define __mitkPointSet_h

#include "mitkObject.h"
#include "mitkMesh.h"
#include "mitkGlobalRegistrationFramework.h"


template<class T> class mitkPoint;
template<class T> class mitkBoundingBox;

class MITK_REGISTRATION_API mitkPointSet
{
public:
    mitkPointSet(); 

    virtual ~mitkPointSet();

    mitkPointSet(unsigned int numberOfPoints,unsigned int pointSetDimension,double fillElement=0);

    mitkPointSet(unsigned int num, unsigned int dim, double* p);

    mitkPointSet(mitkPointSet* pointSet); 

    void GetPoint(unsigned int i, ScalarParameterType* p);
    void GetPoint(unsigned int i, ScalarIndexType* p);
    void GetPoint(unsigned int i, ScalarPixelType* p);
    void GetPoint(unsigned int i, mitkPoint<double> & p);
    void GetPoint(unsigned int i, VectorParameterType & p);
    void GetPoint(unsigned int i, VectorIndexType & p);
    void GetPoint(unsigned int i, VectorPixelType & p);

    double GetPoint(unsigned int i, unsigned int dim);

    void SetPoint(unsigned int i, ScalarParameterType* p);
    void SetPoint(unsigned int i, ScalarIndexType* p);
    void SetPoint(unsigned int i, ScalarPixelType* p);
    void SetPoint(unsigned int i, mitkPoint<double> & p);
    void SetPoint(unsigned int i, VectorParameterType & p);
    void SetPoint(unsigned int i, VectorIndexType & p);
    void SetPoint(unsigned int i, VectorPixelType & p);
    
    vector<double>* GetPointsInSequence(vector<int>* sequence);

    vector<double>* GetPointData();

    double* GetPointDataPtr() {return m_PointSet;}  

    void SetPointData(vector<double>* pointSet, unsigned int dimension);

    void SetPointData(double* pointSet, unsigned int dimension, unsigned int numberOfPoints);

    vector<double>* GetPointsInDimension(unsigned int sDimension);

    void SetPointsInDimension(vector<double>* points, unsigned int sDimension);
    
    void SetPointsInDimension(double fillElement, unsigned int sDimension);

    void SetPointData(mitkMesh* mesh);
    
    void GetPointData(mitkMesh* mesh);

    void GetWeight(unsigned int i, vector<double>* w);

    void SetWeight(unsigned int i, vector<double>* w);

    vector<double>* GetWeightsInSequence(vector<int>* sequence);

    void SetWeightData(vector<double>* weights, unsigned int dimension);

    void SetWeightData(double* weights, unsigned int dimension, unsigned int numberOfPoints);

    vector<double>* GetWeightData();

    double* GetWeightDataPtr() {return m_Weights;}

    void SetPointSetMask(vector<int>* pointSetMask);

    vector<int>* GetPointSetMask() {return m_PointSetMask;}

    void SetPointSetMaskFlag(bool maskFlag) {m_FlagUseMask = maskFlag;}

    bool GetPointSetMaskFlag() {return m_FlagUseMask;}

    void ReadPointDataFromFile(char *inFileName,int pointSetDimension,int dataType);
    
    void WritePointDataToFile(char* outFileName);

    unsigned int GetNumberOfPoints();

    void SetNumberOfPoints(unsigned int num) {m_NumberOfPoints = num;}

    int GetPointSetDimension() {return m_PointSetDimension;}

    void SetPointSetDimension(unsigned int pointSetDimension) {m_PointSetDimension = pointSetDimension;}

    int GetWeightsDimension() {return m_WeightsDimension;}

    void GetBoundingBox(double b[6]);

    void GetBoundingBox(float b[6]);

    void GetCentroid(double p[3]);

    virtual void GetNearestPoint(double* query_point, double* result_point, int &i);

    virtual void GetNearestIndex(double* query_point, int &i);

    virtual void GetKNearestPoints(double* query_point, unsigned int k, vector<int> &indices);

    virtual void GetKNearestPoints(vector<double>* query_point, unsigned int k, vector<int> &indices);

    virtual void GetNearestPoints(mitkPointSet* query_pointset, vector<int> &indices); 

    virtual void GetPointsInBoundingBox(mitkBoundingBox<double>& box, vector<int> &indices);

    virtual void GetPointsInRadius(double* query_point, double radius, vector<int> &indices);

    void Fill(double fillElement);

    bool RemoveDuplicates();

    void ConvertPointSet2Dto3D();

    //int GetPointLabel(double x,double y,double z);

    virtual void Initialize();

    virtual void Clear();
    
    bool IsPointSetEmpty();
    
    bool IsWeightsEmpty();

    bool IsInitialized() {return m_FlagInitialized;}

    void SetInitializedFlag(bool f) { m_FlagInitialized = f; }

    void ShallowCopy(mitkPointSet* pointSet);

    void DeepCopy(mitkPointSet* pointSet);

    void Merge(mitkPointSet* pointSet);

    void Reserve(unsigned int s, unsigned int d = 3);

    void SetOriginAndDirection(double o[3], int d[3]);

    void SetSpacing(double s[3]);

    double & operator()(unsigned int i, unsigned int d)
    {
        return m_PointSet[i * m_PointSetDimension + d];
    }

protected:
    void Allocate();
    void AllocateWeights();

    void _init();

    unsigned int                m_PointSetDimension;
    unsigned int                m_WeightsDimension;
    unsigned int                m_NumberOfPoints;

    vector<int>*                m_PointSetMask;
    double*                     m_PointSet;
    double*                     m_Weights;

    double                      m_Centroid[3];
    double                      m_BoundingBox[6];

    double                      m_Spacings[3];
    double                      m_Origins[3];
    int                         m_Directions[3];

    bool                        m_FlagComputeCentroid;
    bool                        m_FlagComputeBoundingBox;
    bool                        m_FlagUseMask;
    bool                        m_FlagInitialized;

private:
    mitkPointSet(const mitkPointSet&);
    void operator=(const mitkPointSet&);
};

//----------------------------------------------------------------------------
template<class T>
class mitkPoint
{
public:
    mitkPoint() :m_dim(0), m_point_data(0) {}

    ~mitkPoint() {m_point_data.clear();}

    mitkPoint(unsigned int dim)
    {
        m_dim = dim;
        m_point_data.resize(m_dim);
    }

    mitkPoint(unsigned int pointDimension, T value)
    {
        m_dim = pointDimension;
        m_point_data.assign(m_dim,value);
    }

/*  mitkPoint(T x, T y)
    {
        m_dim = 2;
        m_point_data.resize(m_dim);
        m_point_data[0] = x;
        m_point_data[1] = y;
    }
*/
    mitkPoint(T x, T y, T z)
    {
        m_dim = 3;
        m_point_data.resize(m_dim);
        m_point_data[0] = x;
        m_point_data[1] = y;
        m_point_data[2] = z;
    }

    mitkPoint(T a, T b, T c, T d)
    {
        m_dim = 4;
        m_point_data.resize(m_dim);
        m_point_data[0] = a; m_point_data[1] = b; m_point_data[2] = c; m_point_data[3] = d;
    }

    mitkPoint(unsigned int pointDimension, T* value)
    {
        m_dim = pointDimension;
        m_point_data.resize(m_dim);
        for(int i=0;i<m_dim;i++) m_point_data[i] = value[i];
    }

    mitkPoint(const vector<T>& value)
    {   
        m_dim = value.size();
        m_point_data.resize(m_dim);
        m_point_data = value;
    }

    mitkPoint(const mitkPoint& p)
    {
        m_dim = p.m_dim;
        m_point_data = p.m_point_data;
    }

    mitkPoint<T>& operator= (const mitkPoint<T>& p)
    {
        m_dim = p.m_dim;
        m_point_data = p.m_point_data;
        return *this;
    }
    
    mitkPoint<T>& operator= (const vector<T>& p)
    {
        if(m_dim == p.size()) m_point_data = p;
        return *this;
    }

    mitkPoint<T>& operator= (const T p[3])
    {
        m_dim = 3;
        m_point_data.resize(3);
        for(int i=0;i<3;i++) m_point_data[i] = p[i];
        return *this;
    }

    mitkPoint<T> operator+ (const mitkPoint<T>& p)
    {
        if(p.m_dim != m_dim)    return *this;

        mitkPoint<T> point(m_dim);
        for(int i=0;i<m_dim;i++)
        {
            point.m_point_data[i] = m_point_data[i] + p.m_point_data[i];
        }
        return point;
    }

    mitkPoint<T> operator- (const mitkPoint<T>& p)
    {
        if(p.m_dim != m_dim)    return *this;

        mitkPoint<T> point(m_dim);
        for(unsigned int i=0;i<m_dim;i++)
        {
            point.m_point_data[i] = m_point_data[i] - p.m_point_data[i];
        }
        return point;
    }

    mitkPoint<T> operator* (const mitkPoint<T>& p)
    {
        if(p.m_dim != m_dim)    return *this;

        mitkPoint<T> point(m_dim);
        for(int i=0;i<m_dim;i++)
        {
            point.m_point_data[i] = m_point_data[i] * p.m_point_data[i];
        }
        return point;
    }

    mitkPoint<T> operator/ (const mitkPoint<T>& p)
    {
        if(p.m_dim != m_dim)    return *this;

        mitkPoint<T> point(m_dim);
        for(unsigned int i=0;i<m_dim;i++)
        {

            point.m_point_data[i] = p.m_point_data[i] == 0 ? 0 : m_point_data[i] / p.m_point_data[i];
        }
        return point;
    }

    bool operator== (const mitkPoint<T>& p)
    {
        if(p.m_dim != m_dim)    return 0;

        bool m = true;
        for(unsigned int i=0;i<m_dim;i++) m &= m_point_data[i] == p.m_point_data[i];
        return m;
    }

    bool operator!= (const mitkPoint<T>& p)
    {
        if(p.m_dim != m_dim)    return 0;

        bool m = true;
        for(int i=0;i<m_dim;i++) m &= m_point_data[i] != p.m_point_data[i];
        return m;
    }

    T operator[] (unsigned int i) const
    {
        return m_point_data[i];
    }

    unsigned int size()
    {
        return m_dim;
    }

    void resize(unsigned int s) 
    {
        m_point_data.resize(s);
        m_dim = s;
    }

    void set(unsigned int i, const T & value)
    {
        if(i>=0 && i<m_dim)
            m_point_data[i] = value;
    }

    T get_norm() const
    {
        T sum = 0;
        for(unsigned int i=0;i<m_dim;i++)
        {
            const T v = m_point_data[i];
            sum += v * v;
        }
        return sqrt(double(sum));
    }

    vector<T>       m_point_data;
    unsigned int    m_dim;
};
//----------------------------------------------------------------------------
template<class T>
class mitkBoundingBox
{
public:
    mitkBoundingBox() :m_dim(3),m_bound_data(0) {}

    ~mitkBoundingBox() {}

    //construct with elements directly
    mitkBoundingBox(T x0, T x1, T y0, T y1)
    {
        m_dim = 2;
        m_bound_data.resize(m_dim*2);
        m_bound_data[0] = x0;
        m_bound_data[1] = x1;
        m_bound_data[2] = y0;
        m_bound_data[3] = y1;
    }

    //construct with elements directly
    mitkBoundingBox(T x0, T x1, T y0, T y1, T z0, T z1)
    {
        m_dim = 3;
        m_bound_data.resize(m_dim*2);
        m_bound_data[0] = x0;
        m_bound_data[1] = x1;
        m_bound_data[2] = y0;
        m_bound_data[3] = y1;
        m_bound_data[4] = z0;
        m_bound_data[5] = z1;
    }
    
    //construct with arrays
    mitkBoundingBox(T* min, T* max, unsigned int dim)
    {
        m_dim = dim;
        m_bound_data.resize(m_dim*2);
        for(unsigned int i=0;i<m_dim;i++)
        {
            m_bound_data[2*i]   = min[i];
            m_bound_data[2*i+1] = max[i];
        }
    }
    
    //construct with vectors
    mitkBoundingBox(vector<T>* min, vector<T>* max)
    {
        m_dim = min->size();
        m_bound_data.resize(m_dim*2);
        for(int i=0;i<m_dim;i++)
        {
            m_bound_data[2*i]   = (*min)[i];
            m_bound_data[2*i+1] = (*max)[i];
        }
    }
    
    //construct with scalar values
    mitkBoundingBox(T min, T max, unsigned int dim)
    {
        m_dim = dim;
        m_bound_data.resize(m_dim*2);
        for(unsigned int i=0;i<m_dim;i++)
        {
            m_bound_data[2*i]   = min;
            m_bound_data[2*i+1] = max;
        }
    }
    
    //construct with mitkBoundingBox
    mitkBoundingBox(const mitkBoundingBox& b)
    {
        m_dim = b.m_dim;
        m_bound_data.resize(m_dim*2);
        m_bound_data = b.m_bound_data;
    }
    
    //construct with an array contains min & max
    mitkBoundingBox(T* b, unsigned int dim)
    {
        m_dim = dim;
        m_bound_data.resize(dim*2);
        for(int i=0;i<m_bound_data.size();i++) m_bound_data[i] = b[i];
    }

    mitkBoundingBox& operator= (const mitkBoundingBox& b)
    {
        m_dim = b.m_dim;
        m_bound_data = b.m_bound_data;
        return *this;
    }

    //T& operator[] (unsigned int i)
    //{
    //  return m_bound_data[i];
    //}

    void GetHighBound(T* bmax) const
    {
        for(unsigned int i=0;i<m_dim;i++) bmax[i] = m_bound_data[2*i+1];
    }

    void GetLowBound(T* bmin) const
    {
        for(unsigned int i=0;i<m_dim;i++) bmin[i] = m_bound_data[2*i];
    }

    void SetHighBound(T* max) 
    {
        for(int i=0;i<m_dim;i++) m_bound_data[i*2+1] = max[i];
    }

    void SetLowBound(T* min)
    {
        for(int i=0;i<m_dim;i++) m_bound_data[i*2] = min[i];
    }
    unsigned int GetDimension() const {return m_dim;}

    vector<T>       m_bound_data;
    unsigned int    m_dim;
};
//----------------------------------------------------------------------------
template<class T1, class T2>
inline T2 mitk_point_dist(T1* p1, T2* p2, unsigned int dim) //distance - point to point
{
    T2 sum_sq = 0;
    for(unsigned int i=0; i < dim; i++ ) 
    {
        T2 tmp = (*p1)[i] - (*p2)[i];
        sum_sq += tmp * tmp;
    }
    sum_sq = sqrt(sum_sq);
    return sum_sq;
}
//----------------------------------------------------------------------------
template<class T>
inline T mitk_point_dist_sq(T* p1, T* p2, unsigned int dim) //distance^2 - point to point
{
    T sum_sq = 0;
    for(unsigned int i=0; i < dim; i++ ) 
    {
        T tmp = p1[i] - p2[i];
        sum_sq += tmp * tmp;
    }
    return sum_sq;
}
//----------------------------------------------------------------------------
template<class T>
inline T mitk_point_dist_sq(vector<T>* p1, T* p2, unsigned int dim) //distance^2 - point to point
{
    T sum_sq = 0;
    for(unsigned int i=0; i < dim; i++ ) 
    {
        T tmp = (*p1)[i] - p2[i];
        sum_sq += tmp * tmp;
    }
    return sum_sq;
}
//----------------------------------------------------------------------------
template<class T>
inline T mitk_point_dist_sq(vector<T>* p, const mitkBoundingBox<T>& b, bool& inside)  //distance^2 - point to boundingbox
{   
    T sum_sq = 0;
    unsigned int dim = b.GetDimension();
    T x0[3], x1[3];

    b.GetLowBound(x0);
    b.GetHighBound(x1);

    inside = true;

    for( unsigned int i=0; i<dim; i++ )
    {
        if( (*p)[i] < x0[i] ) 
        {
            T tmp = x0[i] - (*p)[i];
            sum_sq += tmp * tmp;
            inside = false;
        }
        else if( (*p)[i] > x1[i] ) 
        {
            T tmp = x1[i] - (*p)[i];
            sum_sq += tmp * tmp;
            inside = false;
        }  
    }

    return sum_sq;
}
//----------------------------------------------------------------------------
template<class T>
inline T mitk_point_dist_sq(T* p, const mitkBoundingBox<T>& b, bool& inside)  //distance^2 - point to boundingbox
{   
    T sum_sq = 0;
    unsigned int dim = b.GetDimension();
    T x0[3], x1[3];

    b.GetLowBound(x0);
    b.GetHighBound(x1);

    inside = true;

    for( unsigned int i=0; i<dim; i++ )
    {
        if( p[i] < x0[i] ) 
        {
            T tmp = x0[i] - p[i];
            sum_sq += tmp * tmp;
            inside = false;
        }
        else if( p[i] > x1[i] ) 
        {
            T tmp = x1[i] - p[i];
            sum_sq += tmp * tmp;
            inside = false;
        }  
    }

    return sum_sq;
}
//----------------------------------------------------------------------------


//#define DEFINED_mitkPointSet

#endif

---