1 /*=========================================================================
2 Program: vv http://www.creatis.insa-lyon.fr/rio/vv
5 - University of LYON http://www.universite-lyon.fr/
6 - Léon Bérard cancer center http://oncora1.lyon.fnclcc.fr
7 - CREATIS CNRS laboratory http://www.creatis.insa-lyon.fr
9 This software is distributed WITHOUT ANY WARRANTY; without even
10 the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
11 PURPOSE. See the copyright notices for more information.
13 It is distributed under dual licence
15 - BSD See included LICENSE.txt file
16 - CeCILL-B http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
17 ======================================================================-====*/
20 #include "clitkCropLikeImageFilter.h"
21 #include "clitkSegmentationUtils.h"
22 #include "clitkExtractSliceFilter.h"
23 #include "clitkResampleImageWithOptionsFilter.h"
26 #include <itkJoinSeriesImageFilter.h>
28 //--------------------------------------------------------------------
29 template <class ImageType>
30 clitk::SliceBySliceRelativePositionFilter<ImageType>::
31 SliceBySliceRelativePositionFilter():
32 clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>()
35 UniqueConnectedComponentBySliceFlagOff();
36 SetIgnoreEmptySliceObjectFlag(false);
37 UseTheLargestObjectCCLFlagOff();
38 this->VerboseStepFlagOff();
39 this->WriteStepFlagOff();
40 this->SetCombineWithOrFlag(false);
41 ObjectCCLSelectionFlagOff();
42 SetObjectCCLSelectionDimension(0);
43 SetObjectCCLSelectionDirection(1);
44 ObjectCCLSelectionIgnoreSingleCCLFlagOff();
45 VerboseSlicesFlagOff();
46 this->SetK1(vcl_acos(-1.0)/2);
48 //--------------------------------------------------------------------
51 //--------------------------------------------------------------------
52 template <class ImageType>
54 clitk::SliceBySliceRelativePositionFilter<ImageType>::
55 SetInput(const ImageType * image)
57 // Process object is not const-correct so the const casting is required.
58 this->SetNthInput(0, const_cast<ImageType *>(image));
60 //--------------------------------------------------------------------
63 //--------------------------------------------------------------------
64 template <class ImageType>
66 clitk::SliceBySliceRelativePositionFilter<ImageType>::
67 SetInputObject(const ImageType * image)
69 // Process object is not const-correct so the const casting is required.
70 this->SetNthInput(1, const_cast<ImageType *>(image));
72 //--------------------------------------------------------------------
75 //--------------------------------------------------------------------
76 template <class ImageType>
78 clitk::SliceBySliceRelativePositionFilter<ImageType>::
79 PrintOptions(std::ostream & os)
81 os << "Slice direction = " << this->GetDirection() << std::endl
82 << "BG value = " << (int)this->GetBackgroundValue() << std::endl;
83 for(int i=0; i<this->GetNumberOfAngles(); i++) {
84 os << "Orientation = " << this->GetOrientationTypeString()[i] << std::endl;
85 os << "Angles = " << clitk::rad2deg(this->GetAngle1InRad(i))
86 << " " << clitk::rad2deg(this->GetAngle2InRad(i)) << std::endl;
88 os << "InverseOrientationFlag = " << this->GetInverseOrientationFlag() << std::endl
89 << "SpacingFlag = " << this->GetIntermediateSpacingFlag() << std::endl
90 << "Spacing = " << this->GetIntermediateSpacing() << std::endl
91 << "FuzzyThreshold = " << this->GetFuzzyThreshold() << std::endl
92 << "UniqueConnectedComponentBySliceFlag = " << this->GetUniqueConnectedComponentBySliceFlag() << std::endl
93 << "AutoCropFlag = " << this->GetAutoCropFlag() << std::endl
94 << "RemoveObjectFlag= " << this->GetRemoveObjectFlag() << std::endl
95 << "CombineWithOrFlag = " << this->GetCombineWithOrFlag() << std::endl
96 << "UseTheLargestObjectCCLFlag = " << this->GetUseTheLargestObjectCCLFlag() << std::endl
97 << "ObjectCCLSelectionFlag = " << this->GetObjectCCLSelectionFlag() << std::endl
98 << "ObjectCCLSelectionDimension = " << this->GetObjectCCLSelectionDimension() << std::endl
99 << "ObjectCCLSelectionIgnoreSingleCCLFlag = " << this->GetObjectCCLSelectionIgnoreSingleCCLFlag() << std::endl
100 << "IgnoreEmptySliceObjectFlag = " << this->GetIgnoreEmptySliceObjectFlag() << std::endl
101 << "(RP) FastFlag = " << this->GetFastFlag() << std::endl
102 << "(RP) Radius = " << this->GetRadius() << std::endl
103 << "(RP) K1 = " << this->GetK1() << std::endl;
105 //--------------------------------------------------------------------
108 //--------------------------------------------------------------------
109 template <class ImageType>
111 clitk::SliceBySliceRelativePositionFilter<ImageType>::
112 GenerateInputRequestedRegion()
115 itk::ImageToImageFilter<ImageType, ImageType>::GenerateInputRequestedRegion();
116 // Get input pointers and set requested region to common region
117 ImagePointer input1 = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(0));
118 ImagePointer input2 = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(1));
119 input1->SetRequestedRegion(input1->GetLargestPossibleRegion());
120 input2->SetRequestedRegion(input2->GetLargestPossibleRegion());
122 //--------------------------------------------------------------------
125 //--------------------------------------------------------------------
126 template <class ImageType>
128 clitk::SliceBySliceRelativePositionFilter<ImageType>::
129 GenerateOutputInformation()
131 if (this->GetVerboseOptionFlag()) {
135 // if (this->GetFuzzyMapOnlyFlag()) this->ComputeFuzzyMapFlagOn();
138 input = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(0));
139 object = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(1));
140 m_working_object = object;
141 m_working_input = input;
143 //--------------------------------------------------------------------
144 // Resample object to the same spacing than input
145 if (!clitk::HaveSameSpacing<ImageType, ImageType>(object, input)) {
146 this->StartNewStep("Resample object to the same spacing than input");
147 m_working_object = clitk::ResampleImageSpacing<ImageType>(object, input->GetSpacing());
148 this->template StopCurrentStep<ImageType>(m_working_object);
151 //--------------------------------------------------------------------
152 // Resize image according to common area (except in Z)
153 if (!clitk::HaveSameSizeAndSpacing<ImageType, ImageType>(m_working_object, input)) {
154 this->StartNewStep("Resize images (union in XY and like input in Z)");
157 this->StartNewStep("Pad object to the same size than input");
158 m_working_object = clitk::ResizeImageLike<ImageType>(m_working_object,
160 this->GetObjectBackgroundValue());
161 this->template StopCurrentStep<ImageType>(m_working_object);
164 // Compute union of bounding boxes in X and Y
165 static const unsigned int dim = ImageType::ImageDimension;
166 typedef itk::BoundingBox<unsigned long, dim> BBType;
167 typename BBType::Pointer bb1 = BBType::New();
168 ComputeBBFromImageRegion<ImageType>(m_working_object, m_working_object->GetLargestPossibleRegion(), bb1);
169 typename BBType::Pointer bb2 = BBType::New();
170 ComputeBBFromImageRegion<ImageType>(input, input->GetLargestPossibleRegion(), bb2);
171 typename BBType::Pointer bbo = BBType::New();
172 ComputeBBUnion<dim>(bbo, bb1, bb2);
174 //We set Z BB like input
175 typename ImageType::PointType maxs = bbo->GetMaximum();
176 typename ImageType::PointType mins = bbo->GetMinimum();
177 maxs[2] = bb2->GetMaximum()[2];
178 mins[2] = bb2->GetMinimum()[2];
179 bbo->SetMaximum(maxs);
180 bbo->SetMinimum(mins);
183 m_working_input = clitk::ResizeImageLike<ImageType>(input, bbo, this->GetBackgroundValue());
184 m_working_object = clitk::ResizeImageLike<ImageType>(m_working_object,
186 this->GetObjectBackgroundValue());
188 // Index can be negative in some cases, and lead to problem with
189 // some filter. So we correct it.
190 m_working_input = clitk::RemoveNegativeIndexFromRegion<ImageType>(m_working_input);
191 m_working_object = clitk::RemoveNegativeIndexFromRegion<ImageType>(m_working_object);
194 this->template StopCurrentStep<ImageType>(m_working_input);
197 //--------------------------------------------------------------------
199 - extract vector of slices in input, in object
200 - slice by slice rel position
205 //--------------------------------------------------------------------
206 // Extract input slices
207 this->StartNewStep("Extract input slices");
208 typedef clitk::ExtractSliceFilter<ImageType> ExtractSliceFilterType;
209 typename ExtractSliceFilterType::Pointer extractSliceFilter = ExtractSliceFilterType::New();
210 extractSliceFilter->SetInput(m_working_input);
211 extractSliceFilter->SetDirection(GetDirection());
212 extractSliceFilter->Update();
213 typedef typename ExtractSliceFilterType::SliceType SliceType;
214 std::vector<typename SliceType::Pointer> mInputSlices;
215 extractSliceFilter->GetOutputSlices(mInputSlices);
216 this->template StopCurrentStep<SliceType>(mInputSlices[0]);
218 //--------------------------------------------------------------------
219 // Extract object slices
220 this->StartNewStep("Extract object slices");
221 extractSliceFilter = ExtractSliceFilterType::New();
222 extractSliceFilter->SetInput(m_working_object);
223 extractSliceFilter->SetDirection(GetDirection());
224 extractSliceFilter->Update();
225 std::vector<typename SliceType::Pointer> mObjectSlices;
226 extractSliceFilter->GetOutputSlices(mObjectSlices);
227 this->template StopCurrentStep<SliceType>(mObjectSlices[0]);
229 //--------------------------------------------------------------------
230 // Prepare fuzzy slices (if needed)
231 std::vector<typename FloatSliceType::Pointer> mFuzzyMapSlices;
232 mFuzzyMapSlices.resize(mInputSlices.size());
234 //--------------------------------------------------------------------
235 // Perform slice by slice relative position
236 this->StartNewStep("Perform slice by slice relative position ("+toString(mInputSlices.size())+")");
237 for(unsigned int i=0; i<mInputSlices.size(); i++) {
239 // Count the number of CCL (allow to ignore empty slice)
241 mObjectSlices[i] = LabelizeAndCountNumberOfObjects<SliceType>(mObjectSlices[i], 0, true, 1, nb);
243 if (GetVerboseSlicesFlag()) {
244 std::cout << "slice " << i << " nb = " << nb << std::endl;
247 // If no object and empty slices and if we need the full fuzzy map, create a dummy one.
248 if ((nb==0) && (this->GetFuzzyMapOnlyFlag())) {
249 typename FloatSliceType::Pointer one = FloatSliceType::New();
250 one->CopyInformation(mObjectSlices[0]);
251 one->SetRegions(mObjectSlices[0]->GetLargestPossibleRegion());
253 one->FillBuffer(2.0);
254 mFuzzyMapSlices[i] = one;
255 } // End nb==0 && GetComputeFuzzyMapFlag
257 if ((!GetIgnoreEmptySliceObjectFlag()) || (nb!=0)) {
259 // Select or not a single CCL ?
260 if (GetUseTheLargestObjectCCLFlag()) {
261 mObjectSlices[i] = KeepLabels<SliceType>(mObjectSlices[i], 0, 1, 1, 1, true);
264 // Select a single according to a position if more than one CCL
265 if (GetObjectCCLSelectionFlag()) {
266 // if several CCL, choose the most extrema according a direction,
267 // if not -> should we consider this slice ?
269 if (GetObjectCCLSelectionIgnoreSingleCCLFlag()) {
270 mObjectSlices[i] = SetBackground<SliceType, SliceType>(mObjectSlices[i], mObjectSlices[i],
271 1, this->GetBackgroundValue(),
275 int dim = GetObjectCCLSelectionDimension();
276 int direction = GetObjectCCLSelectionDirection();
277 std::vector<typename SliceType::PointType> centroids;
278 ComputeCentroids<SliceType>(mObjectSlices[i], this->GetBackgroundValue(), centroids);
280 for(uint j=1; j<centroids.size(); j++) {
281 if (direction == 1) {
282 if (centroids[j][dim] > centroids[index][dim]) index = j;
285 if (centroids[j][dim] < centroids[index][dim]) index = j;
288 for(uint v=1; v<centroids.size(); v++) {
290 mObjectSlices[i] = SetBackground<SliceType, SliceType>(mObjectSlices[i], mObjectSlices[i],
291 (char)v, this->GetBackgroundValue(),
295 } // end GetbjectCCLSelectionFlag = true
298 typedef clitk::AddRelativePositionConstraintToLabelImageFilter<SliceType> RelPosFilterType;
299 typename RelPosFilterType::Pointer relPosFilter = RelPosFilterType::New();
300 relPosFilter->VerboseStepFlagOff();
301 if (GetVerboseSlicesFlag()) {
302 std::cout << "Slice " << i << std::endl;
303 relPosFilter->VerboseStepFlagOn();
304 //relPosFilter->WriteStepFlagOn();
306 relPosFilter->WriteStepFlagOff();
307 // relPosFilter->VerboseMemoryFlagOn();
308 relPosFilter->SetCurrentStepBaseId(this->GetCurrentStepId()+"-"+toString(i));
309 relPosFilter->SetBackgroundValue(this->GetBackgroundValue());
310 relPosFilter->SetInput(mInputSlices[i]);
311 relPosFilter->SetInputObject(mObjectSlices[i]);
312 relPosFilter->SetRemoveObjectFlag(this->GetRemoveObjectFlag());
313 // This flag (InverseOrientation) *must* be set before
314 // AddOrientation because AddOrientation can change it.
315 relPosFilter->SetInverseOrientationFlag(this->GetInverseOrientationFlag());
316 for(int j=0; j<this->GetNumberOfAngles(); j++) {
317 relPosFilter->AddAnglesInRad(this->GetAngle1InRad(j), this->GetAngle2InRad(j));
319 relPosFilter->SetIntermediateSpacing(this->GetIntermediateSpacing());
320 relPosFilter->SetIntermediateSpacingFlag(this->GetIntermediateSpacingFlag());
321 relPosFilter->SetFuzzyThreshold(this->GetFuzzyThreshold());
322 relPosFilter->AutoCropFlagOff(); // important ! because we join the slices after this loop
323 relPosFilter->SetCombineWithOrFlag(this->GetCombineWithOrFlag());
324 // should we stop after fuzzy map ?
325 relPosFilter->SetFuzzyMapOnlyFlag(this->GetFuzzyMapOnlyFlag());
326 // relPosFilter->SetComputeFuzzyMapFlag(this->GetComputeFuzzyMapFlag());
327 relPosFilter->SetFastFlag(this->GetFastFlag());
328 relPosFilter->SetRadius(this->GetRadius());
329 relPosFilter->SetK1(this->GetK1());
332 relPosFilter->Update();
334 // If we stop after the fuzzy map, store the fuzzy slices
335 if (this->GetFuzzyMapOnlyFlag()) {
336 mFuzzyMapSlices[i] = relPosFilter->GetFuzzyMap();
337 // writeImage<FloatSliceType>(mFuzzyMapSlices[i], "slice_"+toString(i)+".mha");
341 if (!this->GetFuzzyMapOnlyFlag()) {
342 mInputSlices[i] = relPosFilter->GetOutput();
343 // Select main CC if needed
344 if (GetUniqueConnectedComponentBySliceFlag()) {
345 mInputSlices[i] = Labelize<SliceType>(mInputSlices[i], 0, true, 1);
346 mInputSlices[i] = KeepLabels<SliceType>(mInputSlices[i], 0, 1, 1, 1, true);
353 // Select unique CC according to the most in a given direction
354 if (GetUniqueConnectedComponentBySliceAccordingToADirection()) {
356 mInputSlices[i] = LabelizeAndCountNumberOfObjects<SliceType>(mInputSlices[i], 0, true, 1, nb);
357 std::vector<typename ImageType::PointType> & centroids;
362 } // End nb!=0 || GetComputeFuzzyMapFlagOFF
364 } // end for i mInputSlices
367 m_working_input = clitk::JoinSlices<ImageType>(mInputSlices, m_working_input, GetDirection());
368 this->template StopCurrentStep<ImageType>(m_working_input);
370 // Join the fuzzy map if needed
371 if (this->GetFuzzyMapOnlyFlag()) {
372 this->m_FuzzyMap = clitk::JoinSlices<FloatImageType>(mFuzzyMapSlices, m_working_input, GetDirection());
373 this->template StopCurrentStep<FloatImageType>(this->m_FuzzyMap);
374 if (this->GetFuzzyMapOnlyFlag()) return;
377 //--------------------------------------------------------------------
379 if (this->GetAutoCropFlag()) {
380 this->StartNewStep("Final AutoCrop");
381 typedef clitk::AutoCropFilter<ImageType> CropFilterType;
382 typename CropFilterType::Pointer cropFilter = CropFilterType::New();
383 cropFilter->SetInput(m_working_input);
384 cropFilter->ReleaseDataFlagOff();
385 cropFilter->Update();
386 m_working_input = cropFilter->GetOutput();
387 this->template StopCurrentStep<ImageType>(m_working_input);
390 // Update output info
391 this->GetOutput(0)->SetRegions(m_working_input->GetLargestPossibleRegion());
393 //--------------------------------------------------------------------
396 //--------------------------------------------------------------------
397 template <class ImageType>
399 clitk::SliceBySliceRelativePositionFilter<ImageType>::
403 //--------------------------------------------------------------------
404 //--------------------------------------------------------------------
405 // Final Step -> set output
406 //this->SetNthOutput(0, m_working_input);
407 if (this->GetFuzzyMapOnlyFlag()) return; // no output in this case
408 this->GraftOutput(m_working_input);
411 //--------------------------------------------------------------------