======================================================================-====*/
// clitk
+#include "clitkCropLikeImageFilter.h"
#include "clitkSegmentationUtils.h"
#include "clitkExtractSliceFilter.h"
#include "clitkResampleImageWithOptionsFilter.h"
template <class ImageType>
clitk::SliceBySliceRelativePositionFilter<ImageType>::
SliceBySliceRelativePositionFilter():
- clitk::FilterBase(),
- itk::ImageToImageFilter<ImageType, ImageType>()
+ clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>()
{
- this->SetNumberOfRequiredInputs(2);
SetDirection(2);
- SetObjectBackgroundValue(0);
- SetFuzzyThreshold(0.6);
- SetOrientationTypeString("Left");
- SetIntermediateSpacing(10);
- ResampleBeforeRelativePositionFilterOff();
- UniqueConnectedComponentBySliceOff();
- NotFlagOff();
+ UniqueConnectedComponentBySliceFlagOff();
+ SetIgnoreEmptySliceObjectFlag(false);
+ UseTheLargestObjectCCLFlagOff();
+ this->VerboseStepFlagOff();
+ this->WriteStepFlagOff();
+ this->SetCombineWithOrFlag(false);
+ ObjectCCLSelectionFlagOff();
+ SetObjectCCLSelectionDimension(0);
+ SetObjectCCLSelectionDirection(1);
+ ObjectCCLSelectionIgnoreSingleCCLFlagOff();
+ VerboseSlicesFlagOff();
+ this->SetK1(std::acos(-1.0)/2);
}
//--------------------------------------------------------------------
//--------------------------------------------------------------------
+//--------------------------------------------------------------------
+template <class ImageType>
+void
+clitk::SliceBySliceRelativePositionFilter<ImageType>::
+PrintOptions(std::ostream & os)
+{
+ os << "Slice direction = " << this->GetDirection() << std::endl
+ << "BG value = " << (int)this->GetBackgroundValue() << std::endl;
+ for(int i=0; i<this->GetNumberOfAngles(); 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;
+}
+//--------------------------------------------------------------------
+
+
//--------------------------------------------------------------------
template <class ImageType>
void
clitk::SliceBySliceRelativePositionFilter<ImageType>::
GenerateOutputInformation()
{
+ if (this->GetVerboseOptionFlag()) {
+ PrintOptions();
+ }
+
+ // if (this->GetFuzzyMapOnlyFlag()) this->ComputeFuzzyMapFlagOn();
+
// Get input pointer
input = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(0));
object = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(1));
+ m_working_object = object;
+ m_working_input = input;
//--------------------------------------------------------------------
// Resample object to the same spacing than input
if (!clitk::HaveSameSpacing<ImageType, ImageType>(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<ImageType>(object, input->GetSpacing());
- StopCurrentStep<ImageType>(m_working_object);
- }
- else {
- m_working_object = object;
+ this->template StopCurrentStep<ImageType>(m_working_object);
}
//--------------------------------------------------------------------
- // Pad object to the same size than input
+ // Resize image according to common area (except in Z)
if (!clitk::HaveSameSizeAndSpacing<ImageType, ImageType>(m_working_object, input)) {
- StartNewStep("Pad object to the same size than 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<ImageType>(m_working_object,
- input,
- GetObjectBackgroundValue());
- StopCurrentStep<ImageType>(m_working_object);
- }
- else {
- }
+ input,
+ this->GetObjectBackgroundValue());
+ this->template StopCurrentStep<ImageType>(m_working_object);
+ */
+
+ // Compute union of bounding boxes in X and Y
+ static const unsigned int dim = ImageType::ImageDimension;
+ typedef itk::BoundingBox<unsigned long, dim> BBType;
+ typename BBType::Pointer bb1 = BBType::New();
+ ComputeBBFromImageRegion<ImageType>(m_working_object, m_working_object->GetLargestPossibleRegion(), bb1);
+ typename BBType::Pointer bb2 = BBType::New();
+ ComputeBBFromImageRegion<ImageType>(input, input->GetLargestPossibleRegion(), bb2);
+ typename BBType::Pointer bbo = BBType::New();
+ ComputeBBUnion<dim>(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<ImageType>(input, bbo, this->GetBackgroundValue());
+ m_working_object = clitk::ResizeImageLike<ImageType>(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<ImageType>(m_working_input);
+ m_working_object = clitk::RemoveNegativeIndexFromRegion<ImageType>(m_working_object);
- /*
+ // End
+ this->template StopCurrentStep<ImageType>(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<ImageType> 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;
std::vector<typename SliceType::Pointer> mInputSlices;
extractSliceFilter->GetOutputSlices(mInputSlices);
- StopCurrentStep<SliceType>(mInputSlices[0]);
+ this->template StopCurrentStep<SliceType>(mInputSlices[0]);
//--------------------------------------------------------------------
// Extract object slices
- StartNewStep("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<typename SliceType::Pointer> mObjectSlices;
extractSliceFilter->GetOutputSlices(mObjectSlices);
- StopCurrentStep<SliceType>(mObjectSlices[0]);
+ this->template StopCurrentStep<SliceType>(mObjectSlices[0]);
+
+ //--------------------------------------------------------------------
+ // Prepare fuzzy slices (if needed)
+ std::vector<typename FloatSliceType::Pointer> mFuzzyMapSlices;
+ mFuzzyMapSlices.resize(mInputSlices.size());
//--------------------------------------------------------------------
// Perform slice by slice relative position
- StartNewStep("Perform slice by slice relative position");
+ this->StartNewStep("Perform slice by slice relative position ("+toString(mInputSlices.size())+")");
for(unsigned int i=0; i<mInputSlices.size(); i++) {
- // Select main CC in each object slice (required ?)
- mObjectSlices[i] = Labelize<SliceType>(mObjectSlices[i], 0, true, 1);
- mObjectSlices[i] = KeepLabels<SliceType>(mObjectSlices[i], 0, 1, 1, 1, true);
-
- // Relative position
- typedef clitk::AddRelativePositionConstraintToLabelImageFilter<SliceType> RelPosFilterType;
- typename RelPosFilterType::Pointer relPosFilter = RelPosFilterType::New();
- relPosFilter->VerboseStepOff();
- relPosFilter->WriteStepOff();
- relPosFilter->SetCurrentStepBaseId(this->GetCurrentStepId());
- relPosFilter->SetInput(mInputSlices[i]);
- relPosFilter->SetInputObject(mObjectSlices[i]);
- relPosFilter->SetNotFlag(GetNotFlag());
- relPosFilter->SetOrientationTypeString(this->GetOrientationTypeString());
- relPosFilter->SetIntermediateSpacing(this->GetIntermediateSpacing());
- relPosFilter->SetResampleBeforeRelativePositionFilter(this->GetResampleBeforeRelativePositionFilter());
- relPosFilter->SetFuzzyThreshold(this->GetFuzzyThreshold());
- relPosFilter->AutoCropFlagOff(); // important ! because we join the slices after this loop
- relPosFilter->Update();
- mInputSlices[i] = relPosFilter->GetOutput();
-
- // Select main CC if needed
- if (GetUniqueConnectedComponentBySlice()) {
- mInputSlices[i] = Labelize<SliceType>(mInputSlices[i], 0, true, 1);
- mInputSlices[i] = KeepLabels<SliceType>(mInputSlices[i], 0, 1, 1, 1, true);
+
+ // Count the number of CCL (allow to ignore empty slice)
+ int nb=0;
+ mObjectSlices[i] = LabelizeAndCountNumberOfObjects<SliceType>(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<SliceType>(mObjectSlices[i], 0, 1, 1, 1, true);
+ }
- typedef itk::JoinSeriesImageFilter<SliceType, ImageType> JoinSeriesFilterType;
- typename JoinSeriesFilterType::Pointer joinFilter = JoinSeriesFilterType::New();
- joinFilter->SetOrigin(input->GetOrigin()[GetDirection()]);
- joinFilter->SetSpacing(input->GetSpacing()[GetDirection()]);
- for(unsigned int i=0; i<mInputSlices.size(); i++) {
- joinFilter->PushBackInput(mInputSlices[i]);
+ // 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<SliceType, SliceType>(mObjectSlices[i], mObjectSlices[i],
+ 1, this->GetBackgroundValue(),
+ true);
+ }
+ }
+ int dim = GetObjectCCLSelectionDimension();
+ int direction = GetObjectCCLSelectionDirection();
+ std::vector<typename SliceType::PointType> centroids;
+ ComputeCentroids<SliceType>(mObjectSlices[i], this->GetBackgroundValue(), centroids);
+ uint index=1;
+ for(uint j=1; j<centroids.size(); j++) {
+ if (direction == 1) {
+ if (centroids[j][dim] > centroids[index][dim]) index = j;
+ }
+ else {
+ if (centroids[j][dim] < centroids[index][dim]) index = j;
+ }
+ }
+ for(uint v=1; v<centroids.size(); v++) {
+ if (v != index) {
+ mObjectSlices[i] = SetBackground<SliceType, SliceType>(mObjectSlices[i], mObjectSlices[i],
+ (char)v, this->GetBackgroundValue(),
+ true);
+ }
+ }
+ } // end GetbjectCCLSelectionFlag = true
+
+ // Relative position
+ typedef clitk::AddRelativePositionConstraintToLabelImageFilter<SliceType> 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; j<this->GetNumberOfAngles(); 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<FloatSliceType>(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<SliceType>(mInputSlices[i], 0, true, 1);
+ mInputSlices[i] = KeepLabels<SliceType>(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<SliceType>(mInputSlices[i], 0, true, 1, nb);
+ std::vector<typename ImageType::PointType> & centroids;
+ ComputeCentroids
+ }
+ */
+
+ } // End nb!=0 || GetComputeFuzzyMapFlagOFF
+
+ } // end for i mInputSlices
+
+ // Join the slices
+ m_working_input = clitk::JoinSlices<ImageType>(mInputSlices, m_working_input, GetDirection());
+ this->template StopCurrentStep<ImageType>(m_working_input);
+
+ // Join the fuzzy map if needed
+ if (this->GetFuzzyMapOnlyFlag()) {
+ this->m_FuzzyMap = clitk::JoinSlices<FloatImageType>(mFuzzyMapSlices, m_working_input, GetDirection());
+ this->template StopCurrentStep<FloatImageType>(this->m_FuzzyMap);
+ if (this->GetFuzzyMapOnlyFlag()) return;
}
- joinFilter->Update();
- m_working_input = joinFilter->GetOutput();
- StopCurrentStep<ImageType>(m_working_input);
//--------------------------------------------------------------------
// Step 7: autocrop
- if (GetAutoCropFlag()) {
- StartNewStep("Final AutoCrop");
+ if (this->GetAutoCropFlag()) {
+ this->StartNewStep("Final AutoCrop");
typedef clitk::AutoCropFilter<ImageType> CropFilterType;
typename CropFilterType::Pointer cropFilter = CropFilterType::New();
cropFilter->SetInput(m_working_input);
cropFilter->ReleaseDataFlagOff();
cropFilter->Update();
m_working_input = cropFilter->GetOutput();
- StopCurrentStep<ImageType>(m_working_input);
+ this->template StopCurrentStep<ImageType>(m_working_input);
}
// Update output info
//--------------------------------------------------------------------
// 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;
}
git://git.creatis.insa-lyon.fr / clitk.git / blobdiff