animaGcStremMsLesionsSegmentationFilter.h

Go to the documentation of this file.

#pragma once

#include "animaTLinksFilter.h"
#include "animaGraphCutFilter.h"
#include "animaCheckStructureNeighborFilter.h"
#include "animaRemoveTouchingBorderFilter.h"
#include "animaComputeMahalanobisImagesFilter.h"
#include "animaComputeSolution.h"

#include <itkMaximumImageFilter.h>
#include <itkMinimumImageFilter.h>
#include <itkRelabelComponentImageFilter.h>
#include <itkRescaleIntensityImageFilter.h>
#include <itkImageToImageFilter.h>
#include <itkBinaryThresholdImageFilter.h>
#include <itkIntensityWindowingImageFilter.h>
#include <itkMaskImageFilter.h>

enum LesionSegmentationType
{
    strem = 0,
    gcem,
    gcemAndManualGC,
    manualGC
};

enum InitializationType
{
    Atlas = 0,
    HierarchicalDP,
    HierarchicalFLAIR
};

namespace anima
{

template <typename TInputImage>
class GcStremMsLesionsSegmentationFilter :
        public itk::ImageToImageFilter< TInputImage, itk::Image <unsigned char, 3> >
{
public:
    typedef GcStremMsLesionsSegmentationFilter Self;
    typedef itk::Image <unsigned char, 3> TOutputImage;
    typedef itk::ImageToImageFilter< TInputImage, TOutputImage > Superclass;
    typedef itk::SmartPointer<Self> Pointer;
    typedef itk::SmartPointer<const Self>  ConstPointer;

    itkNewMacro(Self)

    
    itkTypeMacro(GcStremMsLesionsSegmentationFilter, ImageToImageFilter)

    
    typedef TInputImage InputImageType;
    typedef itk::ImageRegionConstIterator <InputImageType> ImageIteratorTypeConstInput;

    typedef TOutputImage OutputImageType;
    typedef typename OutputImageType::Pointer OutputImagePointer;
    typedef typename itk::ImageRegionIterator< OutputImageType > OutputIteratorType;
    typedef typename OutputImageType::PixelType OutputPixelType;

    typedef unsigned char   PixelTypeUC;
    typedef itk::Image <PixelTypeUC,3> ImageTypeUC;
    typedef itk::ImageRegionIterator <ImageTypeUC> ImageIteratorTypeUC;
    typedef itk::ImageRegionConstIterator <ImageTypeUC> ImageIteratorTypeConstUC;

    typedef double  PixelTypeD;
    typedef itk::Image <PixelTypeD,3> ImageTypeD;
    typedef itk::ImageRegionIterator <ImageTypeD> ImageIteratorTypeD;

    typedef int     PixelTypeInt;
    typedef itk::Image <PixelTypeInt,3> ImageTypeInt;
    typedef itk::ImageRegionIterator <ImageTypeInt> ImageIteratorTypeInt;

    typedef itk::VariableLengthVector<double> MeasurementVectorType;
    typedef itk::Statistics::GaussianMembershipFunction< MeasurementVectorType > GaussianFunctionType;

    typedef anima::GraphCutFilter< TInputImage, ImageTypeUC > GraphCutFilterType;
    typedef anima::CheckStructureNeighborFilter< ImageTypeUC,ImageTypeUC > CheckStructureNeighborFilterFilterType;
    typedef anima::RemoveTouchingBorderFilter<ImageTypeUC,ImageTypeUC,ImageTypeUC> RemoveTouchingBorderFilterType;
    typedef anima::TLinksFilter<TInputImage,ImageTypeD> TLinksFilterType;
    typedef anima::ComputeMahalanobisImagesFilter<ImageTypeUC,ImageTypeUC,ImageTypeD> ComputeMahalanobisImagesFilterType;
    typedef anima::ComputeSolution<ImageTypeUC,ImageTypeUC, ImageTypeD> ComputeSolutionType;

    typedef itk::RescaleIntensityImageFilter<TInputImage,ImageTypeUC> RescaleFilterType;
    typedef itk::MaximumImageFilter<ImageTypeD,ImageTypeD,ImageTypeD> MaximumFilterType;
    typedef itk::MinimumImageFilter<ImageTypeD,ImageTypeD,ImageTypeD> MinimumFilterTypeF;
    typedef itk::MinimumImageFilter<ImageTypeUC,ImageTypeUC,ImageTypeUC> MinimumFilterTypeUC;
    typedef itk::BinaryThresholdImageFilter <ImageTypeD, ImageTypeUC> BinaryThresholdImageFilterType_F_UC;
    typedef itk::BinaryThresholdImageFilter <ImageTypeUC, ImageTypeUC> BinaryThresholdImageFilterType_UC_UC;
    typedef itk::IntensityWindowingImageFilter <ImageTypeUC, ImageTypeD> IntensityWindowingImageFilterType;
    typedef itk::MaskImageFilter< ImageTypeUC, ImageTypeUC > MaskFilterType_UC_UC;
    typedef itk::MaskImageFilter< ImageTypeD, ImageTypeUC > MaskFilterType_F_UC;
    typedef itk::ConnectedComponentImageFilter <ImageTypeUC,ImageTypeInt> ConnectedComponentType;
    typedef itk::RelabelComponentImageFilter< ImageTypeInt, ImageTypeInt > RelabelComponentType;

    void GenerateData() ITK_OVERRIDE;
    void WriteOutputs();

    void CheckInputImages();
    void RescaleImages();
    void ComputeAutomaticInitialization();
    void StremThreshold();
    void CreateFuzzyRuleImage();
    void GraphCut();
    void ApplyHeuristicRules();
    void ComputeNABT();

    void SetInputImageT1(const InputImageType* image);
    void SetInputImageT2(const InputImageType* image);
    void SetInputImageDP(const InputImageType* image);
    void SetInputImageFLAIR(const InputImageType* image);
    void SetInputImageT1Gd(const InputImageType* image);

    void SetInputCSFAtlas(const ImageTypeD* image);
    void SetInputGMAtlas(const ImageTypeD* image);
    void SetInputWMAtlas(const ImageTypeD* image);
    void SetInputLesionPrior(const ImageTypeD* image);

    void SetMask(const TInputImage* MaskImage);

    void SetSourcesMask(const ImageTypeUC* MaskImage);
    void SetSinksMask(const ImageTypeUC* MaskImage);

    TOutputImage* GetOutputWholeSeg();
    TOutputImage* GetOutputLesions();
    TOutputImage* GetOutputCSF();
    TOutputImage* GetOutputGM();
    TOutputImage* GetOutputWM();
    TOutputImage* GetOutputStrem();
    TOutputImage* GetOutputStremCSF();
    TOutputImage* GetOutputStremGM();
    TOutputImage* GetOutputStremWM();
    TOutputImage* GetOutputIntensityImage1();
    TOutputImage* GetOutputIntensityImage2();
    itk::Image <double,3>* GetOutputFuzzyObjectImage();
    itk::Image <double,3>* GetOutputMahaCSFImage();
    itk::Image <double,3>* GetOutputMahaGMImage();
    itk::Image <double,3>* GetOutputMahaWMImage();
    itk::Image <double,3>* GetOutputMahaMinimumImage();
    itk::Image <double,3>* GetOutputMahaMaximumImage();
    TOutputImage* GetOutputGraphCut();

    void SetTol(const double tol)
    {
        this->SetCoordinateTolerance(tol);
        this->SetDirectionTolerance(tol);
        m_Tol = tol;
    }

    itkSetMacro(InitMethodType, unsigned int)
    itkGetMacro(InitMethodType, unsigned int)

    itkSetMacro(RejRatioHierar, double)
    itkGetMacro(RejRatioHierar, double)

    itkSetMacro(MinDistance, double)
    itkGetMacro(MinDistance, double)

    itkSetMacro(EmIter, int)
    itkGetMacro(EmIter, int)

    itkSetMacro(RejRatio, double)
    itkGetMacro(RejRatio, double)

    itkSetMacro(EmIter_concentration, int)
    itkGetMacro(EmIter_concentration, int)

    itkSetMacro(EM_before_concentration, bool)
    itkGetMacro(EM_before_concentration, bool)

    itkSetMacro(MahalanobisThCSF, double)
    itkGetMacro(MahalanobisThCSF, double)

    itkSetMacro(MahalanobisThGM, double)
    itkGetMacro(MahalanobisThGM, double)

    itkSetMacro(MahalanobisThWM, double)
    itkGetMacro(MahalanobisThWM, double)

    itkSetMacro(FuzzyRuleMin, double)
    itkGetMacro(FuzzyRuleMin, double)

    itkSetMacro(FuzzyRuleMax, double)
    itkGetMacro(FuzzyRuleMax, double)

    itkSetMacro(UseT2, bool)
    itkGetMacro(UseT2, bool)

    itkSetMacro(UseDP, bool)
    itkGetMacro(UseDP, bool)

    itkSetMacro(UseFLAIR, bool)
    itkGetMacro(UseFLAIR, bool)

    itkSetMacro(Alpha, double)
    itkGetMacro(Alpha, double)

    itkSetMacro(MultiVarSources, double)
    itkGetMacro(MultiVarSources, double)

    itkSetMacro(MultiVarSinks, double)
    itkGetMacro(MultiVarSinks, double)

    itkSetMacro(LesionPriorProportion, double)
    itkGetMacro(LesionPriorProportion, double)

    itkSetMacro(Sigma, double)
    itkGetMacro(Sigma, double)

    itkSetMacro(UseSpecGrad, bool)
    itkGetMacro(UseSpecGrad, bool)

    itkSetMacro(MinLesionsSize, double)
    itkGetMacro(MinLesionsSize, double)

    itkSetMacro(RemoveBorder, bool)
    itkGetMacro(RemoveBorder, bool)

    itkSetMacro(IntensityT2Factor, double)
    itkGetMacro(IntensityT2Factor, double)

    itkSetMacro(IntensityDPFactor, double)
    itkGetMacro(IntensityDPFactor, double)

    itkSetMacro(IntensityFLAIRFactor, double)
    itkGetMacro(IntensityFLAIRFactor, double)

    itkSetMacro(RatioContourWM, double)
    itkGetMacro(RatioContourWM, double)

    itkSetMacro(Verbose, bool)
    itkGetMacro(Verbose, bool)

    itkSetMacro(ThresoldWMmap, double)
    itkGetMacro(ThresoldWMmap, double)

    LesionSegmentationType GetLesionSegmentationType() {return m_LesionSegmentationType;}
    void SetLesionSegmentationType(LesionSegmentationType type) {m_LesionSegmentationType=type;}

    std::string GetOutputLesionFile() {return m_OutputLesionFilename;}
    void SetOutputLesionFilename(std::string fn) {m_OutputLesionFilename=fn;}

    std::string GetOutputCSFFilename() {return m_OutputCSFFilename;}
    void SetOutputCSFFilename(std::string fn) {m_OutputCSFFilename=fn;}

    std::string GetOutputGMFilename() {return m_OutputGMFilename;}
    void SetOutputGMFilename(std::string fn) {m_OutputGMFilename=fn;}

    std::string GetOutputWMFilename() {return m_OutputWMFilename;}
    void SetOutputWMFilename(std::string fn) {m_OutputWMFilename=fn;}

    std::string GetOutputWholeFilename() {return m_OutputWholeFilename;}
    void SetOutputWholeFilename(std::string fn) {m_OutputWholeFilename=fn;}

    std::string GetOutputGCFilename() {return m_OutputGCFilename;}
    void SetOutputGCFilename(std::string fn) {m_OutputGCFilename=fn;}

    std::string GetOutputStremFilename() {return m_OutputStremFilename;}
    void SetOutputStremFilename(std::string fn) {m_OutputStremFilename=fn;}

    std::string GetOutputStremCSFFilename() {return m_OutputStremCSFFilename;}
    void SetOutputStremCSFFilename(std::string fn) {m_OutputStremCSFFilename=fn;}

    std::string GetOutputStremGMFilename() {return m_OutputStremGMFilename;}
    void SetOutputStremGMFilename(std::string fn) {m_OutputStremGMFilename=fn;}

    std::string GetOutputStremWMFilename() {return m_OutputStremWMFilename;}
    void SetOutputStremWMFilename(std::string fn) {m_OutputStremWMFilename=fn;}

    std::string GetFuzzyOutputObjectFilename() {return m_OutputFuzzyObjectFilename;}
    void SetOutputFuzzyObjectFilename(std::string fn) {m_OutputFuzzyObjectFilename=fn;}

    std::string GetOutputMahaCSFFilename() {return m_OutputMahaCSFFilename;}
    void SetOutputMahaCSFFilename(std::string fn) {m_OutputMahaCSFFilename=fn;}

    std::string GetOutputMahaGMFilename() {return m_OutputMahaGMFilename;}
    void SetOutputMahaGMFilename(std::string fn) {m_OutputMahaGMFilename=fn;}

    std::string GetOutputMahaWMFilename() {return m_OutputMahaWMFilename;}
    void SetOutputMahaWMFilename(std::string fn) {m_OutputMahaWMFilename=fn;}

    std::string GetOutputMahaMaximumFilename() {return m_OutputMahaMaximumFilename;}
    void SetOutputMahaMaximumFilename(std::string fn) {m_OutputMahaMaximumFilename=fn;}

    std::string GetOutputMahaMinimumFilename() {return m_OutputMahaMinimumFilename;}
    void SetOutputMahaMinimumFilename(std::string fn) {m_OutputMahaMinimumFilename=fn;}

    std::string GetOutputIntensityImage1Filename() {return m_OutputIntensityImage1Filename;}
    void SetOutputIntensityImage1Filename(std::string fn) {m_OutputIntensityImage1Filename=fn;}

    std::string GetOutputIntensityImage2Filename() {return m_OutputIntensityImage2Filename;}
    void SetOutputIntensityImage2Filename(std::string fn) {m_OutputIntensityImage2Filename=fn;}

    void SetGaussianModel(std::vector<GaussianFunctionType::Pointer> solution) {m_GaussianModel=solution;}
    std::vector<GaussianFunctionType::Pointer> GetSolution() {return m_GaussianModel;}

    void SetSolutionReadFilename(std::string fn) {m_SolutionReadFilename=fn;}
    void SetSolutionWriteFilename(std::string fn) {m_SolutionWriteFilename=fn;}

    std::string GetMatrixGradFilename() {return m_MatrixGradFilename;}
    void SetMatrixGradFilename(std::string fn) {m_MatrixGradFilename=fn;}

protected:
    typename GraphCutFilterType::Pointer m_GraphCutFilter;
    typename TLinksFilterType::Pointer m_TLinksFilter;
    MaximumFilterType::Pointer m_FilterMaxSources;
    MaximumFilterType::Pointer m_FilterMaxSinks;
    ComputeMahalanobisImagesFilterType::Pointer m_MahalanobisFilter;

    GcStremMsLesionsSegmentationFilter()
    {
        this->SetNumberOfRequiredOutputs(18);
        this->SetNumberOfRequiredInputs(1);

        this->SetNthOutput( 0, this->MakeOutput(0) );
        this->SetNthOutput( 1, this->MakeOutput(1) );
        this->SetNthOutput( 2, this->MakeOutput(2) );
        this->SetNthOutput( 3, this->MakeOutput(3) );
        this->SetNthOutput( 4, this->MakeOutput(4) );

        this->SetNthOutput( 5, this->MakeOutput(5) );
        this->SetNthOutput( 6, this->MakeOutput(6) );
        this->SetNthOutput( 7, this->MakeOutput(7) );
        this->SetNthOutput( 8, this->MakeOutput(8) );
        this->SetNthOutput( 9, this->MakeOutput(9) );

        this->SetNthOutput( 10, this->MakeOutput(10) );
        this->SetNthOutput( 11, this->MakeOutput(11) );
        this->SetNthOutput( 12, this->MakeOutput(12) );
        this->SetNthOutput( 13, this->MakeOutput(13) );
        this->SetNthOutput( 14, this->MakeOutput(14) );

        this->SetNthOutput( 15, this->MakeOutput(15) );
        this->SetNthOutput( 16, this->MakeOutput(16) );
        this->SetNthOutput( 17, this->MakeOutput(17) );

        m_GraphCutFilter = GraphCutFilterType::New();
        m_TLinksFilter = TLinksFilterType::New();
        m_FilterMaxSources = MaximumFilterType::New();
        m_FilterMaxSinks = MaximumFilterType::New();
        m_MahalanobisFilter = ComputeMahalanobisImagesFilterType::New();

        m_InitMethodType = 1;
        m_RejRatioHierar = 0.01;

        m_EmIter = 100;
        m_MinDistance = 0.0001;

        m_RejRatio = 0.20;
        m_EmIter_concentration = 100;
        m_EM_before_concentration = false;

        m_MahalanobisThCSF = 0.5;
        m_MahalanobisThGM = 0.5;
        m_MahalanobisThWM = 0.5;

        m_FuzzyRuleMin = 1;
        m_FuzzyRuleMax = 2;

        m_Alpha = 10;
        m_UseSpecGrad = true;
        m_Sigma = 0.6;
        m_MultiVarSources = 1;
        m_MultiVarSinks = 1;
        m_LesionPriorProportion = 0.23;

        m_MinLesionsSize = 0;
        m_RemoveBorder = false;
        m_IntensityT2Factor = 0;
        m_IntensityDPFactor = 0;
        m_IntensityFLAIRFactor = 0;
        m_RatioContourWM = 0;

        m_LabelCSF = 1;
        m_LabelGM = 2;
        m_LabelWM = 3;
        m_LabelLesions = 4;

        m_Verbose = false;

        m_UseT2 = false;
        m_UseDP = false;
        m_UseFLAIR = false;

        m_NbInputs = 0;
        m_Modalities = 0;
        m_MaxNumberOfInputs = 14;

        m_LesionSegmentationType = strem;

        m_ThresoldWMmap = 0.2;
        m_IndexWMinModel = 2;

        m_Tol = 0.0001;

        m_IndexImageT1 = m_MaxNumberOfInputs;
        m_IndexImageT2 = m_MaxNumberOfInputs;
        m_IndexImageDP = m_MaxNumberOfInputs;
        m_IndexImageFLAIR = m_MaxNumberOfInputs;
        m_IndexImageT1Gd = m_MaxNumberOfInputs;

        m_IndexMask = m_MaxNumberOfInputs;
        m_IndexAtlasCSF = m_MaxNumberOfInputs;
        m_IndexAtlasGM = m_MaxNumberOfInputs;
        m_IndexAtlasWM = m_MaxNumberOfInputs;
        m_IndexLesionPrior = m_MaxNumberOfInputs;
        m_IndexSourcesProba = m_MaxNumberOfInputs;
        m_IndexSinksProba = m_MaxNumberOfInputs;
        m_IndexSourcesMask = m_MaxNumberOfInputs;
        m_IndexSinksMask = m_MaxNumberOfInputs;

        this->SetNumberOfWorkUnits(itk::MultiThreaderBase::GetGlobalDefaultNumberOfThreads());
    }

    ~GcStremMsLesionsSegmentationFilter()
    {
    }

    itk::DataObject::Pointer MakeOutput(unsigned int idx);

    typename InputImageType::ConstPointer GetInputImageT1();
    typename InputImageType::ConstPointer GetInputImageT2();
    typename InputImageType::ConstPointer GetInputImageDP();
    typename InputImageType::ConstPointer GetInputImageFLAIR();
    typename InputImageType::ConstPointer GetInputImageT1Gd();

    typename InputImageType::ConstPointer GetMask();

    ImageTypeD::ConstPointer GetInputCSFAtlas();
    ImageTypeD::ConstPointer GetInputGMAtlas();
    ImageTypeD::ConstPointer GetInputWMAtlas();
    ImageTypeD::ConstPointer GetInputLesionPrior();

    ImageTypeUC::ConstPointer GetSourcesMask();
    ImageTypeUC::ConstPointer GetSinksMask();

    void EmitProgress(int prog);

    static void ManageProgress( itk::Object* caller, const itk::EventObject& event, void* clientData );

private:
    ITK_DISALLOW_COPY_AND_ASSIGN(GcStremMsLesionsSegmentationFilter);

    bool m_Verbose;
    double m_Tol; 
    unsigned char m_LabelLesions;
    unsigned char m_LabelCSF;
    unsigned char m_LabelGM;
    unsigned char m_LabelWM;

    unsigned int m_NbInputs, m_MaxNumberOfInputs, m_Modalities;

    unsigned int m_IndexImageT1, m_IndexImageT2, m_IndexImageDP, m_IndexImageFLAIR, m_IndexImageT1Gd, m_IndexMask;
    unsigned int m_IndexAtlasCSF, m_IndexAtlasGM, m_IndexAtlasWM, m_IndexLesionPrior;
    unsigned int m_IndexSourcesProba, m_IndexSinksProba, m_IndexSourcesMask, m_IndexSinksMask;

    LesionSegmentationType m_LesionSegmentationType;

    bool m_UseT2;
    bool m_UseDP;
    bool m_UseFLAIR;

    unsigned int m_InitMethodType;  
    double m_RejRatioHierar;        
    int m_EmIter;                   
    double m_MinDistance;           
    double m_RejRatio;              
    int m_EmIter_concentration;     
    bool m_EM_before_concentration; 
    double m_MahalanobisThCSF;       
    double m_MahalanobisThGM;        
    double m_MahalanobisThWM;        
    double m_FuzzyRuleMin;          
    double m_FuzzyRuleMax;          
    std::vector<double> m_Alphas;
    std::vector<GaussianFunctionType::Pointer> m_GaussianModel;

    bool m_UseSpecGrad;
    double m_Sigma;
    std::vector<double> m_Matrix;    
    std::string m_MatrixGradFilename;
    double m_Alpha;
    double m_MultiVarSources;
    double m_MultiVarSinks;
    double m_LesionPriorProportion; 
    double m_MinLesionsSize;         
    bool m_RemoveBorder;           
    double m_IntensityT2Factor;          
    double m_IntensityDPFactor;          
    double m_IntensityFLAIRFactor;       
    double m_HyperIntensityThreshold1;          
    double m_HyperIntensityThreshold2;       
    unsigned int m_IndexWMinModel;

    double m_RatioContourWM;              
    double m_ThresoldWMmap;

    std::string m_SolutionReadFilename;
    std::string m_SolutionWriteFilename;

    std::string m_OutputLesionFilename;
    std::string m_OutputCSFFilename;
    std::string m_OutputGMFilename;
    std::string m_OutputWMFilename;
    std::string m_OutputWholeFilename;

    std::string m_OutputGCFilename;
    std::string m_OutputStremFilename;
    std::string m_OutputStremCSFFilename;
    std::string m_OutputStremGMFilename;
    std::string m_OutputStremWMFilename;

    std::string m_OutputFuzzyObjectFilename;

    std::string m_OutputMahaCSFFilename;
    std::string m_OutputMahaGMFilename;
    std::string m_OutputMahaWMFilename;
    std::string m_OutputMahaMaximumFilename;
    std::string m_OutputMahaMinimumFilename;

    std::string m_OutputIntensityImage1Filename;
    std::string m_OutputIntensityImage2Filename;

    ImageTypeUC::Pointer m_InputImage_T1_UC;
    ImageTypeUC::Pointer m_InputImage_1_UC;
    ImageTypeUC::Pointer m_InputImage_2_UC;
    ImageTypeUC::Pointer m_MaskUC;
    ImageTypeUC::Pointer m_LesionsDetectionImage;
    ImageTypeInt::Pointer m_LabeledLesions;
    ImageTypeD::Pointer m_FuzzyObject;
};

}

#include "animaGcStremMsLesionsSegmentationFilter.hxx"