Authors belong to:
- University of LYON http://www.universite-lyon.fr/
- - Léon Bérard cancer center http://oncora1.lyon.fnclcc.fr
+ - Léon Bérard cancer center http://www.centreleonberard.fr
- CREATIS CNRS laboratory http://www.creatis.insa-lyon.fr
This software is distributed WITHOUT ANY WARRANTY; without even
- BSD See included LICENSE.txt file
- CeCILL-B http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
- ======================================================================-====*/
+ ===========================================================================**/
// clitk
#include "clitkCommon.h"
#include "clitkAutoCropFilter.h"
#include "clitkResampleImageWithOptionsFilter.h"
#include "clitkBooleanOperatorLabelImageFilter.h"
+#include "clitkCropLikeImageFilter.h"
// itk
#include <deque>
-#include "itkStatisticsLabelMapFilter.h"
-#include "itkLabelImageToStatisticsLabelMapFilter.h"
-#include "itkRegionOfInterestImageFilter.h"
-#include "itkBinaryThresholdImageFilter.h"
-#include "itkBinaryErodeImageFilter.h"
-#include "itkBinaryBallStructuringElement.h"
+#include <itkStatisticsLabelMapFilter.h>
+#include <itkLabelImageToStatisticsLabelMapFilter.h>
+#include <itkRegionOfInterestImageFilter.h>
+#include <itkBinaryThresholdImageFilter.h>
+#include <itkBinaryErodeImageFilter.h>
+#include <itkBinaryBallStructuringElement.h>
+#include <itkAddImageFilter.h>
+#include <itkDivideImageFilter.h>
// itk [Bloch et al]
#include "RelativePositionPropImageFilter.h"
//--------------------------------------------------------------------
-template <class TImageType>
-clitk::AddRelativePositionConstraintToLabelImageFilter<TImageType>::
+template <class ImageType>
+clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
AddRelativePositionConstraintToLabelImageFilter():
clitk::FilterBase(),
- itk::ImageToImageFilter<TImageType, TImageType>()
+ itk::ImageToImageFilter<ImageType, ImageType>()
{
this->SetNumberOfRequiredInputs(2);
SetFuzzyThreshold(0.6);
SetBackgroundValue(0);
SetObjectBackgroundValue(0);
- SetOrientationType(LeftTo);
- ResampleBeforeRelativePositionFilterOn();
+ ClearOrientationType();
+ IntermediateSpacingFlagOn();
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);
+ AutoCropFlagOn();
+ InverseOrientationFlagOff();
+ RemoveObjectFlagOn();
+ CombineWithOrFlagOff();
+ VerboseStepFlagOff();
+ WriteStepFlagOff();
+ FuzzyMapOnlyFlagOff();
+ FastFlagOff();
+ SetRadius(2.0);
+ SetK1(std::acos(-1.0)/2);
}
//--------------------------------------------------------------------
//--------------------------------------------------------------------
-template <class TImageType>
+template <class ImageType>
void
-clitk::AddRelativePositionConstraintToLabelImageFilter<TImageType>::
-SetInput(const ImageType * image) {
+clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
+SetInput(const ImageType * image)
+{
// Process object is not const-correct so the const casting is required.
this->SetNthInput(0, const_cast<ImageType *>(image));
}
//--------------------------------------------------------------------
-template <class TImageType>
+template <class ImageType>
void
-clitk::AddRelativePositionConstraintToLabelImageFilter<TImageType>::
-SetInputObject(const ImageType * image) {
+clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
+SetInputObject(const ImageType * image)
+{
// Process object is not const-correct so the const casting is required.
this->SetNthInput(1, const_cast<ImageType *>(image));
}
//--------------------------------------------------------------------
-template <class TImageType>
+template <class ImageType>
+void
+clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
+ClearOrientationType()
+{
+ m_OrientationTypeString.clear();
+ m_OrientationType.clear();
+ m_Angle1.clear();
+ m_Angle2.clear();
+}
+//--------------------------------------------------------------------
+
+
+//--------------------------------------------------------------------
+template <class ImageType>
+int
+clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
+GetNumberOfAngles()
+{
+ return m_OrientationType.size();
+}
+//--------------------------------------------------------------------
+
+
+//--------------------------------------------------------------------
+template <class ImageType>
+void
+clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
+AddOrientationTypeString(std::string t)
+{
+ m_OrientationTypeString.push_back(t);
+
+ if (t == "LeftTo") { AddOrientationType(LeftTo); return; }
+ if (t == "RightTo") { AddOrientationType(RightTo); return; }
+ if (t == "AntTo") { AddOrientationType(AntTo); return; }
+ if (t == "PostTo") { AddOrientationType(PostTo); return; }
+ if (t == "SupTo") { AddOrientationType(SupTo); return; }
+ if (t == "InfTo") { AddOrientationType(InfTo); return; }
+
+ if (t == "NotLeftTo") { AddOrientationType(LeftTo); InverseOrientationFlagOn(); return; }
+ if (t == "NotRightTo") { AddOrientationType(RightTo); InverseOrientationFlagOn(); return; }
+ if (t == "NotAntTo") { AddOrientationType(AntTo); InverseOrientationFlagOn(); return; }
+ if (t == "NotPostTo") { AddOrientationType(PostTo); InverseOrientationFlagOn(); return; }
+ if (t == "NotSupTo") { AddOrientationType(SupTo); InverseOrientationFlagOn(); return; }
+ if (t == "NotInfTo") { AddOrientationType(InfTo); InverseOrientationFlagOn(); return; }
+
+ if (t == "Angle") return;
+
+ clitkExceptionMacro("Error, you must provide LeftTo,RightTo or AntTo,PostTo or SupTo,InfTo (or NotLeftTo, NotRightTo etc) but you give " << t);
+}
+//--------------------------------------------------------------------
+
+
+//--------------------------------------------------------------------
+template <class ImageType>
void
-clitk::AddRelativePositionConstraintToLabelImageFilter<TImageType>::
-GenerateOutputInformation() {
- ImagePointer input = dynamic_cast<TImageType*>(itk::ProcessObject::GetInput(0));
+clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
+GenerateOutputInformation()
+{
+ ImagePointer input = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(0));
ImagePointer outputImage = this->GetOutput(0);
outputImage->SetRegions(outputImage->GetLargestPossibleRegion());
}
//--------------------------------------------------------------------
-template <class TImageType>
+template <class ImageType>
void
-clitk::AddRelativePositionConstraintToLabelImageFilter<TImageType>::
-GenerateInputRequestedRegion() {
+clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
+GenerateInputRequestedRegion()
+{
// Call default
- itk::ImageToImageFilter<TImageType, TImageType>::GenerateInputRequestedRegion();
+ itk::ImageToImageFilter<ImageType, ImageType>::GenerateInputRequestedRegion();
// Get input pointers and set requested region to common region
- ImagePointer input1 = dynamic_cast<TImageType*>(itk::ProcessObject::GetInput(0));
- ImagePointer input2 = dynamic_cast<TImageType*>(itk::ProcessObject::GetInput(1));
+ ImagePointer input1 = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(0));
+ ImagePointer input2 = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(1));
input1->SetRequestedRegion(input1->GetLargestPossibleRegion());
input2->SetRequestedRegion(input2->GetLargestPossibleRegion());
}
//--------------------------------------------------------------------
-template <class TImageType>
+template <class ImageType>
void
-clitk::AddRelativePositionConstraintToLabelImageFilter<TImageType>::
-SetAngle1(double a) {
- SetOrientationType(Angle);
- m_Angle1 = a;
+clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
+AddAnglesInRad(double a, double b)
+{
+ m_OrientationTypeString.push_back("Angle");
+ m_OrientationType.push_back(Angle);
+ m_Angle1.push_back(a);
+ m_Angle2.push_back(b);
}
//--------------------------------------------------------------------
//--------------------------------------------------------------------
-template <class TImageType>
+template <class ImageType>
void
-clitk::AddRelativePositionConstraintToLabelImageFilter<TImageType>::
-SetAngle2(double a) {
- SetOrientationType(Angle);
- m_Angle2 = a;
+clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
+AddAnglesInDeg(double a, double b)
+{
+ AddAnglesInRad(clitk::deg2rad(a), clitk::deg2rad(b));
}
//--------------------------------------------------------------------
//--------------------------------------------------------------------
-template <class TImageType>
+template <class ImageType>
void
-clitk::AddRelativePositionConstraintToLabelImageFilter<TImageType>::
-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;
+clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
+AddOrientationType(OrientationTypeEnumeration orientation)
+{
+ m_OrientationType.push_back(orientation);
+ switch (orientation) {
+ case RightTo:
+ m_Angle1.push_back(clitk::deg2rad(0));
+ m_Angle2.push_back(clitk::deg2rad(0));
+ break;
+ case LeftTo:
+ m_Angle1.push_back(clitk::deg2rad(180));
+ m_Angle2.push_back(clitk::deg2rad(0));
+ break;
+ case AntTo:
+ m_Angle1.push_back(clitk::deg2rad(90));
+ m_Angle2.push_back(clitk::deg2rad(0));
+ break;
+ case PostTo:
+ m_Angle1.push_back(clitk::deg2rad(-90));
+ m_Angle2.push_back(clitk::deg2rad(0));
+ break;
+ case InfTo:
+ m_Angle1.push_back(clitk::deg2rad(0));
+ m_Angle2.push_back(clitk::deg2rad(90));
+ break;
+ case SupTo:
+ m_Angle1.push_back(clitk::deg2rad(0));
+ m_Angle2.push_back(clitk::deg2rad(-90));
+ break;
+ case Angle: break;
}
/* A1 A2
Left 0 0
//--------------------------------------------------------------------
-template <class TImageType>
+template <class ImageType>
void
-clitk::AddRelativePositionConstraintToLabelImageFilter<TImageType>::
-GenerateData() {
+clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
+PrintOptions()
+{
+ DD((int)this->GetBackgroundValue());
+ DD((int)this->GetObjectBackgroundValue());
+ DDV(this->GetOrientationTypeString(), (uint)this->GetNumberOfAngles());
+ DD(this->GetIntermediateSpacingFlag());
+ DD(this->GetIntermediateSpacing());
+ DD(this->GetFuzzyThreshold());
+ DD(this->GetAutoCropFlag());
+ DD(this->GetInverseOrientationFlag());
+ DD(this->GetRemoveObjectFlag());
+ DD(this->GetCombineWithOrFlag());
+}
+//--------------------------------------------------------------------
+
+
+//--------------------------------------------------------------------
+template <class ImageType>
+void
+clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
+GenerateData()
+{
+ if (GetNumberOfAngles() <1) {
+ clitkExceptionMacro("Add at least one orientation type");
+ }
+
+ if (GetVerboseOptionFlag()) {
+ for(int i=0; i<GetNumberOfAngles(); i++) {
+ std::cout << "Orientation \t:" << GetOrientationTypeString(i) << std::endl;
+ }
+ std::cout << "Interm Spacing \t:" << GetIntermediateSpacingFlag() << " " << GetIntermediateSpacing() << "mm" << std::endl;
+ std::cout << "Fuzzy threshold\t:" << GetFuzzyThreshold() << std::endl;
+ std::cout << "AutoCrop \t:" << GetAutoCropFlag() << std::endl;
+ std::cout << "InverseOrient \t:" << GetInverseOrientationFlag() << std::endl;
+ std::cout << "RemoveObject \t:" << GetRemoveObjectFlag() << std::endl;
+ std::cout << "CombineWithOr \t:" << GetCombineWithOrFlag() << std::endl;
+ }
+
// Get input pointer
input = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(0));
object = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(1));
+ static const unsigned int dim = ImageType::ImageDimension;
+
+ // 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
+ working_image = object;
+ if (!clitk::HaveSameSizeAndSpacing<ImageType, ImageType>(input, working_image)) {
+ StartNewStep("Pad (resize) object to input size");
+
+ if (0) { // OLD VERSION (TO REMOVE)
+ StartNewStep("Pad object to image size");
+ typename ImageType::Pointer output = ImageType::New();
+ SizeType size;
+ for(unsigned int i=0; i<dim; i++) {
+ size[i] = lrint((input->GetLargestPossibleRegion().GetSize()[i]*
+ input->GetSpacing()[i])/(double)working_image->GetSpacing()[i]);
+ }
+
+ // The index of the input is not necessarily zero, so we have to
+ // take it into account (not done)
+ RegionType region;
+ IndexType index = input->GetLargestPossibleRegion().GetIndex();
+ region.SetSize(size);
+ for(unsigned int i=0; i<dim; i++) {
+ if (index[i] != 0) {
+ std::cerr << "Index diff from zero : " << index << ". not done yet !" << std::endl;
+ exit(0);
+ }
+ }
+ // output->SetLargestPossibleRegion(region);
+ output->SetRegions(region);
+ output->SetSpacing(working_image->GetSpacing());
+ PointType origin = input->GetOrigin();
+ for(unsigned int i=0; i<dim; i++) {
+ origin[i] = index[i]*input->GetSpacing()[i] + input->GetOrigin()[i];
+ }
+ output->SetOrigin(origin);
+ // output->SetOrigin(input->GetOrigin());
+
+ output->Allocate();
+ output->FillBuffer(m_BackgroundValue);
+ typename PasteFilterType::Pointer padFilter = PasteFilterType::New();
+ // typename PasteFilterType::InputImageIndexType index;
+ for(unsigned int i=0; i<dim; i++) {
+ index[i] = -index[i]*input->GetSpacing()[i]/(double)working_image->GetSpacing()[i]
+ + lrint((working_image->GetOrigin()[i] - input->GetOrigin()[i])/working_image->GetSpacing()[i]);
+ }
+ padFilter->SetSourceImage(working_image);
+ padFilter->SetDestinationImage(output);
+ padFilter->SetDestinationIndex(index);
+ padFilter->SetSourceRegion(working_image->GetLargestPossibleRegion());
+ padFilter->Update();
+ working_image = padFilter->GetOutput();
+ }
+
+ // Resize object like input
+ working_image = clitk::ResizeImageLike<ImageType>(working_image, input, GetBackgroundValue());
+ StopCurrentStep<ImageType>(working_image);
+ }
//--------------------------------------------------------------------
//--------------------------------------------------------------------
- StartNewStep("Resample And Relative Position Map");
-
- static const unsigned int dim = ImageType::ImageDimension;
// Step 1 : resample
- if (m_ResampleBeforeRelativePositionFilter) {
+ if (m_IntermediateSpacingFlag) {
+ StartNewStep("Resample object to intermediate spacing (" + toString(m_IntermediateSpacing) + ")");
typedef clitk::ResampleImageWithOptionsFilter<ImageType> ResampleFilterType;
typename ResampleFilterType::Pointer resampleFilter = ResampleFilterType::New();
- resampleFilter->SetInput(object);
+ resampleFilter->SetInput(working_image);
+ resampleFilter->SetDefaultPixelValue(0);
resampleFilter->SetOutputIsoSpacing(m_IntermediateSpacing);
resampleFilter->SetGaussianFilteringEnabled(false);
- // resampleFilter->SetVerboseOptions(true);
+ //resampleFilter->SetVerboseOptions(true);
resampleFilter->Update();
working_image = resampleFilter->GetOutput();
+ StopCurrentStep<ImageType>(working_image);
}
- else {
- working_image = object;
- }
- // writeImage<ImageType>(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; i<dim; i++) {
- size[i] = lrint((input->GetLargestPossibleRegion().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; i<dim; i++) {
- index[i] = lrint((working_image->GetOrigin()[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<ImageType>(working_image, "pad.mhd");
// Step 3: compute rel pos in object
+ StartNewStep("Relative Position Map");
typedef itk::RelativePositionPropImageFilter<ImageType, FloatImageType> 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<FloatImageType>(relPos);
+ typename RelPosFilterType::Pointer relPosFilter;
+
+ for(int i=0; i<GetNumberOfAngles(); i++) {
+ // Compute fuzzy map
+ relPosFilter = RelPosFilterType::New();
+ relPosFilter->SetFast(GetFastFlag());
+ relPosFilter->SetRadius(GetRadius());
+ relPosFilter->SetInput(working_image);
+ relPosFilter->SetAlpha1(m_Angle1[i]); // xy plane
+ relPosFilter->SetAlpha2(m_Angle2[i]);
+ relPosFilter->SetK1(GetK1());// M_PI/2.0); // Opening parameter, default = pi/2
+ // relPosFilter->SetVerboseProgress(true);
+
+ relPosFilter->Update();
+ relPos = relPosFilter->GetOutput();
+
+ if (GetNumberOfAngles() != 1) {
+ // Creation of the first m_FuzzyMap
+ if (i==0) {
+ m_FuzzyMap = clitk::NewImageLike<FloatImageType>(relPos, true);
+ m_FuzzyMap->FillBuffer(0.0);
+ }
+
+ // Add to current fuzzy map
+ typedef itk::AddImageFilter<FloatImageType, FloatImageType, FloatImageType> AddImageFilter;
+ typename AddImageFilter::Pointer addFilter = AddImageFilter::New();
+ addFilter->SetInput1(m_FuzzyMap);
+ addFilter->SetInput2(relPos);
+ addFilter->Update();
+ m_FuzzyMap = addFilter->GetOutput();
+ }
+ else m_FuzzyMap = relPos;
+ }
+
+ // Divide by the number of relpos
+ if (GetNumberOfAngles() != 1) {
+ typedef itk::DivideImageFilter<FloatImageType, FloatImageType, FloatImageType> DivideFilter;
+ typename DivideFilter::Pointer divideFilter = DivideFilter::New();
+ divideFilter->SetConstant2(GetNumberOfAngles());
+ divideFilter->SetInput(m_FuzzyMap);
+ divideFilter->Update();
+ m_FuzzyMap = divideFilter->GetOutput();
+ }
+
+ relPos = m_FuzzyMap;
+ StopCurrentStep<FloatImageType>(relPos);
+ if (GetFuzzyMapOnlyFlag()) {
+ // If the spacing is used, recompute fuzzy map with the input size/spacing
+ if (m_IntermediateSpacingFlag) {
+ StartNewStep("Resample FuzzyMap to come back to the same sampling than input");
+ typedef clitk::ResampleImageWithOptionsFilter<FloatImageType> ResampleFilterType;
+ typename ResampleFilterType::Pointer resampleFilter = ResampleFilterType::New();
+ resampleFilter->SetDefaultPixelValue(0.0); //Default fuzzymap value FIXME
+ resampleFilter->SetInput(relPos);
+ resampleFilter->SetOutputSpacing(input->GetSpacing());
+ resampleFilter->SetGaussianFilteringEnabled(false);
+ resampleFilter->Update();
+ relPos = m_FuzzyMap = resampleFilter->GetOutput();
+ StopCurrentStep<FloatImageType>(relPos);
+
+ // Need to put exactly the same size
+ if (relPos->GetLargestPossibleRegion() != input->GetLargestPossibleRegion()) {
+ StartNewStep("Pad to get the same size than input");
+ typename FloatImageType::Pointer temp = FloatImageType::New();
+ temp->CopyInformation(input);
+ temp->SetRegions(input->GetLargestPossibleRegion()); // Do not forget !!
+ temp->Allocate();
+ temp->FillBuffer(0.0); // Default fuzzymap value FIXME
+ typename PasteFloatFilterType::Pointer padFilter2 = PasteFloatFilterType::New();
+ padFilter2->SetSourceImage(relPos);
+ padFilter2->SetDestinationImage(temp);
+ padFilter2->SetDestinationIndex(input->GetLargestPossibleRegion().GetIndex());
+ padFilter2->SetSourceRegion(relPos->GetLargestPossibleRegion());
+ padFilter2->Update();
+ relPos = padFilter2->GetOutput();
+ StopCurrentStep<FloatImageType>(relPos);
+ m_FuzzyMap = relPos;
+ }
+ }
+ else {
+ StopCurrentStep<FloatImageType>(relPos);
+ }
+ return;
+ }
//--------------------------------------------------------------------
//--------------------------------------------------------------------
- StartNewStep("Map Threshold And Post Processing");
-
+ StartNewStep("Map Threshold");
// Step 1: threshold
typedef itk::BinaryThresholdImageFilter<FloatImageType, ImageType> BinaryThresholdImageFilterType;
typename BinaryThresholdImageFilterType::Pointer thresholdFilter = BinaryThresholdImageFilterType::New();
thresholdFilter->SetLowerThreshold(m_FuzzyThreshold);
thresholdFilter->Update();
working_image = thresholdFilter->GetOutput();
+ StopCurrentStep<ImageType>(working_image);
+ //--------------------------------------------------------------------
+ //--------------------------------------------------------------------
+ StartNewStep("Post Processing: erosion with initial mask");
// Step 2 : erosion with initial mask to exclude pixels that were
// inside the resampled version and outside the original mask
- typedef itk::BinaryBallStructuringElement<unsigned int, 3> StructuringElementType;
+ typedef itk::BinaryBallStructuringElement<unsigned int, ImageDimension> StructuringElementType;
StructuringElementType kernel;
kernel.SetRadius(1);
kernel.CreateStructuringElement();
erodeFilter->SetErodeValue(m_BackgroundValue+1);
erodeFilter->Update();
working_image = erodeFilter->GetOutput();
+ StopCurrentStep<ImageType>(working_image);
+ //--------------------------------------------------------------------
+ //--------------------------------------------------------------------
// Step 5: resample to initial spacing
- if (m_ResampleBeforeRelativePositionFilter) {
+ if (m_IntermediateSpacingFlag) {
+ StartNewStep("Resample to come back to the same sampling than input");
typedef clitk::ResampleImageWithOptionsFilter<ImageType> ResampleFilterType;
typename ResampleFilterType::Pointer resampleFilter = ResampleFilterType::New();
resampleFilter->SetDefaultPixelValue(m_BackgroundValue);
// resampleFilter->SetVerboseOptions(true);
resampleFilter->Update();
working_image = resampleFilter->GetOutput();
+ StopCurrentStep<ImageType>(working_image);
}
- // writeImage<ImageType>(working_image, "resinitial.mhd");
-
+ //--------------------------------------------------------------------
+ //--------------------------------------------------------------------
// Pre Step 6: pad if not the same size : it can occur when downsample and upsample
+ //if (!HaveSameSizeAndSpacing(working_image, input)) {
if (working_image->GetLargestPossibleRegion() != input->GetLargestPossibleRegion()) {
+ StartNewStep("Pad to get the same size than input");
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
+ temp->FillBuffer(m_BackgroundValue);
+ typename PasteFilterType::Pointer padFilter2 = PasteFilterType::New();
padFilter2->SetSourceImage(working_image);
padFilter2->SetDestinationImage(temp);
padFilter2->SetDestinationIndex(input->GetLargestPossibleRegion().GetIndex());
padFilter2->SetSourceRegion(working_image->GetLargestPossibleRegion());
padFilter2->Update();
working_image = padFilter2->GetOutput();
+ StopCurrentStep<ImageType>(working_image);
+ }
+ else {
+ //DD("[debug] Rel Pos : no padding after");
}
- // writeImage<ImageType>(working_image, "pad2.mhd");
+ //--------------------------------------------------------------------
+ //--------------------------------------------------------------------
// Step 6: combine input+thresholded relpos
+ StartNewStep("Combine with initial input (boolean And)");
typedef clitk::BooleanOperatorLabelImageFilter<ImageType> BoolFilterType;
typename BoolFilterType::Pointer combineFilter = BoolFilterType::New();
combineFilter->SetBackgroundValue(m_BackgroundValue);
combineFilter->SetForegroundValue(m_BackgroundValue+1);
combineFilter->SetInput1(input);
combineFilter->SetInput2(working_image);
- combineFilter->SetOperationType(BoolFilterType::And);
- combineFilter->InPlaceOn();
+ if (GetInverseOrientationFlag())
+ combineFilter->SetOperationType(BoolFilterType::AndNot);
+ else {
+ if (GetCombineWithOrFlag())
+ combineFilter->SetOperationType(BoolFilterType::Or);
+ else
+ combineFilter->SetOperationType(BoolFilterType::And);
+ }
+ combineFilter->InPlaceOff(); // Do not modify initial input (!)
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<ImageType>(working_image, "combine.mhd");
+ // Remove (if needed the object from the support)
+ if (GetRemoveObjectFlag()) {
+ combineFilter = BoolFilterType::New();
+ combineFilter->SetInput1(working_image);
+ combineFilter->SetInput2(object);
+ combineFilter->SetOperationType(BoolFilterType::AndNot);
+ combineFilter->InPlaceOn();
+ combineFilter->Update();
+ working_image = combineFilter->GetOutput();
+ }
+ // In the other case, we must *add* the initial object to keep it
+ // but not more than the initial support
+ else {
+ combineFilter = BoolFilterType::New();
+ combineFilter->SetInput1(working_image);
+ combineFilter->SetInput2(object);
+ combineFilter->SetOperationType(BoolFilterType::Or);
+ combineFilter->InPlaceOn();
+ combineFilter->Update();
+ working_image = combineFilter->GetOutput(); // not needed because InPlaceOn ?
+ combineFilter = BoolFilterType::New();
+ combineFilter->SetInput1(working_image);
+ combineFilter->SetInput2(input);
+ combineFilter->SetOperationType(BoolFilterType::And);
+ combineFilter->InPlaceOn();
+ combineFilter->Update();
+ working_image = combineFilter->GetOutput();
+ }
+
+ StopCurrentStep<ImageType>(working_image);
+ //--------------------------------------------------------------------
+ //--------------------------------------------------------------------
// Step 7: autocrop
- typedef clitk::AutoCropFilter<ImageType> CropFilterType;
- typename CropFilterType::Pointer cropFilter = CropFilterType::New();
- cropFilter->SetInput(working_image);
- cropFilter->ReleaseDataFlagOff();
- cropFilter->Update();
- working_image = cropFilter->GetOutput();
- this->template StopCurrentStep<ImageType>(working_image);
+ if (GetAutoCropFlag()) {
+ StartNewStep("Final AutoCrop");
+ typedef clitk::AutoCropFilter<ImageType> CropFilterType;
+ typename CropFilterType::Pointer cropFilter = CropFilterType::New();
+ cropFilter->SetInput(working_image);
+ cropFilter->ReleaseDataFlagOff();
+ cropFilter->Update();
+ working_image = cropFilter->GetOutput();
+ StopCurrentStep<ImageType>(working_image);
+ }
//--------------------------------------------------------------------
//--------------------------------------------------------------------
// Final Step -> set output
this->SetNthOutput(0, working_image);
- return;
+ // this->GraftOutput(working_image);
}
//--------------------------------------------------------------------
git://git.creatis.insa-lyon.fr / clitk.git / blobdiff