Skip to content

Commit

Permalink
BUG and NEW: New Module and Bugs fixed
Browse files Browse the repository at this point in the history
BayesianDTI, ClusteringScalar and LSDP Modules are done. NOTE: Small bug in LSDP still occur when running 1mm images (have to used WM mask to fix it).
New Morphological Shape Constraints Module only to label adjustment in the final DTI Lesion Track Python command.
MRF, GeodesicActiveContour and FCM classifiers were removed to further implementations.
  • Loading branch information
acsenrafilho committed Oct 12, 2016
1 parent c250c74 commit 340e91c
Show file tree
Hide file tree
Showing 34 changed files with 3,777 additions and 2,100 deletions.
153 changes: 143 additions & 10 deletions BayesianDTISegmentation/BayesianDTISegmentation.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -3,6 +3,23 @@
#include "itkImageRegionConstIterator.h"
#include "itkImageRegionIterator.h"

// Utils
#include "itkRescaleIntensityImageFilter.h"
#include "itkAbsImageFilter.h"

//Histogram Matching
#include "itkHistogramMatchingImageFilter.h"

//Substract Image Filter
#include "itkSubtractImageFilter.h"

//Threshold Image Filter
#include "itkThresholdImageFilter.h"

//Sigmoid Image Enhancement
#include "itkLogisticContrastEnhancementImageFilter.h"
#include "itkSigmoidImageFilter.h"

//Segmentation Methods
#include "itkBayesianClassifierInitializationImageFilter.h"
#include "itkBayesianClassifierImageFilter.h"
Expand Down Expand Up @@ -59,15 +76,130 @@ int DoIt( int argc, char * argv[], T )
typedef itk::ImageFileReader<InputImageType> ReaderType;
typedef itk::ImageFileReader<VectorInputImageType> PriorsReaderType;

typename ReaderType::Pointer reader = ReaderType::New();
reader->SetFileName( inputVolume.c_str() );
reader->Update();
typename ReaderType::Pointer inputReader = ReaderType::New();
inputReader->SetFileName( inputVolume.c_str() );
inputReader->Update();

typename ReaderType::Pointer referenceReader = ReaderType::New();

string referenceTemplate = "";
if ((mapType == "FractionalAnisotropy") & (mapResolution == "1mm")) {
referenceTemplate="USP-ICBM-FAmean-131-1mm.nii.gz";
}else if ((mapType == "MeanDiffusivity") & (mapResolution == "1mm")) {
referenceTemplate="USP-ICBM-MDmean-131-1mm.nii.gz";
}else if ((mapType == "RelativeAnisotropy") & (mapResolution == "1mm")) {
referenceTemplate="USP-ICBM-RAmean-131-1mm.nii.gz";
}else if ((mapType == "PerpendicularDiffusivity") & (mapResolution == "1mm")) {
referenceTemplate="USP-ICBM-PerpDiffmean-131-1mm.nii.gz";
}else if ((mapType == "VolumeRatio") & (mapResolution == "1mm")) {
referenceTemplate="USP-ICBM-VRmean-131-1mm.nii.gz";
}else if ((mapType == "FractionalAnisotropy") & (mapResolution == "2mm")) {
referenceTemplate="USP-ICBM-FAmean-131-2mm.nii.gz";
}else if ((mapType == "MeanDiffusivity") & (mapResolution == "2mm")) {
referenceTemplate="USP-ICBM-MDmean-131-2mm.nii.gz";
}else if ((mapType == "RelativeAnisotropy") & (mapResolution == "2mm")) {
referenceTemplate="USP-ICBM-RAmean-131-2mm.nii.gz";
}else if ((mapType == "PerpendicularDiffusivity") & (mapResolution == "2mm")) {
referenceTemplate="USP-ICBM-PerpDiffmean-131-2mm.nii.gz";
}else if ((mapType == "VolumeRatio") & (mapResolution == "2mm")) {
referenceTemplate="USP-ICBM-VRmean-131-2mm.nii.gz";
}

// Read the DTI Template
stringstream referenceTEMPLATE_path;
referenceTEMPLATE_path<<HOME_DIR<<STATISTICALTEMPLATESFOLDER<<PATH_SEPARATOR<<referenceTemplate;
referenceReader->SetFileName(referenceTEMPLATE_path.str().c_str());
referenceReader->Update();
std::cout<<"Reference image: "<<referenceTemplate<<std::endl;

//Histogram matching step
typedef itk::HistogramMatchingImageFilter<InputImageType, InputImageType> HistogramMatchType;
typename HistogramMatchType::Pointer histogramMatch = HistogramMatchType::New();
histogramMatch->SetSourceImage(inputReader->GetOutput());
histogramMatch->SetReferenceImage(referenceReader->GetOutput());
histogramMatch->SetNumberOfHistogramLevels(128);
histogramMatch->SetNumberOfMatchPoints(10000);

