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();
46 //--------------------------------------------------------------------
49 //--------------------------------------------------------------------
50 template <class ImageType>
52 clitk::SliceBySliceRelativePositionFilter<ImageType>::
53 SetInput(const ImageType * image)
55 // Process object is not const-correct so the const casting is required.
56 this->SetNthInput(0, const_cast<ImageType *>(image));
58 //--------------------------------------------------------------------
61 //--------------------------------------------------------------------
62 template <class ImageType>
64 clitk::SliceBySliceRelativePositionFilter<ImageType>::
65 SetInputObject(const ImageType * image)
67 // Process object is not const-correct so the const casting is required.
68 this->SetNthInput(1, const_cast<ImageType *>(image));
70 //--------------------------------------------------------------------
73 //--------------------------------------------------------------------
74 template <class ImageType>
76 clitk::SliceBySliceRelativePositionFilter<ImageType>::
77 PrintOptions(std::ostream & os)
79 os << "Slice direction = " << this->GetDirection() << std::endl
80 << "BG value = " << this->GetBackgroundValue() << std::endl;
81 for(int i=0; i<this->GetNumberOfAngles(); i++) {
82 os << "Orientation = " << this->GetOrientationTypeString()[i] << std::endl;
83 os << "Angles = " << clitk::rad2deg(this->GetAngle1InRad(i))
84 << " " << clitk::rad2deg(this->GetAngle2InRad(i)) << std::endl;
86 os << "InverseOrientationFlag = " << this->GetInverseOrientationFlag() << std::endl
87 << "SpacingFlag = " << this->GetIntermediateSpacingFlag() << std::endl
88 << "Spacing = " << this->GetIntermediateSpacing() << std::endl
89 << "FuzzyThreshold = " << this->GetFuzzyThreshold() << std::endl
90 << "UniqueConnectedComponentBySliceFlag = " << this->GetUniqueConnectedComponentBySliceFlag() << std::endl
91 << "AutoCropFlag = " << this->GetAutoCropFlag() << std::endl
92 << "RemoveObjectFlag= " << this->GetRemoveObjectFlag() << std::endl
93 << "CombineWithOrFlag = " << this->GetCombineWithOrFlag() << std::endl
94 << "UseTheLargestObjectCCLFlag = " << this->GetUseTheLargestObjectCCLFlag() << std::endl
95 << "ObjectCCLSelectionFlag = " << this->GetObjectCCLSelectionFlag() << std::endl
96 << "ObjectCCLSelectionDimension = " << this->GetObjectCCLSelectionDimension() << std::endl
97 << "ObjectCCLSelectionIgnoreSingleCCLFlag = " << this->GetObjectCCLSelectionIgnoreSingleCCLFlag() << std::endl
98 << "IgnoreEmptySliceObjectFlag = " << this->GetIgnoreEmptySliceObjectFlag() << std::endl;
100 //--------------------------------------------------------------------
103 //--------------------------------------------------------------------
104 template <class ImageType>
106 clitk::SliceBySliceRelativePositionFilter<ImageType>::
107 GenerateInputRequestedRegion()
110 itk::ImageToImageFilter<ImageType, ImageType>::GenerateInputRequestedRegion();
111 // Get input pointers and set requested region to common region
112 ImagePointer input1 = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(0));
113 ImagePointer input2 = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(1));
114 input1->SetRequestedRegion(input1->GetLargestPossibleRegion());
115 input2->SetRequestedRegion(input2->GetLargestPossibleRegion());
117 //--------------------------------------------------------------------
120 //--------------------------------------------------------------------
121 template <class ImageType>
123 clitk::SliceBySliceRelativePositionFilter<ImageType>::
124 GenerateOutputInformation()
126 if (this->GetVerboseOptionFlag()) {
131 input = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(0));
132 object = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(1));
134 //--------------------------------------------------------------------
135 // Resample object to the same spacing than input
136 if (!clitk::HaveSameSpacing<ImageType, ImageType>(object, input)) {
137 this->StartNewStep("Resample object to the same spacing than input");
138 m_working_object = clitk::ResampleImageSpacing<ImageType>(object, input->GetSpacing());
139 this->template StopCurrentStep<ImageType>(m_working_object);
142 m_working_object = object;
145 //--------------------------------------------------------------------
146 // Pad object to the same size than input
147 if (!clitk::HaveSameSizeAndSpacing<ImageType, ImageType>(m_working_object, input)) {
148 this->StartNewStep("Pad object to the same size than input");
149 m_working_object = clitk::ResizeImageLike<ImageType>(m_working_object,
151 this->GetObjectBackgroundValue());
152 this->template StopCurrentStep<ImageType>(m_working_object);
158 - extract vector of slices in input, in object
159 - slice by slice rel position
165 //--------------------------------------------------------------------
166 // Extract input slices
167 this->StartNewStep("Extract input slices");
168 typedef clitk::ExtractSliceFilter<ImageType> ExtractSliceFilterType;
169 typename ExtractSliceFilterType::Pointer extractSliceFilter = ExtractSliceFilterType::New();
170 extractSliceFilter->SetInput(input);
171 extractSliceFilter->SetDirection(GetDirection());
172 extractSliceFilter->Update();
173 typedef typename ExtractSliceFilterType::SliceType SliceType;
174 std::vector<typename SliceType::Pointer> mInputSlices;
175 extractSliceFilter->GetOutputSlices(mInputSlices);
176 this->template StopCurrentStep<SliceType>(mInputSlices[0]);
178 //--------------------------------------------------------------------
179 // Extract object slices
180 this->StartNewStep("Extract object slices");
181 extractSliceFilter = ExtractSliceFilterType::New();
182 extractSliceFilter->SetInput(m_working_object);//object);
183 extractSliceFilter->SetDirection(GetDirection());
184 extractSliceFilter->Update();
185 std::vector<typename SliceType::Pointer> mObjectSlices;
186 extractSliceFilter->GetOutputSlices(mObjectSlices);
187 this->template StopCurrentStep<SliceType>(mObjectSlices[0]);
189 //--------------------------------------------------------------------
190 // Prepare fuzzy slices (if needed)
191 std::vector<typename FloatSliceType::Pointer> mFuzzyMapSlices;
192 mFuzzyMapSlices.resize(mInputSlices.size());
194 //--------------------------------------------------------------------
195 // Perform slice by slice relative position
196 this->StartNewStep("Perform slice by slice relative position");
197 for(unsigned int i=0; i<mInputSlices.size(); i++) {
199 // Count the number of CCL (allow to ignore empty slice)
201 mObjectSlices[i] = LabelizeAndCountNumberOfObjects<SliceType>(mObjectSlices[i], 0, true, 1, nb);
203 // If no object and empty slices :
204 if ((nb==0) && (this->GetFuzzyMapOnlyFlag())) {
205 typename FloatSliceType::Pointer one = FloatSliceType::New();
206 one->CopyInformation(mObjectSlices[0]);
207 one->SetRegions(mObjectSlices[0]->GetLargestPossibleRegion());
209 one->FillBuffer(2.0);
210 mFuzzyMapSlices[i] = one;
213 if ((!GetIgnoreEmptySliceObjectFlag()) || (nb!=0)) {
215 // Select or not a single CCL ?
216 if (GetUseTheLargestObjectCCLFlag()) {
217 mObjectSlices[i] = KeepLabels<SliceType>(mObjectSlices[i], 0, 1, 1, 1, true);
220 // Select a single according to a position if more than one CCL
221 if (GetObjectCCLSelectionFlag()) {
222 // if several CCL, choose the most extrema according a direction,
223 // if not -> should we consider this slice ?
225 if (GetObjectCCLSelectionIgnoreSingleCCLFlag()) {
226 mObjectSlices[i] = SetBackground<SliceType, SliceType>(mObjectSlices[i], mObjectSlices[i],
227 1, this->GetBackgroundValue(),
231 int dim = GetObjectCCLSelectionDimension();
232 int direction = GetObjectCCLSelectionDirection();
233 std::vector<typename SliceType::PointType> centroids;
234 ComputeCentroids<SliceType>(mObjectSlices[i], this->GetBackgroundValue(), centroids);
236 for(uint j=1; j<centroids.size(); j++) {
237 if (direction == 1) {
238 if (centroids[j][dim] > centroids[index][dim]) index = j;
241 if (centroids[j][dim] < centroids[index][dim]) index = j;
244 for(uint v=1; v<centroids.size(); v++) {
246 mObjectSlices[i] = SetBackground<SliceType, SliceType>(mObjectSlices[i], mObjectSlices[i],
247 (char)v, this->GetBackgroundValue(),
251 } // end GetbjectCCLSelectionFlag = true
254 typedef clitk::AddRelativePositionConstraintToLabelImageFilter<SliceType> RelPosFilterType;
255 typename RelPosFilterType::Pointer relPosFilter = RelPosFilterType::New();
257 relPosFilter->VerboseStepFlagOff();
258 relPosFilter->WriteStepFlagOff();
259 relPosFilter->SetBackgroundValue(this->GetBackgroundValue());
260 relPosFilter->SetInput(mInputSlices[i]);
261 relPosFilter->SetInputObject(mObjectSlices[i]);
262 relPosFilter->SetRemoveObjectFlag(this->GetRemoveObjectFlag());
263 // This flag (InverseOrientation) *must* be set before
264 // AddOrientation because AddOrientation can change it.
265 relPosFilter->SetInverseOrientationFlag(this->GetInverseOrientationFlag());
266 for(int j=0; j<this->GetNumberOfAngles(); j++) {
267 // relPosFilter->AddOrientationTypeString(this->GetOrientationTypeString(j));
268 relPosFilter->AddAnglesInRad(this->GetAngle1InRad(j), this->GetAngle2InRad(j));
269 // DD(this->GetOrientationTypeString(j));
271 // DD(this->GetInverseOrientationFlag());
272 //relPosFilter->SetOrientationType(this->GetOrientationType());
273 relPosFilter->SetIntermediateSpacing(this->GetIntermediateSpacing());
274 relPosFilter->SetIntermediateSpacingFlag(this->GetIntermediateSpacingFlag());
275 relPosFilter->SetFuzzyThreshold(this->GetFuzzyThreshold());
276 relPosFilter->AutoCropFlagOff(); // important ! because we join the slices after this loop
277 relPosFilter->SetCombineWithOrFlag(this->GetCombineWithOrFlag());
279 // should we stop after fuzzy map ?
280 relPosFilter->SetFuzzyMapOnlyFlag(this->GetFuzzyMapOnlyFlag());
283 relPosFilter->Update();
285 // If we stop after the fuzzy map, store the fuzzy slices
286 if (this->GetFuzzyMapOnlyFlag()) {
287 mFuzzyMapSlices[i] = relPosFilter->GetFuzzyMap();
288 // writeImage<FloatSliceType>(mFuzzyMapSlices[i], "slice_"+toString(i)+".mha");
291 mInputSlices[i] = relPosFilter->GetOutput();
292 // Select main CC if needed
293 if (GetUniqueConnectedComponentBySliceFlag()) {
294 mInputSlices[i] = Labelize<SliceType>(mInputSlices[i], 0, true, 1);
295 mInputSlices[i] = KeepLabels<SliceType>(mInputSlices[i], 0, 1, 1, 1, true);
302 // Select unique CC according to the most in a given direction
303 if (GetUniqueConnectedComponentBySliceAccordingToADirection()) {
305 mInputSlices[i] = LabelizeAndCountNumberOfObjects<SliceType>(mInputSlices[i], 0, true, 1, nb);
306 std::vector<typename ImageType::PointType> & centroids;
313 // Join the fuzzy map if needed
314 if (this->GetFuzzyMapOnlyFlag()) {
315 this->m_FuzzyMap = clitk::JoinSlices<FloatImageType>(mFuzzyMapSlices, input, GetDirection());
316 this->template StopCurrentStep<FloatImageType>(this->m_FuzzyMap);
321 m_working_input = clitk::JoinSlices<ImageType>(mInputSlices, input, GetDirection());
322 this->template StopCurrentStep<ImageType>(m_working_input);
324 //--------------------------------------------------------------------
326 if (this->GetAutoCropFlag()) {
327 this->StartNewStep("Final AutoCrop");
328 typedef clitk::AutoCropFilter<ImageType> CropFilterType;
329 typename CropFilterType::Pointer cropFilter = CropFilterType::New();
330 cropFilter->SetInput(m_working_input);
331 cropFilter->ReleaseDataFlagOff();
332 cropFilter->Update();
333 m_working_input = cropFilter->GetOutput();
334 this->template StopCurrentStep<ImageType>(m_working_input);
337 // Update output info
338 this->GetOutput(0)->SetRegions(m_working_input->GetLargestPossibleRegion());
340 //--------------------------------------------------------------------
343 //--------------------------------------------------------------------
344 template <class ImageType>
346 clitk::SliceBySliceRelativePositionFilter<ImageType>::
350 //--------------------------------------------------------------------
351 //--------------------------------------------------------------------
352 // Final Step -> set output
353 //this->SetNthOutput(0, m_working_input);
354 if (this->GetFuzzyMapOnlyFlag()) return; // no output in this case
355 this->GraftOutput(m_working_input);
358 //--------------------------------------------------------------------