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 "clitkCommon.h"
21 #include "clitkBooleanOperatorLabelImageFilter.h"
22 #include "clitkAutoCropFilter.h"
23 #include "clitkResampleImageWithOptionsFilter.h"
24 #include "clitkBooleanOperatorLabelImageFilter.h"
28 #include "itkStatisticsLabelMapFilter.h"
29 #include "itkLabelImageToStatisticsLabelMapFilter.h"
30 #include "itkRegionOfInterestImageFilter.h"
31 #include "itkBinaryThresholdImageFilter.h"
32 #include "itkBinaryErodeImageFilter.h"
33 #include "itkBinaryBallStructuringElement.h"
36 #include "RelativePositionPropImageFilter.h"
38 //--------------------------------------------------------------------
39 template <class ImageType>
40 clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
41 AddRelativePositionConstraintToLabelImageFilter():
43 itk::ImageToImageFilter<ImageType, ImageType>()
45 this->SetNumberOfRequiredInputs(2);
46 SetFuzzyThreshold(0.6);
47 SetBackgroundValue(0);
48 SetObjectBackgroundValue(0);
49 SetOrientationType(LeftTo);
50 ResampleBeforeRelativePositionFilterOn();
51 SetIntermediateSpacing(10);
54 //--------------------------------------------------------------------
57 //--------------------------------------------------------------------
58 template <class ImageType>
60 clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
61 SetInput(const ImageType * image)
63 // Process object is not const-correct so the const casting is required.
64 this->SetNthInput(0, const_cast<ImageType *>(image));
66 //--------------------------------------------------------------------
69 //--------------------------------------------------------------------
70 template <class ImageType>
72 clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
73 SetInputObject(const ImageType * image)
75 // Process object is not const-correct so the const casting is required.
76 this->SetNthInput(1, const_cast<ImageType *>(image));
78 //--------------------------------------------------------------------
81 //--------------------------------------------------------------------
82 template <class ImageType>
84 clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
85 GenerateOutputInformation()
87 ImagePointer input = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(0));
88 ImagePointer outputImage = this->GetOutput(0);
89 outputImage->SetRegions(outputImage->GetLargestPossibleRegion());
91 //--------------------------------------------------------------------
94 //--------------------------------------------------------------------
95 template <class ImageType>
97 clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
98 GenerateInputRequestedRegion()
101 itk::ImageToImageFilter<ImageType, ImageType>::GenerateInputRequestedRegion();
102 // Get input pointers and set requested region to common region
103 ImagePointer input1 = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(0));
104 ImagePointer input2 = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(1));
105 input1->SetRequestedRegion(input1->GetLargestPossibleRegion());
106 input2->SetRequestedRegion(input2->GetLargestPossibleRegion());
108 //--------------------------------------------------------------------
111 //--------------------------------------------------------------------
112 template <class ImageType>
114 clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
117 SetOrientationType(Angle);
120 //--------------------------------------------------------------------
123 //--------------------------------------------------------------------
124 template <class ImageType>
126 clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
129 SetOrientationType(Angle);
132 //--------------------------------------------------------------------
135 //--------------------------------------------------------------------
136 template <class ImageType>
138 clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
139 SetOrientationType(OrientationTypeEnumeration orientation)
141 m_OrientationType = orientation;
142 switch (m_OrientationType) {
143 case LeftTo: m_Angle1 = clitk::deg2rad(0); m_Angle2 = clitk::deg2rad(0); break;
144 case RightTo: m_Angle1 = clitk::deg2rad(180); m_Angle2 = clitk::deg2rad(0); break;
145 case AntTo: m_Angle1 = clitk::deg2rad(90); m_Angle2 = clitk::deg2rad(0); break;
146 case PostTo: m_Angle1 = clitk::deg2rad(-90); m_Angle2 = clitk::deg2rad(0); break;
147 case InfTo: m_Angle1 = clitk::deg2rad(0); m_Angle2 = clitk::deg2rad(90); break;
148 case SupTo: m_Angle1 = clitk::deg2rad(0); m_Angle2 = clitk::deg2rad(-90); break;
160 //--------------------------------------------------------------------
163 //--------------------------------------------------------------------
164 template <class ImageType>
166 clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
170 input = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(0));
171 object = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(1));
173 //--------------------------------------------------------------------
174 //--------------------------------------------------------------------
175 static const unsigned int dim = ImageType::ImageDimension;
176 StartNewStep("Initial resample and pad");
178 if (m_ResampleBeforeRelativePositionFilter) {
179 typedef clitk::ResampleImageWithOptionsFilter<ImageType> ResampleFilterType;
180 typename ResampleFilterType::Pointer resampleFilter = ResampleFilterType::New();
181 resampleFilter->SetInput(object);
182 resampleFilter->SetOutputIsoSpacing(m_IntermediateSpacing);
183 resampleFilter->SetGaussianFilteringEnabled(false);
184 // resampleFilter->SetVerboseOptions(true);
185 resampleFilter->Update();
186 working_image = resampleFilter->GetOutput();
189 working_image = object;
192 // Step 2: object pad to input image -> we want to compute the
193 // relative position for each point belonging to the input image
194 // domain, so we have to extend (pad) the object image to fit the
196 if (!clitk::HaveSameSizeAndSpacing<ImageType, ImageType>(input, working_image)) {
197 typename ImageType::Pointer output = ImageType::New();
199 for(unsigned int i=0; i<dim; i++) {
200 size[i] = lrint((input->GetLargestPossibleRegion().GetSize()[i]*input->GetSpacing()[i])/(double)working_image->GetSpacing()[i]);
203 region.SetSize(size);
204 // output->SetLargestPossibleRegion(region);
205 output->SetRegions(region);
206 output->SetSpacing(working_image->GetSpacing());
207 output->SetOrigin(input->GetOrigin());
209 output->FillBuffer(m_BackgroundValue);
210 typename PadFilterType::Pointer padFilter = PadFilterType::New();
211 typename PadFilterType::InputImageIndexType index;
212 for(unsigned int i=0; i<dim; i++) {
213 index[i] = lrint((working_image->GetOrigin()[i] - input->GetOrigin()[i])/working_image->GetSpacing()[i]);
215 padFilter->SetSourceImage(working_image);
216 padFilter->SetDestinationImage(output);
217 padFilter->SetDestinationIndex(index);
218 padFilter->SetSourceRegion(working_image->GetLargestPossibleRegion());
220 working_image = padFilter->GetOutput();
223 // DD("[debug] RelPos : same size and spacing : no padding");
225 // Keep object image (with resampline and pad)
226 object_resampled = working_image;
227 StopCurrentStep<ImageType>(working_image);
229 // Step 3: compute rel pos in object
230 StartNewStep("Relative Position Map");
231 typedef itk::RelativePositionPropImageFilter<ImageType, FloatImageType> RelPosFilterType;
232 typename RelPosFilterType::Pointer relPosFilter = RelPosFilterType::New();
233 relPosFilter->SetInput(working_image);
234 relPosFilter->SetAlpha1(m_Angle1); // xy plane
235 relPosFilter->SetAlpha2(m_Angle2);
236 relPosFilter->SetK1(M_PI/2.0); // Opening parameter, default = pi/2
237 relPosFilter->SetFast(true);
238 relPosFilter->SetRadius(1); // seems sufficient in this cas
239 // relPosFilter->SetVerboseProgress(true);
240 relPosFilter->Update();
241 relPos = relPosFilter->GetOutput();
242 StopCurrentStep<FloatImageType>(relPos);
244 //--------------------------------------------------------------------
245 //--------------------------------------------------------------------
246 StartNewStep("Map Threshold");
248 typedef itk::BinaryThresholdImageFilter<FloatImageType, ImageType> BinaryThresholdImageFilterType;
249 typename BinaryThresholdImageFilterType::Pointer thresholdFilter = BinaryThresholdImageFilterType::New();
250 thresholdFilter->SetInput(relPos);
251 thresholdFilter->SetOutsideValue(m_BackgroundValue);
252 thresholdFilter->SetInsideValue(m_BackgroundValue+1);
253 thresholdFilter->SetLowerThreshold(m_FuzzyThreshold);
254 thresholdFilter->Update();
255 working_image = thresholdFilter->GetOutput();
256 StopCurrentStep<ImageType>(working_image);
258 //--------------------------------------------------------------------
259 //--------------------------------------------------------------------
260 StartNewStep("Post Processing: erosion with initial mask");
261 // Step 2 : erosion with initial mask to exclude pixels that were
262 // inside the resampled version and outside the original mask
263 typedef itk::BinaryBallStructuringElement<unsigned int, ImageDimension> StructuringElementType;
264 StructuringElementType kernel;
266 kernel.CreateStructuringElement();
267 typedef itk::BinaryErodeImageFilter<ImageType, ImageType, StructuringElementType> ErodeFilterType;
268 typename ErodeFilterType::Pointer erodeFilter = ErodeFilterType::New();
269 erodeFilter->SetInput(working_image);
270 erodeFilter->SetKernel(kernel);
271 erodeFilter->SetBackgroundValue(m_BackgroundValue);
272 erodeFilter->SetErodeValue(m_BackgroundValue+1);
273 erodeFilter->Update();
274 working_image = erodeFilter->GetOutput();
275 StopCurrentStep<ImageType>(working_image);
277 //--------------------------------------------------------------------
278 //--------------------------------------------------------------------
279 // Step 5: resample to initial spacing
280 if (m_ResampleBeforeRelativePositionFilter) {
281 StartNewStep("Resample to get the same sampling than input");
282 typedef clitk::ResampleImageWithOptionsFilter<ImageType> ResampleFilterType;
283 typename ResampleFilterType::Pointer resampleFilter = ResampleFilterType::New();
284 resampleFilter->SetDefaultPixelValue(m_BackgroundValue);
285 resampleFilter->SetInput(working_image);
286 resampleFilter->SetOutputSpacing(input->GetSpacing());
287 resampleFilter->SetGaussianFilteringEnabled(false);
288 // resampleFilter->SetVerboseOptions(true);
289 resampleFilter->Update();
290 working_image = resampleFilter->GetOutput();
291 StopCurrentStep<ImageType>(working_image);
294 //--------------------------------------------------------------------
295 //--------------------------------------------------------------------
296 // Pre Step 6: pad if not the same size : it can occur when downsample and upsample
297 if (working_image->GetLargestPossibleRegion() != input->GetLargestPossibleRegion()) {
298 StartNewStep("Pad to get the same size than input");
299 typename ImageType::Pointer temp = ImageType::New();
300 temp->CopyInformation(input);
301 temp->SetRegions(input->GetLargestPossibleRegion()); // Do not forget !!
303 temp->FillBuffer(m_BackgroundValue);
304 typename PadFilterType::Pointer padFilter2 = PadFilterType::New(); // if yes : redo relpos
305 padFilter2->SetSourceImage(working_image);
306 padFilter2->SetDestinationImage(temp);
307 padFilter2->SetDestinationIndex(input->GetLargestPossibleRegion().GetIndex());
308 padFilter2->SetSourceRegion(working_image->GetLargestPossibleRegion());
309 padFilter2->Update();
310 working_image = padFilter2->GetOutput();
311 StopCurrentStep<ImageType>(working_image);
314 //DD("[debug] Rel Pos : no padding after");
317 //--------------------------------------------------------------------
318 //--------------------------------------------------------------------
319 // Step 6: combine input+thresholded relpos
320 StartNewStep("Combine with initial input (boolean And)");
321 typedef clitk::BooleanOperatorLabelImageFilter<ImageType> BoolFilterType;
322 typename BoolFilterType::Pointer combineFilter = BoolFilterType::New();
323 writeImage<ImageType>(input, "i.mhd");
324 writeImage<ImageType>(working_image, "w.mhd");
325 combineFilter->SetBackgroundValue(m_BackgroundValue);
326 combineFilter->SetBackgroundValue1(m_BackgroundValue);
327 combineFilter->SetBackgroundValue2(m_BackgroundValue);
328 combineFilter->SetForegroundValue(m_BackgroundValue+1);
329 combineFilter->SetInput1(input);
330 combineFilter->SetInput2(working_image);
331 combineFilter->SetOperationType(BoolFilterType::And);
332 combineFilter->InPlaceOn();
333 combineFilter->Update();
334 working_image = combineFilter->GetOutput();
335 // writeImage<ImageType>(working_image, "res.mhd");
337 combineFilter = BoolFilterType::New();
338 combineFilter->SetInput1(working_image);
339 combineFilter->SetInput2(object);
340 combineFilter->SetOperationType(BoolFilterType::AndNot);
341 combineFilter->InPlaceOn();
342 combineFilter->Update();
344 working_image = combineFilter->GetOutput();
345 StopCurrentStep<ImageType>(working_image);
347 //--------------------------------------------------------------------
348 //--------------------------------------------------------------------
351 StartNewStep("Final AutoCrop");
352 typedef clitk::AutoCropFilter<ImageType> CropFilterType;
353 typename CropFilterType::Pointer cropFilter = CropFilterType::New();
354 cropFilter->SetInput(working_image);
355 cropFilter->ReleaseDataFlagOff();
356 cropFilter->Update();
357 working_image = cropFilter->GetOutput();
358 StopCurrentStep<ImageType>(working_image);
361 //--------------------------------------------------------------------
362 //--------------------------------------------------------------------
364 // Final Step -> set output
365 this->SetNthOutput(0, working_image);
368 //--------------------------------------------------------------------