X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=blobdiff_plain;f=itk%2FclitkAddRelativePositionConstraintToLabelImageFilter.txx;h=8e12c225e59852544637c5db492d6b26afb532df;hb=72375037f90c596a034b2ebe5e54e209e7b45511;hp=2b159bf9e9a98874f2a5f167ad41cdd5a27850e9;hpb=302ee878f46ba0f96a4913b1715ae768bdb665b0;p=clitk.git diff --git a/itk/clitkAddRelativePositionConstraintToLabelImageFilter.txx b/itk/clitkAddRelativePositionConstraintToLabelImageFilter.txx index 2b159bf..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,10 +48,13 @@ AddRelativePositionConstraintToLabelImageFilter(): SetFuzzyThreshold(0.6); SetBackgroundValue(0); SetObjectBackgroundValue(0); - SetOrientationType(LeftTo); + ClearOrientationType(); ResampleBeforeRelativePositionFilterOn(); SetIntermediateSpacing(10); - AutoCropOn(); + AutoCropFlagOn(); + InverseOrientationFlagOff(); + RemoveObjectFlagOn(); + CombineWithOrFlagOff(); } //-------------------------------------------------------------------- @@ -78,6 +83,51 @@ SetInputObject(const ImageType * image) //-------------------------------------------------------------------- +//-------------------------------------------------------------------- +template +void +clitk::AddRelativePositionConstraintToLabelImageFilter:: +ClearOrientationType() +{ + 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"); + } +} +//-------------------------------------------------------------------- + + //-------------------------------------------------------------------- template void @@ -112,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); } //-------------------------------------------------------------------- @@ -124,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 @@ -166,6 +223,10 @@ void clitk::AddRelativePositionConstraintToLabelImageFilter:: GenerateData() { + if (GetNumberOfAngles() <1) { + clitkExceptionMacro("Add at least one orientation type"); + } + // Get input pointer input = dynamic_cast(itk::ProcessObject::GetInput(0)); object = dynamic_cast(itk::ProcessObject::GetInput(1)); @@ -199,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 @@ -229,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); @@ -244,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); //-------------------------------------------------------------------- @@ -308,22 +405,24 @@ GenerateData() // resampleFilter->SetVerboseOptions(true); resampleFilter->Update(); working_image = resampleFilter->GetOutput(); - StopCurrentStep(working_image); + StopCurrentStep(working_image); } //-------------------------------------------------------------------- //-------------------------------------------------------------------- // 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 PadFilterType::Pointer padFilter2 = PadFilterType::New(); padFilter2->SetSourceImage(working_image); padFilter2->SetDestinationImage(temp); + // DD(input->GetLargestPossibleRegion().GetIndex()); padFilter2->SetDestinationIndex(input->GetLargestPossibleRegion().GetIndex()); padFilter2->SetSourceRegion(working_image->GetLargestPossibleRegion()); padFilter2->Update(); @@ -340,34 +439,59 @@ GenerateData() StartNewStep("Combine with initial input (boolean And)"); typedef clitk::BooleanOperatorLabelImageFilter BoolFilterType; typename BoolFilterType::Pointer combineFilter = BoolFilterType::New(); - writeImage(input, "i.mhd"); - writeImage(working_image, "w.mhd"); 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(); + 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(); - // writeImage(working_image, "res.mhd"); - - combineFilter = BoolFilterType::New(); - combineFilter->SetInput1(working_image); - combineFilter->SetInput2(object); - combineFilter->SetOperationType(BoolFilterType::AndNot); - combineFilter->InPlaceOn(); - combineFilter->Update(); - working_image = combineFilter->GetOutput(); + // 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(working_image); //-------------------------------------------------------------------- //-------------------------------------------------------------------- // Step 7: autocrop - if (GetAutoCrop()) { + if (GetAutoCropFlag()) { StartNewStep("Final AutoCrop"); typedef clitk::AutoCropFilter CropFilterType; typename CropFilterType::Pointer cropFilter = CropFilterType::New(); @@ -383,7 +507,7 @@ GenerateData() // Final Step -> set output this->SetNthOutput(0, working_image); - return; + // this->GraftOutput(working_image); } //--------------------------------------------------------------------