zxhMetricMI.h

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

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

=========================================================================*/
#ifndef zxhMetricMI_h
#define zxhMetricMI_h

#include "zxh.h"
#include "zxhMetricBase.h"
//#include "zxhTransformFFDBase.h"

#ifndef ZXHMMIComputeKernelFunctionInLocalRegion
#define ZXHMMIComputeKernelFunctionInLocalRegion 1.0  // 0.5
#endif


class zxhMetricMI :
    public zxhMetricBase
{
public:
    zxhMetricMI(void);
    ~zxhMetricMI(void);

    virtual std::string GetMetricType() {if(GetIsMMI()) return "_SEMI_MI_" ; else return "_MI_"; } ;

    const static int m_stc_iPaddedBins = BSplineHistogramPaddingBins;

    virtual zxhMetricBase* Clone(zxhMetricBase*&pRet); 

    virtual bool PreComputeImageGraylevel();

    virtual bool SetImage(zxhImageData*pImageTest,zxhImageData*pImageRef);
 
    virtual int ComputeHistogramFromPDF4D(zxhlfloat * &afPDF,zxhlfloat&fTotalFreq,
                                            float xfTestFrom,float xfTestTo,float yfTestFrom,float yfTestTo,
                                            float zfTestFrom,float zfTestTo,float tfTestFrom,float tfTestTo);
    virtual int     ComputeHistogramMaskRegionOnTest( zxhImageData*pMask, zxhlfloat * &afPDF, zxhlfloat &fFreq );
    virtual int     ComputeHistogram(zxhlfloat*&afHistogram, zxhlfloat&fTotalFreq);
    virtual int     GetValue(zxhlfloat& fValue);
    virtual zxhlfloat  GetMIValue( const zxhlfloat*afHistogram, const zxhlfloat fTotalFreq );
    virtual void    AddHistogram(zxhlfloat* afHistogramDes, const zxhlfloat* afHistogramSrc, int times=1 );
    virtual void    GetAddHistogram(zxhlfloat* afHistogramDes, zxhlfloat* afHistogramSrc1, zxhlfloat* afHistogramSrc2 );
    
    void    SetNumberOfHistogramBins(int i=32)          {m_iNumberOfHistogramBins=i;};

    int     GetNumberOfHistogramBins()                  {return m_iNumberOfHistogramBins;};

    void    SetPDFType( enum zxh::TypeHistogramCompute i )              {m_iPDFType=i;};
    void    SetMMILambda(float f)                                       {m_fMMILambda=f;};
    float   GetMMILambda()                                              {return m_fMMILambda;};
    void    SetMMISamplePhysical(float*f)
    {
        for( int id=0; id<m_iDimension; ++id )
            m_afMMISamplePhysical[id] = f[id];
    }
    const float * GetMMISamplePhysical()                                {return m_afMMISamplePhysical;};
    void    SetMMILocalRegionSizeInMM(float f)                          { m_fMMILocalRegionSizeInMM=f; };
    void    SetMMIKernelType( enum zxh::TypeSEMIKernelFunction type )   { m_iMMIKernelType = type;}
    void    SetMMI( bool b )    
    {   if(b) m_iComputePDType=zxh::MixtureProbabilityModel;
        else m_iComputePDType=zxh::NormalizedMutualInformation;};
    bool    GetIsMMI()
    {   if(m_iComputePDType==zxh::MixtureProbabilityModel) return true; 
        else return false;};
    void    SetMIType( enum zxh::TypeMutualInformation i )              {m_iMIType=i;};
    enum zxh::TypeMutualInformation GetMIType()                         {return m_iMIType;};

    virtual std::string GetPrintString();

    virtual bool        SetMetricFromStream(std::ifstream & ifs);
 
    virtual int ComputeConstanceForGradient() ;  

protected:

        inline float SEMIKernelDistanceToWeight(float dis )
        {
            if( dis*2>=m_fMMILocalRegionSizeInMM ) return 0;
            switch( m_iMMIKernelType )
            {
            case zxh::SEMIGaussian: return zxh::NGaussian(dis*6.0/m_fMMILocalRegionSizeInMM); break;
            case zxh::SEMIBSpline: return zxh::BSpline(dis*4/m_fMMILocalRegionSizeInMM); break;
            case zxh::SEMIBSplineZeroOrder: return 1; break;
            default: std::cerr<<"error: wrong entry calling semi distance to weight\n"; break;
            }
            return 0 ;
        }

    enum zxh::TypeMutualInformation m_iComputePDType;
    /* 0~1 : mixture pdf (E(x)=1)
     * -1~0: pixel-wise weighted pdf (E(x)=Gaussian)
     * 2   : 0.5LMI+0.5GMI
     */
    float   m_fMMILambda ; 
    float   m_afMMISamplePhysical[ImageDimensionMax] ;

    float   m_fMMILocalRegionSizeInMM ;
    enum zxh::TypeSEMIKernelFunction m_iMMIKernelType ;
    zxhImageData m_imgMMItestImage;
    zxhImageData m_imgMMIrefImage;
    zxhImageDataT<float> m_imgMMIWeighttestref ;

    enum zxh::TypeHistogramCompute m_iPDFType;
    
        enum zxh::TypeMutualInformation m_iMIType;

protected:
    int m_iNumberOfHistogramBins;
    PixelTypeDefault m_aiImageTestGraySpan[2]; 
    PixelTypeDefault m_aiImageRefGraySpan[2];
    zxhlfloat   m_fBinSizeTestImage ;
    zxhlfloat   m_fBinSizeRefImage ;

    zxhlfloat   m_fEntropyT; 
    zxhlfloat   m_fEntropyR; 
    zxhlfloat   m_fEntropyJ;
    zxhlfloat   m_fMutualInformation ;
    zxhlfloat   m_fNormalisedMutualInformation ;
    zxhlfloat   *m_afPDFt;
    zxhlfloat   *m_afPDFr;
    zxhlfloat   *m_afPDFj;

    zxhlfloat   *m_afHistogram; 
    zxhlfloat   m_fTotalFreq;
    zxhlfloat   *m_afC_Ht;
    zxhlfloat   *m_afC_Hr;
    zxhlfloat   *m_afC_Hj;

    inline int      GetLengthOfHistogramDim()
    {
        if(m_iPDFType==zxh::PDFType_ParzenWindow_BSpline)
            return (m_iNumberOfHistogramBins+m_stc_iPaddedBins);
        return m_iNumberOfHistogramBins;
    } 
    inline void     IntensityBoundBk(float & gray, const short aiImageGraySpan[] ) const
    {
        gray = MAX(gray, aiImageGraySpan[0]);
        gray = gray>(aiImageGraySpan[1]-1)?aiImageGraySpan[0]:gray;
    }
    inline void     IntensityBoundBk(zxhlfloat &gray, const short aiImageGraySpan[] ) const
    { 
        gray = MAX(gray, aiImageGraySpan[0]) ;
        gray = gray>(aiImageGraySpan[1]-1)?aiImageGraySpan[0]:gray;
    }
    inline void     IntensityBound(zxhlfloat &gray, const short aiImageGraySpan[] ) const
    {
        gray = MAX(gray, aiImageGraySpan[0]) ;
        gray = MIN(gray,(aiImageGraySpan[1]-1));
    }
    void        IntensityBound(float &gray, const short aiImageGraySpan[] ) const
    {
        gray = MAX(gray, aiImageGraySpan[0]) ;
        gray = MIN(gray,(aiImageGraySpan[1]-1));
    }
    virtual void EvaluateOnJointHistogram(zxhlfloat &fTotalFreq,zxhlfloat *afPDF, int iLengthOfHistogramDim,
                                zxhlfloat fGrayTest, zxhlfloat fGrayRef , zxhlfloat weight=1);
    

    //virtual void  EvaluateOnSumingEntropyDerivative( zxhlfloat &dEntropyRef, zxhlfloat &dEntropyJoint,
    //              float fGrayTest, float fGrayRef, float fGrayRefForward, float fGrayRefBackward, float length2) ;
    //virtual void  EvaluateOnPixelForce(   zxhlfloat *afForce, int iLengthOfHistogramDim,
    //                              float * gridCoTest, float * gridCoRef );
    virtual void    EvaluateOnPixelSource2TargetForce(  zxhlfloat *afForce, int iLengthOfHistogramDim,
                                                const float * gridCoTest, const float * gridCoRef ) const ;
            
    inline float    BSplinei(int iOrd, float u)
    {   
        zxhlfloat v;
        switch(iOrd)
        {
        case -1:
            v=1.0-u;
            return v*v*v/6.0;
        case 0:
            v=u*u;
            return 0.5*u*v-v+2.0/3.0;
        case 1:
            v=u*u;
            return (-3.0*u*v+3.0*v+3.0f*u+1)/6.0;
        case 2:return u*u*u/6.0;
        }
        return -1;
    } 
    inline float    DerivativeOfBSplinei(int iOrd,float u)
    { // -1,0,1,2
        zxhlfloat v;
        switch(iOrd)
        {
        case -1:
            v=1.0-u;
            return v*v*(-0.5);
        case 0:return (3.0*u-4.0)*u*0.5;
        case 1:return (u+0.5-1.5*u*u);
        case 2:return u*u*0.5;
        }
        return -1;
    };
    /* if m_Lambda < 0, then regards using lMI +Lambda*gMI
     * **/
    zxhlfloat   ComputMMILocalRegionPDFValue( zxhlfloat*afHistL, zxhlfloat&fVl, 
        const zxhlfloat *afHistG, const zxhlfloat fVg );  
    void EvaluateMMIpositionWeightedHistogram(int *roiFrom, int*roiTo, float* WorldCentre, zxhlfloat *afHistogram, zxhlfloat &fTotalFreq,
                                              int *excludeFrom, int*excludeTo );
    void EvaluateMMIpositionWeightedHistogram(int *roiFrom, int*roiTo, float* WorldCentre, zxhlfloat *afHistogram, zxhlfloat &fTotalFreq,
                                              zxhImageData *pExcludeImageMask );
};

#endif