/*========================================================================= Program: vv http://www.creatis.insa-lyon.fr/rio/vv Authors belong to: - University of LYON http://www.universite-lyon.fr/ - Léon Bérard cancer center http://oncora1.lyon.fnclcc.fr - CREATIS CNRS laboratory http://www.creatis.insa-lyon.fr This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the copyright notices for more information. It is distributed under dual licence - BSD See included LICENSE.txt file - CeCILL-B http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html ======================================================================-====*/ // clitk #include "clitkCropLikeImageFilter.h" #include "clitkSegmentationUtils.h" #include "clitkExtractSliceFilter.h" #include "clitkResampleImageWithOptionsFilter.h" // itk #include //-------------------------------------------------------------------- template clitk::SliceBySliceRelativePositionFilter:: SliceBySliceRelativePositionFilter(): clitk::AddRelativePositionConstraintToLabelImageFilter() { SetDirection(2); UniqueConnectedComponentBySliceFlagOff(); SetIgnoreEmptySliceObjectFlag(false); UseTheLargestObjectCCLFlagOff(); this->VerboseStepFlagOff(); this->WriteStepFlagOff(); this->SetCombineWithOrFlag(false); ObjectCCLSelectionFlagOff(); SetObjectCCLSelectionDimension(0); SetObjectCCLSelectionDirection(1); ObjectCCLSelectionIgnoreSingleCCLFlagOff(); } //-------------------------------------------------------------------- //-------------------------------------------------------------------- template void clitk::SliceBySliceRelativePositionFilter:: SetInput(const ImageType * image) { // Process object is not const-correct so the const casting is required. this->SetNthInput(0, const_cast(image)); } //-------------------------------------------------------------------- //-------------------------------------------------------------------- template void clitk::SliceBySliceRelativePositionFilter:: SetInputObject(const ImageType * image) { // Process object is not const-correct so the const casting is required. this->SetNthInput(1, const_cast(image)); } //-------------------------------------------------------------------- //-------------------------------------------------------------------- template void clitk::SliceBySliceRelativePositionFilter:: PrintOptions(std::ostream & os) { os << "Slice direction = " << this->GetDirection() << std::endl << "BG value = " << 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; } //-------------------------------------------------------------------- //-------------------------------------------------------------------- template void clitk::SliceBySliceRelativePositionFilter:: GenerateInputRequestedRegion() { // Call default itk::ImageToImageFilter::GenerateInputRequestedRegion(); // Get input pointers and set requested region to common region ImagePointer input1 = dynamic_cast(itk::ProcessObject::GetInput(0)); ImagePointer input2 = dynamic_cast(itk::ProcessObject::GetInput(1)); input1->SetRequestedRegion(input1->GetLargestPossibleRegion()); input2->SetRequestedRegion(input2->GetLargestPossibleRegion()); } //-------------------------------------------------------------------- //-------------------------------------------------------------------- template void clitk::SliceBySliceRelativePositionFilter:: GenerateOutputInformation() { if (this->GetVerboseOptionFlag()) { PrintOptions(); } // Get input pointer input = dynamic_cast(itk::ProcessObject::GetInput(0)); object = dynamic_cast(itk::ProcessObject::GetInput(1)); //-------------------------------------------------------------------- // Resample object to the same spacing than input if (!clitk::HaveSameSpacing(object, input)) { this->StartNewStep("Resample object to the same spacing than input"); 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 if (!clitk::HaveSameSizeAndSpacing(m_working_object, input)) { 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); } else { } /* - 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->SetDirection(GetDirection()); extractSliceFilter->Update(); typedef typename ExtractSliceFilterType::SliceType SliceType; std::vector mInputSlices; extractSliceFilter->GetOutputSlices(mInputSlices); this->template StopCurrentStep(mInputSlices[0]); //-------------------------------------------------------------------- // Extract object slices this->StartNewStep("Extract object slices"); extractSliceFilter = ExtractSliceFilterType::New(); extractSliceFilter->SetInput(m_working_object);//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"); for(unsigned int i=0; i(mObjectSlices[i], 0, true, 1, nb); // If no object and empty slices : 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; } 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(); 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)); relPosFilter->AddAnglesInRad(this->GetAngle1InRad(j), this->GetAngle2InRad(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()); // should we stop after fuzzy map ? relPosFilter->SetFuzzyMapOnlyFlag(this->GetFuzzyMapOnlyFlag()); // 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"); } else { 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 } */ } } // Join the fuzzy map if needed if (this->GetFuzzyMapOnlyFlag()) { this->m_FuzzyMap = clitk::JoinSlices(mFuzzyMapSlices, input, GetDirection()); this->template StopCurrentStep(this->m_FuzzyMap); return; } // Join the slices m_working_input = clitk::JoinSlices(mInputSlices, input, GetDirection()); this->template StopCurrentStep(m_working_input); //-------------------------------------------------------------------- // 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); } // 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); if (this->GetFuzzyMapOnlyFlag()) return; // no output in this case this->GraftOutput(m_working_input); return; } //--------------------------------------------------------------------