X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=blobdiff_plain;f=itk%2FclitkSliceBySliceRelativePositionFilter.txx;h=16571ea7a224df6b17d87425d7cf8bde619a66b0;hb=f05cc2cf37d7ec960b73fa9d2393d5a888e1f87a;hp=af9c4c5c35f3378e5d0c6e90c1da399c224acfab;hpb=f34fdef46fa656a1eb32271c0b9bf7d81bc7fc35;p=clitk.git diff --git a/itk/clitkSliceBySliceRelativePositionFilter.txx b/itk/clitkSliceBySliceRelativePositionFilter.txx index af9c4c5..16571ea 100644 --- a/itk/clitkSliceBySliceRelativePositionFilter.txx +++ b/itk/clitkSliceBySliceRelativePositionFilter.txx @@ -17,35 +17,33 @@ ======================================================================-====*/ // clitk +#include "clitkCropLikeImageFilter.h" #include "clitkSegmentationUtils.h" -// #include "clitkBooleanOperatorLabelImageFilter.h" -// #include "clitkAutoCropFilter.h" -// #include "clitkResampleImageWithOptionsFilter.h" -// #include "clitkBooleanOperatorLabelImageFilter.h" #include "clitkExtractSliceFilter.h" +#include "clitkResampleImageWithOptionsFilter.h" -// // itk -// #include -// #include "itkStatisticsLabelMapFilter.h" -// #include "itkLabelImageToStatisticsLabelMapFilter.h" -// #include "itkRegionOfInterestImageFilter.h" -// #include "itkBinaryThresholdImageFilter.h" -// #include "itkBinaryErodeImageFilter.h" -// #include "itkBinaryBallStructuringElement.h" - -// // itk [Bloch et al] -// #include "RelativePositionPropImageFilter.h" +// itk +#include //-------------------------------------------------------------------- template clitk::SliceBySliceRelativePositionFilter:: SliceBySliceRelativePositionFilter(): - clitk::FilterBase(), - itk::ImageToImageFilter() + clitk::AddRelativePositionConstraintToLabelImageFilter() { - this->SetNumberOfRequiredInputs(2); SetDirection(2); - SetObjectBackgroundValue(0); + UniqueConnectedComponentBySliceFlagOff(); + SetIgnoreEmptySliceObjectFlag(false); + UseTheLargestObjectCCLFlagOff(); + this->VerboseStepFlagOff(); + this->WriteStepFlagOff(); + this->SetCombineWithOrFlag(false); + ObjectCCLSelectionFlagOff(); + SetObjectCCLSelectionDimension(0); + SetObjectCCLSelectionDirection(1); + ObjectCCLSelectionIgnoreSingleCCLFlagOff(); + VerboseSlicesFlagOff(); + this->SetK1(vcl_acos(-1.0)/2); } //-------------------------------------------------------------------- @@ -54,7 +52,8 @@ SliceBySliceRelativePositionFilter(): template void clitk::SliceBySliceRelativePositionFilter:: -SetInput(const ImageType * image) { +SetInput(const ImageType * image) +{ // Process object is not const-correct so the const casting is required. this->SetNthInput(0, const_cast(image)); } @@ -65,7 +64,8 @@ SetInput(const ImageType * image) { template void clitk::SliceBySliceRelativePositionFilter:: -SetInputObject(const ImageType * image) { +SetInputObject(const ImageType * image) +{ // Process object is not const-correct so the const casting is required. this->SetNthInput(1, const_cast(image)); } @@ -76,10 +76,31 @@ SetInputObject(const ImageType * image) { template void clitk::SliceBySliceRelativePositionFilter:: -GenerateOutputInformation() { - ImagePointer input = dynamic_cast(itk::ProcessObject::GetInput(0)); - ImagePointer outputImage = this->GetOutput(0); - outputImage->SetRegions(input->GetLargestPossibleRegion()); +PrintOptions(std::ostream & os) +{ + 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; } //-------------------------------------------------------------------- @@ -88,7 +109,8 @@ GenerateOutputInformation() { template void clitk::SliceBySliceRelativePositionFilter:: -GenerateInputRequestedRegion() { +GenerateInputRequestedRegion() +{ // Call default itk::ImageToImageFilter::GenerateInputRequestedRegion(); // Get input pointers and set requested region to common region @@ -99,104 +121,291 @@ GenerateInputRequestedRegion() { } //-------------------------------------------------------------------- - + //-------------------------------------------------------------------- template void clitk::SliceBySliceRelativePositionFilter:: -GenerateData() { +GenerateOutputInformation() +{ + if (this->GetVerboseOptionFlag()) { + 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 if (!clitk::HaveSameSpacing(object, input)) { - StartNewStep("Resample object to the same spacing than input"); + this->StartNewStep("Resample object to the same spacing than input"); m_working_object = clitk::ResampleImageSpacing(object, input->GetSpacing()); - StopCurrentStep(m_working_object); - } - else { - DD("no resampling"); - m_working_object = object; + this->template StopCurrentStep(m_working_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)) { - StartNewStep("Pad object to the same size than input"); - m_working_object = clitk::EnlargeImageLike(m_working_object, - input, - GetObjectBackgroundValue()); - StopCurrentStep(m_working_object); - } - else { - DD("no pad"); + 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()); + this->template StopCurrentStep(m_working_object); + */ + + // 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 - StartNewStep("Extract input slices"); + this->StartNewStep("Extract input slices"); typedef clitk::ExtractSliceFilter ExtractSliceFilterType; - typename ExtractSliceFilterType::Pointer extractSliceFilter = Extractslicefilter::New(); - extractSliceFilter->SetInput(input); + typename ExtractSliceFilterType::Pointer extractSliceFilter = ExtractSliceFilterType::New(); + extractSliceFilter->SetInput(m_working_input); extractSliceFilter->SetDirection(GetDirection()); extractSliceFilter->Update(); typedef typename ExtractSliceFilterType::SliceType SliceType; - std::vector & mInputSlices = extractSliceFilter->GetOutput(); - DD(mInputSlices.size()); - StopCurrentStep(mInputSlices[5]); + std::vector mInputSlices; + extractSliceFilter->GetOutputSlices(mInputSlices); + this->template StopCurrentStep(mInputSlices[0]); //-------------------------------------------------------------------- // Extract object slices - - StartNewStep("Extract object slices"); - extractSliceFilter = Extractslicefilter::New(); - extractSliceFilter->SetInput(input); + this->StartNewStep("Extract object slices"); + extractSliceFilter = ExtractSliceFilterType::New(); + extractSliceFilter->SetInput(m_working_object); extractSliceFilter->SetDirection(GetDirection()); extractSliceFilter->Update(); - std::vector & mObjectSlices = extractSliceFilter->GetOutput(); - DD(mObjectSlices.size()); - StopCurrentStep(mInputSlices[5]); + 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 - StartNewStep("Perform slice by slice relative position"); - for(int i=0; i RelPosFilterType; - typename RelPosFilterType::Pointer relPosFilter = RelPosFilterType::New(); - relPosFilter->VerboseStepOff(); - relPosFilter->WriteStepOff(); - relPosFilter->SetInput(mInputSlices[i]); - relPosFilter->SetInputObject(mObjectSlices[i]); - relPosFilter->SetOrientationType(GetOrientationType()); - relPosFilter->SetIntermediateSpacing(GetIntermediateSpacing()); - relPosFilter->SetFuzzyThreshold(GetFuzzyThreshold()); - relPosFilter->Update(); - mInputSlices[i] = relPosFilter->GetOutput(); + this->StartNewStep("Perform slice by slice relative position ("+toString(mInputSlices.size())+")"); + for(unsigned int i=0; i(mObjectSlices[i], 0, true, 1, nb); + + if (GetVerboseSlicesFlag()) { + std::cout << "slice " << i << " nb = " << nb << std::endl; + } + + // 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); + } + + // 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); + } + } + 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); + } + } + } // 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(); + } + 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 + } + */ + + } // End nb!=0 || GetComputeFuzzyMapFlagOFF + + } // end for i mInputSlices + + // Join the slices + 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; } - typedef itk::JoinSeriesImageFilter JointSeriesFilterType; - typename JointSeriesFilterType::Pointer jointFilter = JointSeriesFilterType::New(); - for(int i=0; iSetInput(i, mInputSlices[i]); + + //-------------------------------------------------------------------- + // Step 7: autocrop + if (this->GetAutoCropFlag()) { + this->StartNewStep("Final AutoCrop"); + typedef clitk::AutoCropFilter CropFilterType; + typename CropFilterType::Pointer cropFilter = CropFilterType::New(); + cropFilter->SetInput(m_working_input); + cropFilter->ReleaseDataFlagOff(); + cropFilter->Update(); + m_working_input = cropFilter->GetOutput(); + this->template StopCurrentStep(m_working_input); } - m_working_input = jointFilter->GetOutput(); - DD(m_working_input->GetSpacing()); - DD(m_working_input->GetOrigin()); - StopCurrentStep(m_working_input); - */ + // Update output info + this->GetOutput(0)->SetRegions(m_working_input->GetLargestPossibleRegion()); +} +//-------------------------------------------------------------------- +//-------------------------------------------------------------------- +template +void +clitk::SliceBySliceRelativePositionFilter:: +GenerateData() +{ + // Get input pointer //-------------------------------------------------------------------- //-------------------------------------------------------------------- // Final Step -> set output - this->SetNthOutput(0, m_working_input); + //this->SetNthOutput(0, m_working_input); + if (this->GetFuzzyMapOnlyFlag()) return; // no output in this case + this->GraftOutput(m_working_input); return; } //--------------------------------------------------------------------