zxhRegistration.h

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

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

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

#ifndef zxhRegistration_h
#define zxhRegistration_h

#include <string>
#include <iostream>
#include <fstream>
#include <stdio.h>
#include <stdlib.h>
#include "zxhImageData.h"
#include "zxhImageParRec.h"
#include "zxhImageGipl.h"
#include "zxhTransformBase.h"
#include "zxhTransformAffine.h"
#include "zxhTransformFFD.h"
#include "zxhTransformFields.h"
#include "zxhImageModelingLinear.h"
#include "zxhOptimizerGradient.h"
#include "zxhMetricFFDBase.h"
#include "zxhTransformLocalAffines.h"
#include "zxhRegistrationStruct.h"

class zxhOptimizerGradient;

namespace zxh
{

void RegHELP() ;
void RegArg( int argc, char* argv[], zxhRegistrationStruct& StructRegistration,
            int &dimension );
void RegArgSetRoiForDim2Reg( zxhRegistrationStruct& StructRegistration, int &dimension ) ;
void RegNewOptimiser( zxhRegistrationStruct& StructRegistration, zxhOptimizerGradient* &optimizer ) ;
void RegOptimiseArg( int argc, char*argv[], zxhRegistrationStruct& StructRegistration );
int RegMetricArg( int argc, char*argv[], zxhRegistrationStruct& StructRegistration, zxhMetricBase*&pMetric ) ;

//void RegFFDArg( int argc, char* argv[],zxhRegistrationStruct& StructRegistration, zxhImageModelingLinearT<float,float> &ModelNVRef ) ;


void SetDilateLocalRegionForIRBReg( zxhTransformLocalAffines* pTransform, float* afLength, bool Dilation ) ;
void RegAffineArg( int argc, char*argv[],zxhRegistrationStruct& StructRegistration  );

//void RegLocalAffinesArg( int argc, char*argv[], zxhRegistrationStruct& StructRegistration , zxhOptimizerGradient* &optimizer );

void RegLocalAffinesBackwardTransformAllTestImages( zxhRegistrationStruct& StructRegistration, zxhTransformBase* &pInverseTransformTestImages ) ;
void RegSet2DInflucence( int iTransType, int dimension,
                        zxhTransformBase* & pTransform,
                        zxhTransformBase* & pDifferential,
                        zxhTransformBase* & pInfluence );
void SetInfluence2D( zxhTransformAffine * aff );
void SetRegOptimizer(zxhOptimizerGradient* optimizer,zxhRegistrationStruct& StructRegistration);
void RegSetTransformInfo( zxhRegistrationStruct& StructRegistration );
void RegPreSetGlobalAffine( zxhRegistrationStruct& StructRegistration ) ;
void InactivateFFDBorders( zxhTransformFFD* pDifferential );
void RegSetDifferentialInfo( zxhRegistrationStruct& StructRegistration ) ; 
bool RegComputeSimilarity(int argc, char * argv[]);
int reg_option(int argc, char * argv[]);
void SaveRegistrationMiddleResults( int argc, char* argv[], zxhRegistrationStruct&StructRegistration );
void SaveRegistrationResult( zxhOptimizerBase*pOptimiser, zxhRegistrationStruct & StructRegistration,
                                bool bComputeInverse, float fInvJacWeighting ); 
void SaveTransform(zxhTransformBase*pTransform,const char*pFilename);
void SaveString2File(std::string str2save,char*pfilename);
void SaveString2File(std::string str2save,std::string filename);
void BackwardTransformImageByList(zxhImageData*pSave, zxhImageData*pImg,
                                  zxhTransformBase*pList[],int nTransform,int interpolation,
                                  bool bBackground=true, PixelTypeDefault background=ZXH_Background);
void ForwardNearestTransformImageByList(zxhImageModelingBase*pTest,zxhImageModelingBase*pRef,zxhImageData*pSave,
                                        zxhTransformBase*pList[],int nTransform,int empty);
void ForwardBSplineTransformImageByList(zxhImageData*pTest,zxhImageData*pSave,
                                        zxhTransformBase*pList[],int nTransform,int empty);
void FillEmptyVoxelByMean3D(zxhImageData*pImage,int empyvoxel,int nSearchSize, int background=0);
bool ReadTransformFromFile(zxhTransformBase*&pTransform,std::string file);

bool NewTransformFromTypeName(zxhTransformBase*&pTransform,std::string type);
bool NewMetricFromTypeName(zxhMetricBase*&pMetric,std::string type);

zxhTransformBase *  ConcatenateTransformsByExponent( zxhTransformBase*pTrans, int eK, const zxhImageInfo*pImageInfo) ;
zxhTransformBase *  ConcatenateTransforms( zxhTransformBase* pTransList[], int nTrans, int dimension,
                                    const float * spacing, const float * extend,
                                    const float  *roiWorldFrom, const float  *roiWorldTo, const zxhImageInfo*pImgInfo ) ;
zxhTransformBase *  RegriddingLocalAffines( zxhTransformBase* pTjAfter, zxhTransformLocalAffines* pLocalAffinesFirst,
                    zxhImageData*pTestImageTemplate,
                    float *roiWorldFrom=0, float *roiWorldTo=0 ) ;
;
template <class PixelType>
bool SetImageIntensityBetween(zxhImageDataT<PixelType>*pIn,
                              PixelType RangeFrom,PixelType RangeTo,
                              PixelType value, int op, PixelType background=ZXH_Background)
{
    const int* size = pIn->GetImageSize();
    for( int it=0; it<size[3]; ++it )
    for(int iz=0; iz<size[2]; ++iz)
    for(int iy=0; iy<size[1]; ++iy)
    for(int ix=0; ix<size[0]; ++ix)
    {
        PixelType inten = pIn->GetPixelGreyscale( ix,iy,iz,it ) ;
        switch(op)
        {
        case 11:  //va
            if( inten>= RangeFrom && inten <= RangeTo )
                pIn->SetPixelByGreyscale( ix,iy,iz,it, inten+value ) ;
            break;
        case 12:   //vs
            if( inten>= RangeFrom && inten <= RangeTo )
                pIn->SetPixelByGreyscale( ix,iy,iz,it, value ) ;
            break;
        case 13:  //VS
            if( inten>= RangeFrom && inten <= RangeTo )
                pIn->SetPixelByGreyscale( ix,iy,iz,it, value ) ;
            else pIn->SetPixelByGreyscale( ix,iy,iz,it, background ) ;
            break;
        default:break;
        }
    }
    return true;
}

template <class PixelType,class MaskType>
bool MaskTestImageByTrans3D(zxhImageDataT<PixelType>*pIn, zxhImageDataT<MaskType>*pMask,
                          zxhTransformBase*pTrans,float value,int compare, float background,
                          int nRanges, float fRangeFrom[], float fRangeTo[])
{
    int* size   = pIn->GetImageSize();
    zxhImageModelingLinearT<MaskType,float> model;
    model.SetImage( pMask );
    bool inside=false;
    bool samedimension= ImageNvoxelDimensionSame( pIn, pMask );
    PixelType Pixel;
    for(int iz=0; iz<size[2]; ++iz)
    for(int iy=0; iy<size[1]; ++iy)
    for(int ix=0; ix<size[0]; ++ix)
    {
        inside=false;
        if( nRanges== 0 )
            inside = true ;
        else
        {
            Pixel = pIn->GetPixelGreyscale(ix,iy,iz);
            for( int i=0; i<nRanges; ++i )
                if( Pixel>=fRangeFrom[i] && Pixel<=fRangeTo[i] )
                {inside = true;break;}
        }
        if( inside == false )
            continue ;
        float gray ;
        if( pTrans==0 && samedimension )
        {
            gray = pMask->GetPixelGreyscale(ix,iy,iz) ;
        }else
        {
            float fco [ImageDimensionMax]   = { float(ix), float(iy), float(iz), 0} ;
            float fco2[ImageDimensionMax]   = {0,0,0,0} ;
            pIn->AdjustToWorldical(fco);
            if(pTrans!=0)
                pTrans->TransformWorldToWorld(fco,fco2);
            else
            {fco2[0]=fco[0]; fco2[1]=fco[1]; fco2[2]=fco[2]; fco2[3]=fco[3];}

            pMask->WorldToImage(fco2);

            gray = model.GetPixelFloatValueWithCheckBk(background, fco2[0],fco2[1],fco2[2] );
        }

        switch(compare)
        {
        case 0:  if( gray == value ) pIn->SetPixelByGreyscale(ix,iy,iz,0, PixelType(background) ); //=
            break;
        case 1:  if( gray >= value ) pIn->SetPixelByGreyscale(ix,iy,iz,0, PixelType(background) ); //>=
            break;
        case -1: if( gray <= value ) pIn->SetPixelByGreyscale(ix,iy,iz,0, PixelType(background) ); //<=
            break;
        case 2:  if( gray != value ) pIn->SetPixelByGreyscale(ix,iy,iz,0, PixelType(background) ); 
        default:break;
        }

    }
    return true;
}
template <class PixelType,class MaskType>
bool MaskRefImageByTrans3D(zxhImageDataT<PixelType>*pIn, zxhImageDataT<MaskType>*pMask,
                          zxhTransformBase*pTrans,float value,int compare, float background )
{
    zxhImageData RefMask ;
    RefMask.NewImage( pIn->GetDimension(), pIn->GetImageSize(), pIn->GetImageSpacing() ); CopyImageInfo( pIn, &RefMask ) ;
    int * size  = pMask->GetImageSize();
    for(int iz=0;iz<size[2];++iz)
    for(int iy=0;iy<size[1];++iy)
    for(int ix=0;ix<size[0];++ix)
    {
        PixelType gray = pMask->GetPixelGreyscale( ix, iy, iz, 0 ) ;
        bool inside = false;
        switch(compare)
        {
        case 0:if( gray == value ) inside = true ;//=
            break;
        case 1:if( gray >= value ) inside = true ;//>=
            break;
        case -1:if(gray <= value ) inside = true ;//<=
            break;
        default:break;
        }
        if ( inside )
        {
            float fco [ImageDimensionMax]   = {float(ix), float(iy), float(iz), 0} ;
            float fco2[ImageDimensionMax]   = {0,0,0,0} ;
            pMask->AdjustToWorldical(fco);
            if(pTrans!=0)
                pTrans->TransformWorldToWorld(fco,fco2);
            else
                {fco2[0]=fco[0]; fco2[1]=fco[1]; fco2[2]=fco[2]; fco2[3]=fco[3];}
            pIn->WorldToImage( &fco2[0] ) ;
            for ( int iix = int(fco2[0]) ; iix <= int(fco2[0])+1 ; ++iix )
            for ( int iiy = int(fco2[1]) ; iiy <= int(fco2[1])+1 ; ++iiy )
            for ( int iiz = int(fco2[2]) ; iiz <= int(fco2[2])+1 ; ++iiz )
            {
                if ( pIn->InsideImage( iix, iiy, iiz, 0 ) )
                    RefMask.SetPixelByGreyscale( iix, iiy, iiz, 0, 111 );
            }
        }
    }
    size    = pIn->GetImageSize();
    for(int iz=0;iz<size[2];++iz)
    for(int iy=0;iy<size[1];++iy)
    for(int ix=0;ix<size[0];++ix)
    {
        if ( RefMask.GetPixelGreyscale ( ix, iy, iz, 0 ) == 111 )
            pIn->SetPixelByGreyscale ( ix, iy, iz, 0, PixelType(background) );
    }
    return true;
}
bool MaskFFDCtrPntStatus( zxhTransformFFDBase* pInfl, //influence of the optimiser for the FFD control points
                          zxhTransformBase* pTrans, // transform on the pInfl control points, could be null which means the mask is on the test image spacing before transformation
                          zxhImageData * pMask,
                          bool bActivate=false); // also to activate the control points within the mask

float ComputeAverageDisplacement(
                zxhImageData* pTest, int *iRoiTestFrom, int *iRoiTestTo, zxhImageData* pMaskTest,
                zxhTransformBase** pTransformList, int nTransforms,
                zxhImageData* pMaskRef );


float GetMaskImageVolumeByResample(zxhImageData*pImg, float spacex, float spacey, float spacez );
}// end namespace zxh
#endif