//Image subtraction with the DTI atlas
typedef itk::SubtractImageFilter<InputImageType,InputImageType, InputImageType> SubtractType;
typename SubtractType::Pointer subtract = SubtractType::New();
subtract->SetInput1(referenceReader->GetOutput());
subtract->SetInput2(histogramMatch->GetOutput());

//Cleaning unrelated differences from the difference image
typedef itk::ThresholdImageFilter<InputImageType> ThresholderType;
typename ThresholderType::Pointer cleanVoxels = ThresholderType::New();
cleanVoxels->SetInput(subtract->GetOutput());
if (mapType=="FractionalAnisotropy") {
cleanVoxels->SetLower(static_cast<InputPixelType>(0.0));
cleanVoxels->SetUpper(static_cast<InputPixelType>(1.0));
cleanVoxels->SetOutsideValue(static_cast<InputPixelType>(0.0));
}else if (mapType=="MeanDiffusivity") {
cleanVoxels->SetLower(static_cast<InputPixelType>(-1.0));
cleanVoxels->SetUpper(static_cast<InputPixelType>(0.0));
cleanVoxels->SetOutsideValue(static_cast<InputPixelType>(0.0));
}else if (mapType=="RelativeAnisotropy") {
cleanVoxels->SetLower(static_cast<InputPixelType>(0.0));
cleanVoxels->SetUpper(static_cast<InputPixelType>(1.0));
cleanVoxels->SetOutsideValue(static_cast<InputPixelType>(0.0));
}else if (mapType=="ParallelDiffusivity") {
cleanVoxels->SetLower(static_cast<InputPixelType>(-1.0));
cleanVoxels->SetUpper(static_cast<InputPixelType>(0.0));
cleanVoxels->SetOutsideValue(static_cast<InputPixelType>(0.0));
}else if (mapType=="VolumeRatio") {
cleanVoxels->SetLower(static_cast<InputPixelType>(-1.0));
cleanVoxels->SetUpper(static_cast<InputPixelType>(0.0));
cleanVoxels->SetOutsideValue(static_cast<InputPixelType>(0.0));
}
cleanVoxels->Update();

//Return absolute pixel values image
typedef itk::AbsImageFilter<InputImageType, InputImageType> AbsoluteImageType;
typename AbsoluteImageType::Pointer abs = AbsoluteImageType::New();
abs->SetInput(cleanVoxels->GetOutput());

//Rescale pixel intensity to [0,1] range
typedef itk::RescaleIntensityImageFilter<InputImageType, InputImageType> RescalerType;
typename RescalerType::Pointer rescaler = RescalerType::New();
rescaler->SetOutputMinimum(0.0);
rescaler->SetOutputMaximum(1.0);
rescaler->SetInput(abs->GetOutput());

//Find Sigmoid optimum parameters
typedef itk::LogisticContrastEnhancementImageFilter<InputImageType, InputImageType> OptimumSigmoidParametersType;
typename OptimumSigmoidParametersType::Pointer optSigmoid = OptimumSigmoidParametersType::New();
optSigmoid->SetInput(rescaler->GetOutput());
if (thrMethod == "MaxEntropy") {
optSigmoid->SetThresholdMethod(OptimumSigmoidParametersType::MAXENTROPY);
}else if (thrMethod == "Otsu") {
optSigmoid->SetThresholdMethod(OptimumSigmoidParametersType::OTSU);
}else if (thrMethod == "Renyi") {
optSigmoid->SetThresholdMethod(OptimumSigmoidParametersType::RENYI);
}else if (thrMethod == "Moments") {
optSigmoid->SetThresholdMethod(OptimumSigmoidParametersType::MOMENTS);
}else if (thrMethod == "IsoData") {
optSigmoid->SetThresholdMethod(OptimumSigmoidParametersType::ISODATA);
}else if (thrMethod == "Yen") {
optSigmoid->SetThresholdMethod(OptimumSigmoidParametersType::YEN);
}
optSigmoid->Update();
std::cout<<"Optimum [Alpha,Beta] = [ "<<optSigmoid->GetAlpha()<<" , "<<optSigmoid->GetBeta()<<" ]"<<std::endl;

//Sigmoid lesion enhancement step
typedef itk::SigmoidImageFilter<InputImageType,InputImageType> SigmoidType;
typename SigmoidType::Pointer sigmoid = SigmoidType::New();
sigmoid->SetInput(rescaler->GetOutput());
sigmoid->SetOutputMinimum(0.0);
sigmoid->SetOutputMaximum(1.0);
sigmoid->SetAlpha(optSigmoid->GetAlpha());
sigmoid->SetBeta(optSigmoid->GetBeta());


