X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=blobdiff_plain;f=itk%2FclitkAddRelativePositionConstraintToLabelImageFilter.txx;h=8e12c225e59852544637c5db492d6b26afb532df;hb=72375037f90c596a034b2ebe5e54e209e7b45511;hp=83993831fa4f2f86606b89025378ea27b1bed884;hpb=cd0efd1d1f63d972c449b0550ab0b3c8ed500a15;p=clitk.git 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; } //--------------------------------------------------------------------