X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=blobdiff_plain;f=itk%2FclitkSliceBySliceRelativePositionFilter.txx;h=16571ea7a224df6b17d87425d7cf8bde619a66b0;hb=f05cc2cf37d7ec960b73fa9d2393d5a888e1f87a;hp=d75cc889be2fb6a39b06fceb688e69c438425b65;hpb=32539f6a33f6a97a8493a1fa11360d09be843a4c;p=clitk.git diff --git a/itk/clitkSliceBySliceRelativePositionFilter.txx b/itk/clitkSliceBySliceRelativePositionFilter.txx index d75cc88..16571ea 100644 --- a/itk/clitkSliceBySliceRelativePositionFilter.txx +++ b/itk/clitkSliceBySliceRelativePositionFilter.txx @@ -17,6 +17,7 @@ ======================================================================-====*/ // clitk +#include "clitkCropLikeImageFilter.h" #include "clitkSegmentationUtils.h" #include "clitkExtractSliceFilter.h" #include "clitkResampleImageWithOptionsFilter.h" @@ -41,6 +42,8 @@ SliceBySliceRelativePositionFilter(): SetObjectCCLSelectionDimension(0); SetObjectCCLSelectionDirection(1); ObjectCCLSelectionIgnoreSingleCCLFlagOff(); + VerboseSlicesFlagOff(); + this->SetK1(vcl_acos(-1.0)/2); } //-------------------------------------------------------------------- @@ -73,23 +76,31 @@ SetInputObject(const ImageType * image) template void clitk::SliceBySliceRelativePositionFilter:: -PrintOptions() +PrintOptions(std::ostream & os) { - DD(this->GetDirection()); - DD((int)this->GetObjectBackgroundValue()); - DDV(this->GetOrientationTypeString(), (uint)this->GetNumberOfAngles()); - DD(this->GetIntermediateSpacingFlag()); - DD(this->GetIntermediateSpacing()); - DD(this->GetFuzzyThreshold()); - DD(this->GetUniqueConnectedComponentBySliceFlag()); - DD(this->GetAutoCropFlag()); - DD(this->GetInverseOrientationFlag()); - DD(this->GetRemoveObjectFlag()); - DD(this->GetCombineWithOrFlag()); - DD(this->GetUseTheLargestObjectCCLFlag()); - DD(this->GetObjectCCLSelectionFlag()); - DD(this->GetObjectCCLSelectionDimension()); - DD(this->GetObjectCCLSelectionIgnoreSingleCCLFlag()); + os << "Slice direction = " << this->GetDirection() << std::endl + << "BG value = " << (int)this->GetBackgroundValue() << std::endl; + for(int i=0; iGetNumberOfAngles(); i++) { + os << "Orientation = " << this->GetOrientationTypeString()[i] << std::endl; + os << "Angles = " << clitk::rad2deg(this->GetAngle1InRad(i)) + << " " << clitk::rad2deg(this->GetAngle2InRad(i)) << std::endl; + } + os << "InverseOrientationFlag = " << this->GetInverseOrientationFlag() << std::endl + << "SpacingFlag = " << this->GetIntermediateSpacingFlag() << std::endl + << "Spacing = " << this->GetIntermediateSpacing() << std::endl + << "FuzzyThreshold = " << this->GetFuzzyThreshold() << std::endl + << "UniqueConnectedComponentBySliceFlag = " << this->GetUniqueConnectedComponentBySliceFlag() << std::endl + << "AutoCropFlag = " << this->GetAutoCropFlag() << std::endl + << "RemoveObjectFlag= " << this->GetRemoveObjectFlag() << std::endl + << "CombineWithOrFlag = " << this->GetCombineWithOrFlag() << std::endl + << "UseTheLargestObjectCCLFlag = " << this->GetUseTheLargestObjectCCLFlag() << std::endl + << "ObjectCCLSelectionFlag = " << this->GetObjectCCLSelectionFlag() << std::endl + << "ObjectCCLSelectionDimension = " << this->GetObjectCCLSelectionDimension() << std::endl + << "ObjectCCLSelectionIgnoreSingleCCLFlag = " << this->GetObjectCCLSelectionIgnoreSingleCCLFlag() << std::endl + << "IgnoreEmptySliceObjectFlag = " << this->GetIgnoreEmptySliceObjectFlag() << std::endl + << "(RP) FastFlag = " << this->GetFastFlag() << std::endl + << "(RP) Radius = " << this->GetRadius() << std::endl + << "(RP) K1 = " << this->GetK1() << std::endl; } //-------------------------------------------------------------------- @@ -121,9 +132,13 @@ GenerateOutputInformation() PrintOptions(); } + // if (this->GetFuzzyMapOnlyFlag()) this->ComputeFuzzyMapFlagOn(); + // Get input pointer input = dynamic_cast(itk::ProcessObject::GetInput(0)); object = dynamic_cast(itk::ProcessObject::GetInput(1)); + m_working_object = object; + m_working_input = input; //-------------------------------------------------------------------- // Resample object to the same spacing than input @@ -132,36 +147,67 @@ GenerateOutputInformation() m_working_object = clitk::ResampleImageSpacing(object, input->GetSpacing()); this->template StopCurrentStep(m_working_object); } - else { - m_working_object = object; - } //-------------------------------------------------------------------- - // Pad object to the same size than input + // Resize image according to common area (except in Z) if (!clitk::HaveSameSizeAndSpacing(m_working_object, input)) { + this->StartNewStep("Resize images (union in XY and like input in Z)"); + + /* OLD STUFF this->StartNewStep("Pad object to the same size than input"); m_working_object = clitk::ResizeImageLike(m_working_object, - input, - this->GetObjectBackgroundValue()); + input, + this->GetObjectBackgroundValue()); this->template StopCurrentStep(m_working_object); - } - else { - } + */ + + // Compute union of bounding boxes in X and Y + static const unsigned int dim = ImageType::ImageDimension; + typedef itk::BoundingBox BBType; + typename BBType::Pointer bb1 = BBType::New(); + ComputeBBFromImageRegion(m_working_object, m_working_object->GetLargestPossibleRegion(), bb1); + typename BBType::Pointer bb2 = BBType::New(); + ComputeBBFromImageRegion(input, input->GetLargestPossibleRegion(), bb2); + typename BBType::Pointer bbo = BBType::New(); + ComputeBBUnion(bbo, bb1, bb2); + + //We set Z BB like input + typename ImageType::PointType maxs = bbo->GetMaximum(); + typename ImageType::PointType mins = bbo->GetMinimum(); + maxs[2] = bb2->GetMaximum()[2]; + mins[2] = bb2->GetMinimum()[2]; + bbo->SetMaximum(maxs); + bbo->SetMinimum(mins); + + // Crop + m_working_input = clitk::ResizeImageLike(input, bbo, this->GetBackgroundValue()); + m_working_object = clitk::ResizeImageLike(m_working_object, + m_working_input, + this->GetObjectBackgroundValue()); + + // Index can be negative in some cases, and lead to problem with + // some filter. So we correct it. + m_working_input = clitk::RemoveNegativeIndexFromRegion(m_working_input); + m_working_object = clitk::RemoveNegativeIndexFromRegion(m_working_object); - /* + // End + this->template StopCurrentStep(m_working_input); + } + + //-------------------------------------------------------------------- + /* Steps : - extract vector of slices in input, in object - slice by slice rel position - joint result - post process */ - //-------------------------------------------------------------------- // Extract input slices this->StartNewStep("Extract input slices"); typedef clitk::ExtractSliceFilter ExtractSliceFilterType; typename ExtractSliceFilterType::Pointer extractSliceFilter = ExtractSliceFilterType::New(); - extractSliceFilter->SetInput(input); + extractSliceFilter->SetInput(m_working_input); extractSliceFilter->SetDirection(GetDirection()); extractSliceFilter->Update(); typedef typename ExtractSliceFilterType::SliceType SliceType; @@ -173,110 +219,161 @@ GenerateOutputInformation() // Extract object slices this->StartNewStep("Extract object slices"); extractSliceFilter = ExtractSliceFilterType::New(); - extractSliceFilter->SetInput(m_working_object);//object); + extractSliceFilter->SetInput(m_working_object); extractSliceFilter->SetDirection(GetDirection()); extractSliceFilter->Update(); std::vector mObjectSlices; extractSliceFilter->GetOutputSlices(mObjectSlices); this->template StopCurrentStep(mObjectSlices[0]); + //-------------------------------------------------------------------- + // Prepare fuzzy slices (if needed) + std::vector mFuzzyMapSlices; + mFuzzyMapSlices.resize(mInputSlices.size()); + //-------------------------------------------------------------------- // Perform slice by slice relative position - this->StartNewStep("Perform slice by slice relative position"); + this->StartNewStep("Perform slice by slice relative position ("+toString(mInputSlices.size())+")"); for(unsigned int i=0; i(mObjectSlices[i], 0, true, 1, nb); - if ((!GetIgnoreEmptySliceObjectFlag()) || (nb!=0)) { - // Select or not a single CCL ? - if (GetUseTheLargestObjectCCLFlag()) { - mObjectSlices[i] = KeepLabels(mObjectSlices[i], 0, 1, 1, 1, true); - } + if (GetVerboseSlicesFlag()) { + std::cout << "slice " << i << " nb = " << nb << std::endl; + } - // Select a single according to a position if more than one CCL - if (GetObjectCCLSelectionFlag()) { - // if several CCL, choose the most extrema according a direction, - // if not -> should we consider this slice ? - if (nb<2) { - if (GetObjectCCLSelectionIgnoreSingleCCLFlag()) { - mObjectSlices[i] = SetBackground(mObjectSlices[i], mObjectSlices[i], - 1, this->GetBackgroundValue(), - true); - } + // If no object and empty slices and if we need the full fuzzy map, create a dummy one. + if ((nb==0) && (this->GetFuzzyMapOnlyFlag())) { + typename FloatSliceType::Pointer one = FloatSliceType::New(); + one->CopyInformation(mObjectSlices[0]); + one->SetRegions(mObjectSlices[0]->GetLargestPossibleRegion()); + one->Allocate(); + one->FillBuffer(2.0); + mFuzzyMapSlices[i] = one; + } // End nb==0 && GetComputeFuzzyMapFlag + else { + if ((!GetIgnoreEmptySliceObjectFlag()) || (nb!=0)) { + + // Select or not a single CCL ? + if (GetUseTheLargestObjectCCLFlag()) { + mObjectSlices[i] = KeepLabels(mObjectSlices[i], 0, 1, 1, 1, true); } - int dim = GetObjectCCLSelectionDimension(); - int direction = GetObjectCCLSelectionDirection(); - std::vector centroids; - ComputeCentroids(mObjectSlices[i], this->GetBackgroundValue(), centroids); - uint index=1; - for(uint j=1; j centroids[index][dim]) index = j; + + // Select a single according to a position if more than one CCL + if (GetObjectCCLSelectionFlag()) { + // if several CCL, choose the most extrema according a direction, + // if not -> should we consider this slice ? + if (nb<2) { + if (GetObjectCCLSelectionIgnoreSingleCCLFlag()) { + mObjectSlices[i] = SetBackground(mObjectSlices[i], mObjectSlices[i], + 1, this->GetBackgroundValue(), + true); + } } - else { - if (centroids[j][dim] < centroids[index][dim]) index = j; + int dim = GetObjectCCLSelectionDimension(); + int direction = GetObjectCCLSelectionDirection(); + std::vector centroids; + ComputeCentroids(mObjectSlices[i], this->GetBackgroundValue(), centroids); + uint index=1; + for(uint j=1; j centroids[index][dim]) index = j; + } + else { + if (centroids[j][dim] < centroids[index][dim]) index = j; + } } - } - for(uint v=1; v(mObjectSlices[i], mObjectSlices[i], - (char)v, this->GetBackgroundValue(), - true); + for(uint v=1; v(mObjectSlices[i], mObjectSlices[i], + (char)v, this->GetBackgroundValue(), + true); + } } + } // end GetbjectCCLSelectionFlag = true + + // Relative position + typedef clitk::AddRelativePositionConstraintToLabelImageFilter RelPosFilterType; + typename RelPosFilterType::Pointer relPosFilter = RelPosFilterType::New(); + relPosFilter->VerboseStepFlagOff(); + if (GetVerboseSlicesFlag()) { + std::cout << "Slice " << i << std::endl; + relPosFilter->VerboseStepFlagOn(); + //relPosFilter->WriteStepFlagOn(); } - } // end GetbjectCCLSelectionFlag = true - - // Relative position - typedef clitk::AddRelativePositionConstraintToLabelImageFilter RelPosFilterType; - typename RelPosFilterType::Pointer relPosFilter = RelPosFilterType::New(); - - relPosFilter->VerboseStepFlagOff(); - relPosFilter->WriteStepFlagOff(); - relPosFilter->SetBackgroundValue(this->GetBackgroundValue()); - relPosFilter->SetInput(mInputSlices[i]); - relPosFilter->SetInputObject(mObjectSlices[i]); - relPosFilter->SetRemoveObjectFlag(this->GetRemoveObjectFlag()); - // This flag (InverseOrientation) *must* be set before - // AddOrientation because AddOrientation can change it. - relPosFilter->SetInverseOrientationFlag(this->GetInverseOrientationFlag()); - for(int j=0; jGetNumberOfAngles(); j++) { - relPosFilter->AddOrientationTypeString(this->GetOrientationTypeString(j)); - //DD(this->GetOrientationTypeString(j)); - } - //DD(this->GetInverseOrientationFlag()); - //relPosFilter->SetOrientationType(this->GetOrientationType()); - relPosFilter->SetIntermediateSpacing(this->GetIntermediateSpacing()); - relPosFilter->SetIntermediateSpacingFlag(this->GetIntermediateSpacingFlag()); - relPosFilter->SetFuzzyThreshold(this->GetFuzzyThreshold()); - relPosFilter->AutoCropFlagOff(); // important ! because we join the slices after this loop - relPosFilter->SetCombineWithOrFlag(this->GetCombineWithOrFlag()); - relPosFilter->Update(); - mInputSlices[i] = relPosFilter->GetOutput(); - - // Select main CC if needed - if (GetUniqueConnectedComponentBySliceFlag()) { - mInputSlices[i] = Labelize(mInputSlices[i], 0, true, 1); - mInputSlices[i] = KeepLabels(mInputSlices[i], 0, 1, 1, 1, true); + relPosFilter->WriteStepFlagOff(); + // relPosFilter->VerboseMemoryFlagOn(); + relPosFilter->SetCurrentStepBaseId(this->GetCurrentStepId()+"-"+toString(i)); + relPosFilter->SetBackgroundValue(this->GetBackgroundValue()); + relPosFilter->SetInput(mInputSlices[i]); + relPosFilter->SetInputObject(mObjectSlices[i]); + relPosFilter->SetRemoveObjectFlag(this->GetRemoveObjectFlag()); + // This flag (InverseOrientation) *must* be set before + // AddOrientation because AddOrientation can change it. + relPosFilter->SetInverseOrientationFlag(this->GetInverseOrientationFlag()); + for(int j=0; jGetNumberOfAngles(); j++) { + relPosFilter->AddAnglesInRad(this->GetAngle1InRad(j), this->GetAngle2InRad(j)); + } + relPosFilter->SetIntermediateSpacing(this->GetIntermediateSpacing()); + relPosFilter->SetIntermediateSpacingFlag(this->GetIntermediateSpacingFlag()); + relPosFilter->SetFuzzyThreshold(this->GetFuzzyThreshold()); + relPosFilter->AutoCropFlagOff(); // important ! because we join the slices after this loop + relPosFilter->SetCombineWithOrFlag(this->GetCombineWithOrFlag()); + // should we stop after fuzzy map ? + relPosFilter->SetFuzzyMapOnlyFlag(this->GetFuzzyMapOnlyFlag()); + // relPosFilter->SetComputeFuzzyMapFlag(this->GetComputeFuzzyMapFlag()); + relPosFilter->SetFastFlag(this->GetFastFlag()); + relPosFilter->SetRadius(this->GetRadius()); + relPosFilter->SetK1(this->GetK1()); + + // Go ! + relPosFilter->Update(); + + // If we stop after the fuzzy map, store the fuzzy slices + if (this->GetFuzzyMapOnlyFlag()) { + mFuzzyMapSlices[i] = relPosFilter->GetFuzzyMap(); + // writeImage(mFuzzyMapSlices[i], "slice_"+toString(i)+".mha"); + } + + // Set input slices + if (!this->GetFuzzyMapOnlyFlag()) { + mInputSlices[i] = relPosFilter->GetOutput(); + // Select main CC if needed + if (GetUniqueConnectedComponentBySliceFlag()) { + mInputSlices[i] = Labelize(mInputSlices[i], 0, true, 1); + mInputSlices[i] = KeepLabels(mInputSlices[i], 0, 1, 1, 1, true); + } + } + } /* // Select unique CC according to the most in a given direction if (GetUniqueConnectedComponentBySliceAccordingToADirection()) { - int nb; - mInputSlices[i] = LabelizeAndCountNumberOfObjects(mInputSlices[i], 0, true, 1, nb); - std::vector & centroids; - ComputeCentroids - } + int nb; + mInputSlices[i] = LabelizeAndCountNumberOfObjects(mInputSlices[i], 0, true, 1, nb); + std::vector & centroids; + ComputeCentroids + } */ - } - } + + } // End nb!=0 || GetComputeFuzzyMapFlagOFF + + } // end for i mInputSlices // Join the slices - m_working_input = clitk::JoinSlices(mInputSlices, input, GetDirection()); + m_working_input = clitk::JoinSlices(mInputSlices, m_working_input, GetDirection()); this->template StopCurrentStep(m_working_input); + // Join the fuzzy map if needed + if (this->GetFuzzyMapOnlyFlag()) { + this->m_FuzzyMap = clitk::JoinSlices(mFuzzyMapSlices, m_working_input, GetDirection()); + this->template StopCurrentStep(this->m_FuzzyMap); + if (this->GetFuzzyMapOnlyFlag()) return; + } + //-------------------------------------------------------------------- // Step 7: autocrop if (this->GetAutoCropFlag()) { @@ -307,6 +404,7 @@ GenerateData() //-------------------------------------------------------------------- // Final Step -> set output //this->SetNthOutput(0, m_working_input); + if (this->GetFuzzyMapOnlyFlag()) return; // no output in this case this->GraftOutput(m_working_input); return; }