zxhTransformFFDBase.h

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

  Program:   ZXH Registration Software
  Author:    Xiahai Zhuang
  Module:    $RCSfle: zxhTransformFFDBase.h    $
  Language:  C++
  Date:      $Date: From  2004-01 $
  Version:   $Revision: 1.0, 2.0, 2.1 $

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

#ifndef ZXHTRANSFORMFFDBASE_H
#define ZXHTRANSFORMFFDBASE_H

#include "zxhTransformBase.h"
#include "zxhImageInfo.h"


class zxhTransformFFDBase :
    public zxhTransformBase
{
public:
    zxhTransformFFDBase(void);
    virtual ~zxhTransformFFDBase(void);
    virtual std::string GetTransformType(){return "FFD";};

    virtual bool    SetDimension(int i)                 {m_ImageInfo.Dimension=i; return zxhTransformBase::SetDimension(i);};

    void    MallocForControlPoints(int ni,float fi,int nj=1,float fj=0,  int nk=1,float fk=0,int nl=1,float fl=0);
    void    MallocForControlPoints(float fi, float fj=1, float fk=1, float fl=1);
    void    MallocForControlPoints(int ni,int nj=1,int nk=1,int nl=1);
    void    FreeControlPoints();

    bool    ConvertBetweenImage( zxhImageDataT<float> & ffdimg, bool fromFFDtoImage ) ;

    inline float*   GetCtrPntsWithCheck(int i,int j=0,int k=0,int l=0) const
    {
        if(i<0||j<0||k<0||l<0||
            i>=m_ImageInfo.Size[0]||
            j>=m_ImageInfo.Size[1]||
            k>=m_ImageInfo.Size[2]||
            l>=m_ImageInfo.Size[3])
        {
            return 0;//m_fZeroCtrPnt;
        }
        return (m_pfCtrPntOffset+(l*m_nCtrPntVolume+k*m_nCtrPntResolution+j*m_ImageInfo.Size[0]+i)*m_iDimension);
    };
    inline const float* GetCtrPntsDisplacement(int i,int j=0,int k=0,int l=0)
    {
        if(i<0||j<0||k<0||l<0||
            i>=m_ImageInfo.Size[0]||
            j>=m_ImageInfo.Size[1]||
            k>=m_ImageInfo.Size[2]||
            l>=m_ImageInfo.Size[3])
        {
            return m_fZeroCtrPnt;
        }
        return (m_pfCtrPntOffset+(l*m_nCtrPntVolume+k*m_nCtrPntResolution+j*m_ImageInfo.Size[0]+i)*m_iDimension);
    };
    bool    GetCtrPntsDisplacement(float pdis[], int i,int j=0,int k=0,int l=0)
    {
        bool ret=false ; float *p = &m_fZeroCtrPnt[0] ;
        if(i>=0&&j>=0&&k>=0&&l>=0&& i<m_ImageInfo.Size[0]&& j<m_ImageInfo.Size[1]&& k<m_ImageInfo.Size[2]&& l<m_ImageInfo.Size[3])
        {   ret=true ; p=(m_pfCtrPntOffset+(l*m_nCtrPntVolume+k*m_nCtrPntResolution+j*m_ImageInfo.Size[0]+i)*m_iDimension); }
        for( int i=0;i<m_iDimension;++i) pdis[i] = p[i] ;
        return ret ;
    };
    inline float*   GetCtrPntsWithoutCheck(int i,int j=0,int k=0,int l=0)
    {
        return (m_pfCtrPntOffset+(l*m_nCtrPntVolume+k*m_nCtrPntResolution+j*m_ImageInfo.Size[0]+i)*m_iDimension);
    };
    inline void SetOffsetOfCtrPnts(float pOffset[ImageDimensionMax],int i,int j=0,int k=0,int l=0)
    {
        for(int index=0;index<m_iDimension;++index)
            (m_pfCtrPntOffset+(l*m_nCtrPntVolume+k*m_nCtrPntResolution+j*m_ImageInfo.Size[0]+i)*m_iDimension)[index]=pOffset[index];
    };

    void GetCtrGridSize(int& ix, int&iy, int&iz, int&it )           { ix=m_ImageInfo.Size[0]; iy=m_ImageInfo.Size[1]; iz=m_ImageInfo.Size[2]; it=m_ImageInfo.Size[3]; } ;
    void GetCtrGridSpacing(float& fx, float&fy, float&fz, float&ft )
    { fx=m_ImageInfo.Spacing[0]; fy=m_ImageInfo.Spacing[1]; fz=m_ImageInfo.Spacing[2]; ft=m_ImageInfo.Spacing[3]; } ;
    const float*    GetCtrGridSpacing(void)         {return &m_ImageInfo.Spacing[0];};

    const int * GetCtrGridSize(void)            {return &m_ImageInfo.Size[0];} ;
    //const float*  GetCtrGridSpacing(void)         {return &m_ImageInfo.Spacing[0];};
    void GetCtrGridExtent(float e[4])
    {   for( int id=0; id<m_iDimension; ++id ) e[id] = m_ImageInfo.Spacing[id]*float(m_ImageInfo.Size[id]-1) ;} ;
    void WorldToGrid( float fv[] ) const
    {   m_ImageInfo.WorldToImage( fv ) ; } ;
    void GridToWorld( float fv[] ) const
    {   m_ImageInfo.ImageToWorld( fv ) ; } ;
    void GridToPhysical( float fv[] ) const
    {   m_ImageInfo.ImageToPhysical(fv); }
    void PhysicalToGrid( float fv[] ) const
    {   m_ImageInfo.PhysicalToImage(fv); }

    virtual void    GetOrientationAsImageInfo( zxhImageInfo*p )
    {if(p==0){ std::cerr<<"error: get orientation as image info\n"; return ;} p->CopyOrientationInfoFrom( & m_ImageInfo ) ; };
    virtual const zxhImageInfo*     GetOrientationImageInfo()   {return &m_ImageInfo ; } ;
    virtual void    SetOrientationUsingImageInfo(const zxhImageInfo*p ) ;

    virtual std::string GetPrintString() ;

    virtual bool    SetTransformFromFile(const char* pFile);
    virtual bool    SetTransformFromStream(std::ifstream& ifs);
    virtual bool    SetTransformIdentity();

    virtual bool    SetImage(zxhImageData*pTest,zxhImageData*pRef) ;

    virtual zxhTransformBase*   Clone(zxhTransformBase*&pRet);
    virtual void    TransformPhysToPhys(const float fVectorFromPhys[ImageDimensionMax],float fVectorToPhys[ImageDimensionMax]);
    virtual void    TransformWorldToWorld(const float fVectorFromWorld[ImageDimensionMax],float fVectorToWorld[ImageDimensionMax]);

    virtual bool    GetJacobianMatrixOnGrid(float fGrid[ImageDimensionMax],float pJacobian[ImageDimensionMax*ImageDimensionMax]);

    ;
    virtual bool    ResampleControlGrid(int ix,int iy,int iz,int it,  int level=1, float ERROR_E=0.001);
    virtual bool    ResampleControlGrid(float sx, float sy, float sz, float st=1,  int level=1, float ERROR_E=0.001);
    virtual bool    ResampleControlGrid( int* size, float* spacing,  int level=1, float ERROR_E=0.001) ;

    virtual float   ApproximateFFD(zxhTransformBase*trans, int level=1, float ERROR_E=0.001) ;
    virtual float   ApproximateFFD(zxhTransformBase*trans, zxhImageInfo &imageinfo, int level=1, float ERROR_E=0.001) ;

    virtual bool    ResampleDownHalf() ;

    virtual bool    SetTransformPara(float) ;
    virtual bool    GetAdd(zxhTransformBase*,zxhTransformBase*) ;
    virtual bool    Add(zxhTransformBase*) ;
    virtual bool    GetSubtract(zxhTransformBase*,zxhTransformBase*) ;
    virtual bool    Subtract(zxhTransformBase*) ;
    virtual bool    GetMultiply(zxhTransformBase*,float) ;
    virtual bool    Multiply(float) ;
    virtual bool    GetMultiplyByPara(zxhTransformBase*,zxhTransformBase*) ;
    virtual bool    MultiplyByPara(zxhTransformBase*) ;

    virtual float   GetMagnitudeAsVector(void);
    virtual float   GetMaxLocalMagnitudeAndNormalisation( bool bNormalisationByLocalMagnitude ) ;
    virtual float   GetMaxAbsValueFromParameters(void) ;

    virtual float   Guarantee1To1(void);

    virtual int GetDegreeOfFreedom () ;
    virtual int GetDOFWithValue( float f )   ;
    virtual int     GetNoParameters()               { return GetDegreeOfFreedom(); } ;
    virtual float   GetParameters(int index)
    { if(index>=0&&index<GetDegreeOfFreedom()) return m_pfCtrPntOffset[index]; return 0 ; };
    virtual bool    SetParameterValue( int index, float f )
    {
        if( index<0 || index>=GetDegreeOfFreedom( ) ) return false ;
        m_pfCtrPntOffset[index]=f; return true ;
    };

protected:
    const static unsigned int static_buffersize_aline=16384;

    zxhImageInfo    m_ImageInfo ;

    float       *m_pfCtrPntOffset;

    int         m_nCtrPntVolume;
    int         m_nCtrPntResolution;

    zxhImageDataT<float>    m_imgPreComputeFFDWeight ;
    float       m_afPortionPreSpacingImageSpacing[ImageDimensionMax] ;
    float       m_afPreComputeTransformSpacing[ImageDimensionMax] ;
    float       m_fZeroCtrPnt[4] ;

    inline float    BSplinei(int iOrd,float u) const
    {
        /*if(m_pBSpline)
            return *(m_pBSpline+(iOrd+1)*ZXH_BSPLINE_PLUS+int(u*ZXH_BSPLINE));*/
        float v;
        switch(iOrd)
        {
        case -1:
            v=1-u;
            return v*v*v/6.0f;
        case 0:
            v=u*u*3.0f;
            return (v*u-v*2.0f+4.0f)/6.0f;
        case 1:
            v=3.0f*u*u;
            return (-v*u+v+3.0f*u+1.0f)/6.0f;
        case 2:return u*u*u/6.0f;
        }
        return 0;
    };
    inline float    DerivativeOfBSplinei(int iOrd,float u) const
    { // -1,0,1,2
        float v;
        switch(iOrd)
        {
        case -1:
            v=1.0f-u;
            return v*v*(-0.5f);
        case 0:return (3.0f*u-4.0f)*u*0.5f;
        case 1:return (u+0.5f-1.5f*u*u);
        case 2:return u*u*0.5f;
        }
        return 0;
    };
    inline float    SecondDerivativeOfBSplinei(int iOrd,float u)
    { // -1,0,1,2 
        switch(iOrd)
        {
        case -1: 
            return 1.0f-u;
        case 0:
            return (3.0f*u-2.0f);
        case 1:return (1.0f-3.0f*u);
        case 2:return u;
        }
        return 0;
    };

};


#endif