From a24b0a699298efe54b53c53cb215455fecd633fe Mon Sep 17 00:00:00 2001 From: dsarrut Date: Tue, 15 Feb 2011 10:51:03 +0000 Subject: [PATCH] correction orientation, add 'keep object' option --- ...tivePositionConstraintToLabelImageFilter.h | 42 +-- ...vePositionConstraintToLabelImageFilter.txx | 247 ++++++++++++------ 2 files changed, 193 insertions(+), 96 deletions(-) diff --git a/itk/clitkAddRelativePositionConstraintToLabelImageFilter.h b/itk/clitkAddRelativePositionConstraintToLabelImageFilter.h index 3a82ee2..c365571 100644 --- a/itk/clitkAddRelativePositionConstraintToLabelImageFilter.h +++ b/itk/clitkAddRelativePositionConstraintToLabelImageFilter.h @@ -88,15 +88,13 @@ namespace clitk { void SetInputObject(const ImageType * image); // Options - void SetOrientationType(OrientationTypeEnumeration orientation); - itkGetConstMacro(OrientationType, OrientationTypeEnumeration); - void SetOrientationTypeString(std::string s); - itkGetConstMacro(OrientationTypeString, std::string); - - void SetAngle1(double a); - void SetAngle2(double a); - itkGetConstMacro(Angle1, double); - itkGetConstMacro(Angle2, double); + void AddOrientationType(OrientationTypeEnumeration orientation); + void AddOrientationTypeString(std::string s); + void ClearOrientationType(); + void AddAngles(double a, double b); + int GetNumberOfAngles(); + std::string GetOrientationTypeString(int i) { return m_OrientationTypeString[i]; } + std::vector & GetOrientationTypeString() { return m_OrientationTypeString; } itkGetConstMacro(ResampleBeforeRelativePositionFilter, bool); itkSetMacro(ResampleBeforeRelativePositionFilter, bool); @@ -118,25 +116,35 @@ namespace clitk { itkSetMacro(AutoCropFlag, bool); itkBooleanMacro(AutoCropFlag); - itkGetConstMacro(NotFlag, bool); - itkSetMacro(NotFlag, bool); - itkBooleanMacro(NotFlag); + itkGetConstMacro(InverseOrientationFlag, bool); + itkSetMacro(InverseOrientationFlag, bool); + itkBooleanMacro(InverseOrientationFlag); + + itkGetConstMacro(RemoveObjectFlag, bool); + itkSetMacro(RemoveObjectFlag, bool); + itkBooleanMacro(RemoveObjectFlag); + + itkGetConstMacro(CombineWithOrFlag, bool); + itkSetMacro(CombineWithOrFlag, bool); + itkBooleanMacro(CombineWithOrFlag); protected: AddRelativePositionConstraintToLabelImageFilter(); virtual ~AddRelativePositionConstraintToLabelImageFilter() {} - OrientationTypeEnumeration m_OrientationType; - std::string m_OrientationTypeString; + std::vector m_OrientationType; + std::vector m_OrientationTypeString; double m_IntermediateSpacing; double m_FuzzyThreshold; PixelType m_BackgroundValue; PixelType m_ObjectBackgroundValue; - double m_Angle1; - double m_Angle2; + std::vector m_Angle1; + std::vector m_Angle2; bool m_ResampleBeforeRelativePositionFilter; bool m_AutoCropFlag; - bool m_NotFlag; + bool m_InverseOrientationFlag; + bool m_RemoveObjectFlag; + bool m_CombineWithOrFlag; virtual void GenerateOutputInformation(); virtual void GenerateInputRequestedRegion(); diff --git a/itk/clitkAddRelativePositionConstraintToLabelImageFilter.txx b/itk/clitkAddRelativePositionConstraintToLabelImageFilter.txx index 8399383..8e12c22 100644 --- a/itk/clitkAddRelativePositionConstraintToLabelImageFilter.txx +++ b/itk/clitkAddRelativePositionConstraintToLabelImageFilter.txx @@ -25,12 +25,14 @@ // itk #include -#include "itkStatisticsLabelMapFilter.h" -#include "itkLabelImageToStatisticsLabelMapFilter.h" -#include "itkRegionOfInterestImageFilter.h" -#include "itkBinaryThresholdImageFilter.h" -#include "itkBinaryErodeImageFilter.h" -#include "itkBinaryBallStructuringElement.h" +#include +#include +#include +#include +#include +#include +#include +#include // itk [Bloch et al] #include "RelativePositionPropImageFilter.h" @@ -46,11 +48,13 @@ AddRelativePositionConstraintToLabelImageFilter(): SetFuzzyThreshold(0.6); SetBackgroundValue(0); SetObjectBackgroundValue(0); - SetOrientationType(LeftTo); + ClearOrientationType(); ResampleBeforeRelativePositionFilterOn(); SetIntermediateSpacing(10); AutoCropFlagOn(); - NotFlagOff(); + InverseOrientationFlagOff(); + RemoveObjectFlagOn(); + CombineWithOrFlagOff(); } //-------------------------------------------------------------------- @@ -83,15 +87,43 @@ SetInputObject(const ImageType * image) template void clitk::AddRelativePositionConstraintToLabelImageFilter:: -SetOrientationTypeString(std::string t) +ClearOrientationType() { - SetOrientationType(Angle); - if (t[0] == 'L') SetOrientationType(LeftTo); - if (t[0] == 'R') SetOrientationType(RightTo); - if (t[0] == 'A') SetOrientationType(AntTo); - if (t[0] == 'P') SetOrientationType(PostTo); - if (t[0] == 'S') SetOrientationType(SupTo); - if (t[0] == 'I') SetOrientationType(InfTo); + m_OrientationTypeString.clear(); + m_OrientationType.clear(); + m_Angle1.clear(); + m_Angle2.clear(); +} +//-------------------------------------------------------------------- + + +//-------------------------------------------------------------------- +template +int +clitk::AddRelativePositionConstraintToLabelImageFilter:: +GetNumberOfAngles() +{ + return m_OrientationType.size(); +} +//-------------------------------------------------------------------- + + +//-------------------------------------------------------------------- +template +void +clitk::AddRelativePositionConstraintToLabelImageFilter:: +AddOrientationTypeString(std::string t) +{ + m_OrientationTypeString.push_back(t); + switch (t[0]) { + case 'L' : AddOrientationType(LeftTo); break; + case 'R' : AddOrientationType(RightTo);break; + case 'A' : AddOrientationType(AntTo);break; + case 'P' : AddOrientationType(PostTo);break; + case 'S' : AddOrientationType(SupTo);break; + case 'I' : AddOrientationType(InfTo);break; + default: clitkExceptionMacro("Error, you must provide L,R or A,P or S,I"); + } } //-------------------------------------------------------------------- @@ -130,10 +162,11 @@ GenerateInputRequestedRegion() template void clitk::AddRelativePositionConstraintToLabelImageFilter:: -SetAngle1(double a) +AddAngles(double a, double b) { - SetOrientationType(Angle); - m_Angle1 = a; + AddOrientationTypeString("Angle"); + m_Angle1.push_back(a); + m_Angle2.push_back(b); } //-------------------------------------------------------------------- @@ -142,29 +175,35 @@ SetAngle1(double a) template void clitk::AddRelativePositionConstraintToLabelImageFilter:: -SetAngle2(double a) +AddOrientationType(OrientationTypeEnumeration orientation) { - 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; + m_OrientationType.push_back(orientation); + switch (orientation) { + case LeftTo: + m_Angle1.push_back(clitk::deg2rad(0)); + m_Angle2.push_back(clitk::deg2rad(0)); + break; + case RightTo: + 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 @@ -184,17 +223,9 @@ void clitk::AddRelativePositionConstraintToLabelImageFilter:: GenerateData() { - // Print Option - /* - DD(GetFuzzyThreshold()); - DD((int)GetBackgroundValue()); - DD((int)GetObjectBackgroundValue()); - DD(GetOrientationType()); - DD(GetResampleBeforeRelativePositionFilter()); - DD(GetIntermediateSpacing()); - DD(GetAutoCropFlag()); - DD(GetNotFlag()); - */ + if (GetNumberOfAngles() <1) { + clitkExceptionMacro("Add at least one orientation type"); + } // Get input pointer input = dynamic_cast(itk::ProcessObject::GetInput(0)); @@ -229,7 +260,8 @@ GenerateData() typename ImageType::Pointer output = ImageType::New(); SizeType size; for(unsigned int i=0; iGetLargestPossibleRegion().GetSize()[i]*input->GetSpacing()[i])/(double)working_image->GetSpacing()[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 @@ -259,7 +291,7 @@ GenerateData() // typename PadFilterType::InputImageIndexType index; for(unsigned int i=0; iGetSpacing()[i]/(double)working_image->GetSpacing()[i] - + lrint((working_image->GetOrigin()[i] - input->GetOrigin()[i])/working_image->GetSpacing()[i]); + + lrint((working_image->GetOrigin()[i] - input->GetOrigin()[i])/working_image->GetSpacing()[i]); } padFilter->SetSourceImage(working_image); padFilter->SetDestinationImage(output); @@ -274,21 +306,56 @@ GenerateData() } // Keep object image (with resampline and pad) object_resampled = working_image; - // StopCurrentStep(working_image); + // StopCurrentStep(working_image); // Step 3: compute rel pos in object StartNewStep("Relative Position Map"); 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(); + typename RelPosFilterType::Pointer relPosFilter; + + typename FloatImageType::Pointer m_FuzzyMap; + for(int i=0; iSetInput(working_image); + relPosFilter->SetAlpha1(m_Angle1[i]); // xy plane + relPosFilter->SetAlpha2(m_Angle2[i]); + relPosFilter->SetK1(M_PI/2.0); // Opening parameter, default = pi/2 + relPosFilter->SetFast(true); + relPosFilter->SetRadius(1); // seems sufficient in this case + // relPosFilter->SetVerboseProgress(true); + relPosFilter->Update(); + relPos = relPosFilter->GetOutput(); + + if (GetNumberOfAngles() != 1) { + // Creation of the first m_FuzzyMap + if (i==0) { + m_FuzzyMap = clitk::NewImageLike(relPos, true); + m_FuzzyMap->FillBuffer(0.0); + } + + // Add to current fuzzy map + typedef itk::AddImageFilter 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::DivideByConstantImageFilter DivideFilter; + typename DivideFilter::Pointer divideFilter = DivideFilter::New(); + divideFilter->SetInput(m_FuzzyMap); + divideFilter->SetConstant(GetNumberOfAngles()); + divideFilter->Update(); + m_FuzzyMap = divideFilter->GetOutput(); + } + + relPos = m_FuzzyMap; StopCurrentStep(relPos); //-------------------------------------------------------------------- @@ -344,8 +411,6 @@ GenerateData() //-------------------------------------------------------------------- //-------------------------------------------------------------------- // Pre Step 6: pad if not the same size : it can occur when downsample and upsample - // DD(working_image->GetLargestPossibleRegion()); - // DD(input->GetLargestPossibleRegion()); //if (!HaveSameSizeAndSpacing(working_image, input)) { if (working_image->GetLargestPossibleRegion() != input->GetLargestPossibleRegion()) { StartNewStep("Pad to get the same size than input"); @@ -380,22 +445,47 @@ GenerateData() combineFilter->SetForegroundValue(m_BackgroundValue+1); combineFilter->SetInput1(input); combineFilter->SetInput2(working_image); - if (GetNotFlag()) + if (GetInverseOrientationFlag()) combineFilter->SetOperationType(BoolFilterType::AndNot); - else - combineFilter->SetOperationType(BoolFilterType::And); - combineFilter->InPlaceOn(); + 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(); + // 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(); + } - working_image = combineFilter->GetOutput(); StopCurrentStep(working_image); //-------------------------------------------------------------------- @@ -418,7 +508,6 @@ GenerateData() // Final Step -> set output this->SetNthOutput(0, working_image); // this->GraftOutput(working_image); - return; } //-------------------------------------------------------------------- -- 2.47.1