//Bayesian Segmentation Approach
typedef itk::BayesianClassifierInitializationImageFilter< InputImageType > BayesianInitializerType;
typename BayesianInitializerType::Pointer bayesianInitializer = BayesianInitializerType::New();

bayesianInitializer->SetInput( reader->GetOutput() );
bayesianInitializer->SetInput( sigmoid->GetOutput() );
bayesianInitializer->SetNumberOfClasses( 2 ); // Always set to 2 classes: Lesion and non-Lesions probability
bayesianInitializer->Update();

Expand All @@ -78,16 +210,16 @@ int DoIt( int argc, char * argv[], T )

string priorsTemplate = "";
if (priorsImage == "Multiple Sclerosis DTI Lesions") {
if (resolution == "1mm") {
priorsTemplate="USP-ICBM-MSLesionPriors-30-1mm.nii.gz";
if (mapResolution == "1mm") {
priorsTemplate="USP-ICBM-MSLesionPriors-46-1mm.nii.gz";
}else {
priorsTemplate="USP-ICBM-MSLesionPriors-30-2mm.nii.gz";
priorsTemplate="USP-ICBM-MSLesionPriors-46-2mm.nii.gz";
}
// Read the MS lesion priors probability image
stringstream TEMPLATE_path;
TEMPLATE_path<<HOME_DIR<<STATISTICALTEMPLATESFOLDER<<PATH_SEPARATOR<<priorsTemplate;
stringstream priors_path;
priors_path<<HOME_DIR<<STATISTICALTEMPLATESFOLDER<<PATH_SEPARATOR<<priorsTemplate;
typename PriorsReaderType::Pointer priorsReader = PriorsReaderType::New();
priorsReader->SetFileName(TEMPLATE_path.str().c_str());
priorsReader->SetFileName(priors_path.str().c_str());
priorsReader->Update();

bayesClassifier->SetPriors(priorsReader->GetOutput());
Expand All @@ -102,6 +234,7 @@ int DoIt( int argc, char * argv[], T )
writer->SetUseCompression(1);
writer->Update();


return EXIT_SUCCESS;
}

Expand Down
45 changes: 35 additions & 10 deletions BayesianDTISegmentation/BayesianDTISegmentation.xml
Original file line number Diff line number Diff line change
Expand Up @@ -15,16 +15,37 @@
<name>inputVolume</name>
<label>Input Volume</label>
<channel>input</channel>
<index>0</index>
<index>1</index>
<description><![CDATA[Input volume]]></description>
</image>
<image type="label">
<name>outputLabel</name>
<label>Output Label</label>
<channel>output</channel>
<index>1</index>
<index>2</index>
<description><![CDATA[Output Label]]></description>
</image>
<string-enumeration>
<name>mapType</name>
<longflag>--diffMap</longflag>
<description><![CDATA[Choose the diffusion map inserted in the input image field.]]></description>
<label>Diffusion Map Type</label>
<default>FractionalAnisotropy</default>
<element>FractionalAnisotropy</element>
<element>MeanDiffusivity</element>
<element>RelativeAnisotropy</element>
<element>PerpendicularDiffusivity</element>
<element>VolumeRatio</element>
</string-enumeration>
<string-enumeration>
<name>mapResolution</name>
<longflag>--diffMapRes</longflag>
<description><![CDATA[Choose the diffusion map resolution.]]></description>
<label>Diffusion Map Resolution</label>
<default>2mm</default>
<element>2mm</element>
<element>1mm</element>
</string-enumeration>
</parameters>
<parameters>
<label>Segmentation Parameters</label>
Expand All @@ -38,14 +59,18 @@
<element>Unconditional</element>
<element>Multiple Sclerosis DTI Lesions</element>
</string-enumeration>
<string-enumeration>
<name>resolution</name>
<longflag>--res</longflag>
<description><![CDATA[Choose the image resolution.]]></description>
<label>Image Resolution</label>
<default>2mm</default>
<element>2mm</element>
<element>1mm</element>
<string-enumeration>
<name>thrMethod</name>
<longflag>--thrMethod</longflag>
<description><![CDATA[Choose the threshold method that will optimize the logistic lesion enhancment parameters.]]></description>
<label>Threshold Method</label>
<default>Otsu</default>
<element>MaxEntropy</element>
<element>Otsu</element>
<element>Yen</element>
<element>IsoData</element>
<element>Moments</element>
<element>Renyi</element>
</string-enumeration>
</parameters>
</executable>
2 changes: 2 additions & 0 deletions BayesianDTISegmentation/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -38,6 +38,8 @@ SEMMacroBuildCLI(
TARGET_LIBRARIES ${MODULE_TARGET_LIBRARIES}
INCLUDE_DIRECTORIES ${MODULE_INCLUDE_DIRECTORIES}
ADDITIONAL_SRCS ${MODULE_SRCS}
itkLogisticContrastEnhancementImageFilter.h
itkLogisticContrastEnhancementImageFilter.hxx
)

#-----------------------------------------------------------------------------
Expand Down
119 changes: 119 additions & 0 deletions BayesianDTISegmentation/itkLogisticContrastEnhancementImageFilter.h
Original file line number Diff line number Diff line change
@@ -0,0 +1,119 @@
#ifndef __itkLogisticContrastEnhancementImageFilter_h
#define __itkLogisticContrastEnhancementImageFilter_h
#include "itkImageToImageFilter.h"
#include "itkImage.h"
#include "itkNumericTraits.h"

namespace itk
{

template< typename TInputImage , typename TOutputImage>
class ITK_EXPORT LogisticContrastEnhancementImageFilter:
public ImageToImageFilter< TInputImage, TOutputImage >
{
public:
/** Extract dimension from inputs images, where it is assumed there are with the same type. */
itkStaticConstMacro(InputImageDimension, unsigned int,
TInputImage::ImageDimension);
itkStaticConstMacro(OutputImageDimension, unsigned int,
TOutputImage::ImageDimension);

/** Convenient typedefs for simplifying declarations. */
typedef TInputImage InputImageType;
typedef TOutputImage OutputImageType;


/** Standard class typedefs. */
typedef LogisticContrastEnhancementImageFilter Self;
typedef ImageToImageFilter< TInputImage, TOutputImage > Superclass;
typedef SmartPointer< Self > Pointer;
typedef SmartPointer< const Self > ConstPointer;

/** Method for creation through the object factory. */
itkNewMacro(Self)

/** Run-time type information (and related methods). */
itkTypeMacro(LogisticContrastEnhancementImageFilter, ImageToImageFilter)

typedef typename InputImageType::PixelType InputPixelType;
typedef typename OutputImageType::PixelType OutputPixelType;

/** Set the maximum output. */
itkSetMacro(MaximumOutput, double)

/** Set the minimum output. */
itkSetMacro(MinimumOutput, double)

/** Set the Histogram bins. */
itkSetMacro(NumberOfBins, unsigned int)

/** Set the Tolerance. */
itkSetMacro(Tolerance, char)

/** Set the object area. */
itkSetMacro(FlipObjectArea, bool)
itkBooleanMacro(FlipObjectArea)

/** Set if use manual tolerance. */
itkSetMacro(ManualTolerance, bool)
itkBooleanMacro(ManualTolerance)

/** Set threshold method. */
itkSetMacro(ThresholdMethod, unsigned char)

itkGetMacro(FlipObjectArea, bool)
itkGetMacro(ManualTolerance, bool)
itkGetMacro(Alpha, double)
itkGetMacro(Beta, double)
itkGetMacro(MaximumOutput, double)
itkGetMacro(MinimumOutput, double)
itkGetMacro(NumberOfBins, unsigned int)
itkGetMacro(Tolerance, char)
itkGetMacro(ThresholdMethod, unsigned char)

#ifdef ITK_USE_CONCEPT_CHECKING
// Begin concept checking
itkConceptMacro( InputHasNumericTraitsCheck,
( Concept::HasNumericTraits< InputPixelType > ) );
itkConceptMacro( SameDimensionCheck,
( Concept::SameDimension< InputImageDimension, OutputImageDimension > ) );
#endif

enum ThresholdMethod {
MAXENTROPY=1,
OTSU=2,
RENYI=3,
MOMENTS=4,
YEN=5,
ISODATA=6
};

protected:
LogisticContrastEnhancementImageFilter();
virtual ~LogisticContrastEnhancementImageFilter() {}
bool m_FlipObjectArea;
bool m_ManualTolerance;
char m_Tolerance;
double m_Alpha;
double m_Beta;
double m_MaximumOutput;
double m_MinimumOutput;
unsigned int m_NumberOfBins;
unsigned char m_ThresholdMethod;

void GenerateData();
private:
// const unsigned short MAX_TOLERANCE = 0.99;
LogisticContrastEnhancementImageFilter(const Self &); //purposely not implemented
void operator=(const Self &); //purposely not implemented
double sigmoid(double x, double alpha, double beta);
void checkTolerance(char tolerance);
};

} // end namespace itk

#ifndef ITK_MANUAL_INSTANTIATION
#include "itkLogisticContrastEnhancementImageFilter.hxx"
#endif

#endif
Loading

0 comments on commit 340e91c

Please sign in to comment.