From de2db2051ac9b17a5e3e2b7dc34f42deea3f1526 Mon Sep 17 00:00:00 2001 From: dsarrut Date: Wed, 30 Jun 2010 06:01:31 +0000 Subject: [PATCH] constraint a binary image to be at a relative position of an object --- ...tivePositionConstraintToLabelImageFilter.h | 146 +++++++ ...vePositionConstraintToLabelImageFilter.txx | 358 ++++++++++++++++++ 2 files changed, 504 insertions(+) create mode 100644 itk/clitkAddRelativePositionConstraintToLabelImageFilter.h create mode 100644 itk/clitkAddRelativePositionConstraintToLabelImageFilter.txx diff --git a/itk/clitkAddRelativePositionConstraintToLabelImageFilter.h b/itk/clitkAddRelativePositionConstraintToLabelImageFilter.h new file mode 100644 index 0000000..83400ef --- /dev/null +++ b/itk/clitkAddRelativePositionConstraintToLabelImageFilter.h @@ -0,0 +1,146 @@ +/*========================================================================= + Program: vv http://www.creatis.insa-lyon.fr/rio/vv + + Authors belong to: + - University of LYON http://www.universite-lyon.fr/ + - Léon Bérard cancer center http://oncora1.lyon.fnclcc.fr + - CREATIS CNRS laboratory http://www.creatis.insa-lyon.fr + + This software is distributed WITHOUT ANY WARRANTY; without even + the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR + PURPOSE. See the copyright notices for more information. + + It is distributed under dual licence + + - BSD See included LICENSE.txt file + - CeCILL-B http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html + ======================================================================-====*/ + +#ifndef CLITKADDRELATIVEPOSITIONCONSTRAINTTOLABELIMAGEFILTER_H +#define CLITKADDRELATIVEPOSITIONCONSTRAINTTOLABELIMAGEFILTER_H + +// clitk +#include "clitkFilterBase.h" + +// itk +#include "itkPasteImageFilter.h" + +// itk ENST +#include "RelativePositionPropImageFilter.h" + +namespace clitk { + + //-------------------------------------------------------------------- + /* + Let A be an initial label image. + Let B be a label image with an object. + Let o be an orientation relatively to the B object (for example RightTo, AntTo, InferiorTo ...) + + This filter removes (=set background) from A all points that are + not in the wanted o orientation. It uses downsampled version for + faster processing, and (try to) take into account discretization + problem. Uses [Bloch 1999]. + */ + //-------------------------------------------------------------------- + + template + class ITK_EXPORT AddRelativePositionConstraintToLabelImageFilter: + public clitk::FilterBase, + public itk::ImageToImageFilter + { + + public: + /** Standard class typedefs. */ + typedef itk::ImageToImageFilter Superclass; + typedef AddRelativePositionConstraintToLabelImageFilter Self; + typedef itk::SmartPointer Pointer; + typedef itk::SmartPointer ConstPointer; + + /** Method for creation through the object factory. */ + itkNewMacro(Self); + + /** Run-time type information (and related methods). */ + itkTypeMacro(AddRelativePositionConstraintToLabelImageFilter, ImageToImageFilter); + FILTERBASE_INIT; + + /** Some convenient typedefs. */ + typedef TImageType ImageType; + typedef typename ImageType::ConstPointer ImageConstPointer; + typedef typename ImageType::Pointer ImagePointer; + typedef typename ImageType::RegionType RegionType; + typedef typename ImageType::PixelType PixelType; + typedef typename ImageType::SpacingType SpacingType; + typedef typename ImageType::SizeType SizeType; + + /** ImageDimension constants */ + itkStaticConstMacro(ImageDimension, unsigned int, TImageType::ImageDimension); + typedef itk::Image FloatImageType; + + /** Orientation types */ + typedef enum { RightTo = 0, LeftTo = 1, + AntTo = 2, PostTo = 3, + InfTo = 4, SupTo = 5, Angle = 6 + } OrientationTypeEnumeration; + + /** Input : initial image and object */ + void SetInput(const ImageType * image); + void SetInputObject(const ImageType * image); + + // Options + void SetOrientationType(OrientationTypeEnumeration orientation); + itkGetConstMacro(OrientationType, OrientationTypeEnumeration); + void SetAngle1(double a); + void SetAngle2(double a); + itkGetConstMacro(Angle1, double); + itkGetConstMacro(Angle2, double); + itkGetConstMacro(ResampleBeforeRelativePositionFilter, bool); + itkSetMacro(ResampleBeforeRelativePositionFilter, bool); + itkBooleanMacro(ResampleBeforeRelativePositionFilter); + itkGetConstMacro(IntermediateSpacing, double); + itkSetMacro(IntermediateSpacing, double); + itkGetConstMacro(FuzzyThreshold, double); + itkSetMacro(FuzzyThreshold, double); + itkGetConstMacro(BackgroundValue, PixelType); + itkSetMacro(BackgroundValue, PixelType); + itkGetConstMacro(ObjectBackgroundValue, PixelType); + itkSetMacro(ObjectBackgroundValue, PixelType); + + protected: + AddRelativePositionConstraintToLabelImageFilter(); + virtual ~AddRelativePositionConstraintToLabelImageFilter() {} + + OrientationTypeEnumeration m_OrientationType; + double m_IntermediateSpacing; + double m_FuzzyThreshold; + PixelType m_BackgroundValue; + PixelType m_ObjectBackgroundValue; + double m_Angle1; + double m_Angle2; + bool m_ResampleBeforeRelativePositionFilter; + + virtual void GenerateOutputInformation(); + virtual void GenerateInputRequestedRegion(); + virtual void GenerateData(); + + typedef itk::PasteImageFilter PadFilterType; + typename ImageType::Pointer working_image; + typename ImageType::Pointer object_resampled; + typename FloatImageType::Pointer relPos; + ImagePointer input; + ImagePointer object; + + private: + AddRelativePositionConstraintToLabelImageFilter(const Self&); //purposely not implemented + void operator=(const Self&); //purposely not implemented + + }; // end class + //-------------------------------------------------------------------- + +} // end namespace clitk +//-------------------------------------------------------------------- + +#ifndef ITK_MANUAL_INSTANTIATION +#include "clitkAddRelativePositionConstraintToLabelImageFilter.txx" +#endif + +#endif diff --git a/itk/clitkAddRelativePositionConstraintToLabelImageFilter.txx b/itk/clitkAddRelativePositionConstraintToLabelImageFilter.txx new file mode 100644 index 0000000..4078b84 --- /dev/null +++ b/itk/clitkAddRelativePositionConstraintToLabelImageFilter.txx @@ -0,0 +1,358 @@ +/*========================================================================= + Program: vv http://www.creatis.insa-lyon.fr/rio/vv + + Authors belong to: + - University of LYON http://www.universite-lyon.fr/ + - Léon Bérard cancer center http://oncora1.lyon.fnclcc.fr + - CREATIS CNRS laboratory http://www.creatis.insa-lyon.fr + + This software is distributed WITHOUT ANY WARRANTY; without even + the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR + PURPOSE. See the copyright notices for more information. + + It is distributed under dual licence + + - BSD See included LICENSE.txt file + - CeCILL-B http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html + ======================================================================-====*/ + +// clitk +#include "clitkCommon.h" +#include "clitkBooleanOperatorLabelImageFilter.h" +#include "clitkAutoCropFilter.h" +#include "clitkResampleImageWithOptionsFilter.h" +#include "clitkBooleanOperatorLabelImageFilter.h" + +// itk +#include +#include "itkStatisticsLabelMapFilter.h" +#include "itkLabelImageToStatisticsLabelMapFilter.h" +#include "itkRegionOfInterestImageFilter.h" +#include "itkBinaryThresholdImageFilter.h" +#include "itkBinaryErodeImageFilter.h" +#include "itkBinaryBallStructuringElement.h" + +// itk [Bloch et al] +#include "RelativePositionPropImageFilter.h" + +//-------------------------------------------------------------------- +template +clitk::AddRelativePositionConstraintToLabelImageFilter:: +AddRelativePositionConstraintToLabelImageFilter(): + clitk::FilterBase(), + itk::ImageToImageFilter() +{ + this->SetNumberOfRequiredInputs(2); + SetFuzzyThreshold(0.6); + SetBackgroundValue(0); + SetObjectBackgroundValue(0); + SetOrientationType(LeftTo); + ResampleBeforeRelativePositionFilterOn(); + SetIntermediateSpacing(10); + + // Step 1 : resample (option=sampling) + // Step 2 : Pad (no) + // Step 3 : relative position (option = angle) + // Step 4 : Threshold + // Step 5 : Erode for boundary + // Step 6 : resample to initial spacing + // Step 7 : pad if not the same size : it can occur when downsample and upsample + // Step 6: combine input+thresholded relpos + // Step 7: autocrop + + // Step 1 : resample + ResampleBeforeRelativePositionFilterOn(); + SetIntermediateSpacing(10); + SetOrientationType(LeftTo); + // SegmentationStep * step = new SegmentationStep; + // step->SetName("Initial resampling and relative position map"); + // SetStep(step, &Self::ResampleAndRelativePositionMap); + + // Step 3 : threshold + postprocess + // SetFuzzyThreshold(0.4); + // step = new SegmentationStep; + // step->SetName("Map threshold and post-processing"); + // SetStep(step, &Self::MapThresholdAndPostProcessing); +} +//-------------------------------------------------------------------- + + +//-------------------------------------------------------------------- +template +void +clitk::AddRelativePositionConstraintToLabelImageFilter:: +SetInput(const ImageType * image) { + // Process object is not const-correct so the const casting is required. + this->SetNthInput(0, const_cast(image)); +} +//-------------------------------------------------------------------- + + +//-------------------------------------------------------------------- +template +void +clitk::AddRelativePositionConstraintToLabelImageFilter:: +SetInputObject(const ImageType * image) { + // Process object is not const-correct so the const casting is required. + this->SetNthInput(1, const_cast(image)); +} +//-------------------------------------------------------------------- + + +//-------------------------------------------------------------------- +template +void +clitk::AddRelativePositionConstraintToLabelImageFilter:: +GenerateOutputInformation() { + ImagePointer input = dynamic_cast(itk::ProcessObject::GetInput(0)); + ImagePointer outputImage = this->GetOutput(0); + outputImage->SetRegions(outputImage->GetLargestPossibleRegion()); +} +//-------------------------------------------------------------------- + + +//-------------------------------------------------------------------- +template +void +clitk::AddRelativePositionConstraintToLabelImageFilter:: +GenerateInputRequestedRegion() { + // Call default + itk::ImageToImageFilter::GenerateInputRequestedRegion(); + // Get input pointers and set requested region to common region + ImagePointer input1 = dynamic_cast(itk::ProcessObject::GetInput(0)); + ImagePointer input2 = dynamic_cast(itk::ProcessObject::GetInput(1)); + input1->SetRequestedRegion(input1->GetLargestPossibleRegion()); + input2->SetRequestedRegion(input2->GetLargestPossibleRegion()); +} +//-------------------------------------------------------------------- + + +//-------------------------------------------------------------------- +template +void +clitk::AddRelativePositionConstraintToLabelImageFilter:: +SetAngle1(double a) { + SetOrientationType(Angle); + m_Angle1 = a; +} +//-------------------------------------------------------------------- + + +//-------------------------------------------------------------------- +template +void +clitk::AddRelativePositionConstraintToLabelImageFilter:: +SetAngle2(double a) { + SetOrientationType(Angle); + m_Angle2 = a; +} +//-------------------------------------------------------------------- + + +//-------------------------------------------------------------------- +template +void +clitk::AddRelativePositionConstraintToLabelImageFilter:: +SetOrientationType(OrientationTypeEnumeration orientation) { + m_OrientationType = orientation; + switch (m_OrientationType) { + case LeftTo: m_Angle1 = clitk::deg2rad(0); m_Angle2 = clitk::deg2rad(0); break; + case RightTo: m_Angle1 = clitk::deg2rad(180); m_Angle2 = clitk::deg2rad(0); break; + case AntTo: m_Angle1 = clitk::deg2rad(90); m_Angle2 = clitk::deg2rad(0); break; + case PostTo: m_Angle1 = clitk::deg2rad(-90); m_Angle2 = clitk::deg2rad(0); break; + case InfTo: m_Angle1 = clitk::deg2rad(0); m_Angle2 = clitk::deg2rad(90); break; + case SupTo: m_Angle1 = clitk::deg2rad(0); m_Angle2 = clitk::deg2rad(-90); break; + case Angle: break; + } + /* A1 A2 + Left 0 0 + Right 180 0 + Ant 90 0 + Post -90 0 + Inf 0 90 + Sup 0 -90 + */ +} +//-------------------------------------------------------------------- + + +//-------------------------------------------------------------------- +template +void +clitk::AddRelativePositionConstraintToLabelImageFilter:: +GenerateData() { + // Get input pointer + input = dynamic_cast(itk::ProcessObject::GetInput(0)); + object = dynamic_cast(itk::ProcessObject::GetInput(1)); + + //-------------------------------------------------------------------- + //-------------------------------------------------------------------- + StartNewStep("Resample And Relative Position Map"); + + static const unsigned int dim = ImageType::ImageDimension; + // Step 1 : resample + if (m_ResampleBeforeRelativePositionFilter) { + typedef clitk::ResampleImageWithOptionsFilter ResampleFilterType; + typename ResampleFilterType::Pointer resampleFilter = ResampleFilterType::New(); + resampleFilter->SetInput(object); + resampleFilter->SetOutputIsoSpacing(m_IntermediateSpacing); + resampleFilter->SetGaussianFilteringEnabled(false); + // resampleFilter->SetVerboseOptions(true); + resampleFilter->Update(); + working_image = resampleFilter->GetOutput(); + } + else { + working_image = object; + } + // writeImage(working_image, "res.mhd"); + + // Step 2: object pad to input image -> we want to compute the + // relative position for each point belonging to the input image + // domain, so we have to extend (pad) the object image to fit the + // domain size + typename ImageType::Pointer output = ImageType::New(); + SizeType size; + for(unsigned int i=0; iGetLargestPossibleRegion().GetSize()[i]*input->GetSpacing()[i])/(double)working_image->GetSpacing()[i]); + } + RegionType region; + region.SetSize(size); + // output->SetLargestPossibleRegion(region); + output->SetRegions(region); + output->SetSpacing(working_image->GetSpacing()); + output->SetOrigin(input->GetOrigin()); + output->Allocate(); + output->FillBuffer(m_BackgroundValue); + typename PadFilterType::Pointer padFilter = PadFilterType::New(); + typename PadFilterType::InputImageIndexType index; + for(unsigned int i=0; iGetOrigin()[i] - input->GetOrigin()[i])/(double)m_IntermediateSpacing); + } + padFilter->SetSourceImage(working_image); + padFilter->SetDestinationImage(output); + padFilter->SetDestinationIndex(index); + padFilter->SetSourceRegion(working_image->GetLargestPossibleRegion()); + padFilter->Update(); + working_image = padFilter->GetOutput(); + + // Keep object image (with resampline and pad) + object_resampled = working_image; + // writeImage(working_image, "pad.mhd"); + + // Step 3: compute rel pos in object + typedef itk::RelativePositionPropImageFilter RelPosFilterType; + typename RelPosFilterType::Pointer relPosFilter = RelPosFilterType::New(); + relPosFilter->SetInput(working_image); + relPosFilter->SetAlpha1(m_Angle1); // xy plane + relPosFilter->SetAlpha2(m_Angle2); + relPosFilter->SetK1(M_PI/2.0); // Opening parameter, default = pi/2 + relPosFilter->SetFast(true); + relPosFilter->SetRadius(1); // seems sufficient in this cas + // relPosFilter->SetVerboseProgress(true); + relPosFilter->Update(); + relPos = relPosFilter->GetOutput(); + this->template StopCurrentStep(relPos); + + //-------------------------------------------------------------------- + //-------------------------------------------------------------------- + StartNewStep("Map Threshold And Post Processing"); + + // Step 1: threshold + typedef itk::BinaryThresholdImageFilter BinaryThresholdImageFilterType; + typename BinaryThresholdImageFilterType::Pointer thresholdFilter = BinaryThresholdImageFilterType::New(); + thresholdFilter->SetInput(relPos); + thresholdFilter->SetOutsideValue(m_BackgroundValue); + thresholdFilter->SetInsideValue(m_BackgroundValue+1); + thresholdFilter->SetLowerThreshold(m_FuzzyThreshold); + thresholdFilter->Update(); + working_image = thresholdFilter->GetOutput(); + + // Step 2 : erosion with initial mask to exclude pixels that were + // inside the resampled version and outside the original mask + typedef itk::BinaryBallStructuringElement StructuringElementType; + StructuringElementType kernel; + kernel.SetRadius(1); + kernel.CreateStructuringElement(); + typedef itk::BinaryErodeImageFilter ErodeFilterType; + typename ErodeFilterType::Pointer erodeFilter = ErodeFilterType::New(); + erodeFilter->SetInput(working_image); + erodeFilter->SetKernel(kernel); + erodeFilter->SetBackgroundValue(m_BackgroundValue); + erodeFilter->SetErodeValue(m_BackgroundValue+1); + erodeFilter->Update(); + working_image = erodeFilter->GetOutput(); + + // Step 5: resample to initial spacing + if (m_ResampleBeforeRelativePositionFilter) { + typedef clitk::ResampleImageWithOptionsFilter ResampleFilterType; + typename ResampleFilterType::Pointer resampleFilter = ResampleFilterType::New(); + resampleFilter->SetDefaultPixelValue(m_BackgroundValue); + resampleFilter->SetInput(working_image); + resampleFilter->SetOutputSpacing(input->GetSpacing()); + resampleFilter->SetGaussianFilteringEnabled(false); + // resampleFilter->SetVerboseOptions(true); + resampleFilter->Update(); + working_image = resampleFilter->GetOutput(); + } + + // writeImage(working_image, "resinitial.mhd"); + + // Pre Step 6: pad if not the same size : it can occur when downsample and upsample + if (working_image->GetLargestPossibleRegion() != input->GetLargestPossibleRegion()) { + typename ImageType::Pointer temp = ImageType::New(); + temp->CopyInformation(input); + temp->SetRegions(input->GetLargestPossibleRegion()); // Do not forget !! + temp->Allocate(); + temp->FillBuffer(m_BackgroundValue); + typename PadFilterType::Pointer padFilter2 = PadFilterType::New(); // if yes : redo relpos + padFilter2->SetSourceImage(working_image); + padFilter2->SetDestinationImage(temp); + padFilter2->SetDestinationIndex(input->GetLargestPossibleRegion().GetIndex()); + padFilter2->SetSourceRegion(working_image->GetLargestPossibleRegion()); + padFilter2->Update(); + working_image = padFilter2->GetOutput(); + } + // writeImage(working_image, "pad2.mhd"); + + // Step 6: combine input+thresholded relpos + typedef clitk::BooleanOperatorLabelImageFilter BoolFilterType; + typename BoolFilterType::Pointer combineFilter = BoolFilterType::New(); + combineFilter->SetBackgroundValue(m_BackgroundValue); + combineFilter->SetBackgroundValue1(m_BackgroundValue); + combineFilter->SetBackgroundValue2(m_BackgroundValue); + combineFilter->SetForegroundValue(m_BackgroundValue+1); + combineFilter->SetInput1(input); + combineFilter->SetInput2(working_image); + combineFilter->SetOperationType(BoolFilterType::And); + combineFilter->InPlaceOn(); + combineFilter->Update(); + working_image = combineFilter->GetOutput(); + + combineFilter = BoolFilterType::New(); + combineFilter->SetInput1(working_image); + combineFilter->SetInput2(object); + combineFilter->SetOperationType(BoolFilterType::AndNot); + combineFilter->InPlaceOn(); + combineFilter->Update(); + + working_image = combineFilter->GetOutput(); + // writeImage(working_image, "combine.mhd"); + + // Step 7: autocrop + typedef clitk::AutoCropFilter CropFilterType; + typename CropFilterType::Pointer cropFilter = CropFilterType::New(); + cropFilter->SetInput(working_image); + cropFilter->ReleaseDataFlagOff(); + cropFilter->Update(); + working_image = cropFilter->GetOutput(); + this->template StopCurrentStep(working_image); + + //-------------------------------------------------------------------- + //-------------------------------------------------------------------- + + // Final Step -> set output + this->SetNthOutput(0, working_image); + return; +} +//-------------------------------------------------------------------- + -- 2.47.1