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>
34 #include <itkAddImageFilter.h>
35 #include <itkDivideByConstantImageFilter.h>
38 #include "RelativePositionPropImageFilter.h"
40 //--------------------------------------------------------------------
41 template <class ImageType>
42 clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
43 AddRelativePositionConstraintToLabelImageFilter():
45 itk::ImageToImageFilter<ImageType, ImageType>()
47 this->SetNumberOfRequiredInputs(2);
48 SetFuzzyThreshold(0.6);
49 SetBackgroundValue(0);
50 SetObjectBackgroundValue(0);
51 ClearOrientationType();
52 IntermediateSpacingFlagOn();
53 SetIntermediateSpacing(10);
55 InverseOrientationFlagOff();
57 CombineWithOrFlagOff();
61 //--------------------------------------------------------------------
64 //--------------------------------------------------------------------
65 template <class ImageType>
67 clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
68 SetInput(const ImageType * image)
70 // Process object is not const-correct so the const casting is required.
71 this->SetNthInput(0, const_cast<ImageType *>(image));
73 //--------------------------------------------------------------------
76 //--------------------------------------------------------------------
77 template <class ImageType>
79 clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
80 SetInputObject(const ImageType * image)
82 // Process object is not const-correct so the const casting is required.
83 this->SetNthInput(1, const_cast<ImageType *>(image));
85 //--------------------------------------------------------------------
88 //--------------------------------------------------------------------
89 template <class ImageType>
91 clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
92 ClearOrientationType()
94 m_OrientationTypeString.clear();
95 m_OrientationType.clear();
99 //--------------------------------------------------------------------
102 //--------------------------------------------------------------------
103 template <class ImageType>
105 clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
108 return m_OrientationType.size();
110 //--------------------------------------------------------------------
113 //--------------------------------------------------------------------
114 template <class ImageType>
116 clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
117 AddOrientationTypeString(std::string t)
119 m_OrientationTypeString.push_back(t);
121 case 'L' : AddOrientationType(AtLeftTo); break;
122 case 'R' : AddOrientationType(AtRightTo);break;
123 case 'A' : AddOrientationType(AntTo);break;
124 case 'P' : AddOrientationType(PostTo);break;
125 case 'S' : AddOrientationType(SupTo);break;
126 case 'I' : AddOrientationType(InfTo);break;
128 if (t == "NotLeftTo") { AddOrientationType(AtLeftTo); InverseOrientationFlagOn(); break; }
129 if (t == "NotRightTo") { AddOrientationType(AtRightTo); InverseOrientationFlagOn(); break; }
130 if (t == "NotAntTo") { AddOrientationType(AntTo); InverseOrientationFlagOn(); break; }
131 if (t == "NotPostTo") { AddOrientationType(PostTo); InverseOrientationFlagOn(); break; }
132 if (t == "NotSupTo") { AddOrientationType(SupTo); InverseOrientationFlagOn(); break; }
133 if (t == "NotInfTo") { AddOrientationType(InfTo); InverseOrientationFlagOn(); break; }
134 default: clitkExceptionMacro("Error, you must provide L,R or A,P or S,I (or NotLeftTo, NotRightTo etc)");
137 //--------------------------------------------------------------------
140 //--------------------------------------------------------------------
141 template <class ImageType>
143 clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
144 GenerateOutputInformation()
146 ImagePointer input = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(0));
147 ImagePointer outputImage = this->GetOutput(0);
148 outputImage->SetRegions(outputImage->GetLargestPossibleRegion());
150 //--------------------------------------------------------------------
153 //--------------------------------------------------------------------
154 template <class ImageType>
156 clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
157 GenerateInputRequestedRegion()
160 itk::ImageToImageFilter<ImageType, ImageType>::GenerateInputRequestedRegion();
161 // Get input pointers and set requested region to common region
162 ImagePointer input1 = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(0));
163 ImagePointer input2 = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(1));
164 input1->SetRequestedRegion(input1->GetLargestPossibleRegion());
165 input2->SetRequestedRegion(input2->GetLargestPossibleRegion());
167 //--------------------------------------------------------------------
170 //--------------------------------------------------------------------
171 template <class ImageType>
173 clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
174 AddAngles(double a, double b)
176 AddOrientationTypeString("Angle");
177 m_Angle1.push_back(a);
178 m_Angle2.push_back(b);
180 //--------------------------------------------------------------------
183 //--------------------------------------------------------------------
184 template <class ImageType>
186 clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
187 AddOrientationType(OrientationTypeEnumeration orientation)
189 m_OrientationType.push_back(orientation);
190 switch (orientation) {
192 m_Angle1.push_back(clitk::deg2rad(0));
193 m_Angle2.push_back(clitk::deg2rad(0));
196 m_Angle1.push_back(clitk::deg2rad(180));
197 m_Angle2.push_back(clitk::deg2rad(0));
200 m_Angle1.push_back(clitk::deg2rad(90));
201 m_Angle2.push_back(clitk::deg2rad(0));
204 m_Angle1.push_back(clitk::deg2rad(-90));
205 m_Angle2.push_back(clitk::deg2rad(0));
208 m_Angle1.push_back(clitk::deg2rad(0));
209 m_Angle2.push_back(clitk::deg2rad(90));
212 m_Angle1.push_back(clitk::deg2rad(0));
213 m_Angle2.push_back(clitk::deg2rad(-90));
226 //--------------------------------------------------------------------
229 //--------------------------------------------------------------------
230 template <class ImageType>
232 clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
235 if (GetNumberOfAngles() <1) {
236 clitkExceptionMacro("Add at least one orientation type");
240 input = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(0));
241 object = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(1));
242 static const unsigned int dim = ImageType::ImageDimension;
244 // Step 2: object pad to input image -> we want to compute the
245 // relative position for each point belonging to the input image
246 // domain, so we have to extend (pad) the object image to fit the
248 working_image = object;
249 if (!clitk::HaveSameSizeAndSpacing<ImageType, ImageType>(input, working_image)) {
250 StartNewStep("Pad (resize) object to input size");
252 if (0) { // OLD VERSION (TO REMOVE)
253 StartNewStep("Pad object to image size");
254 typename ImageType::Pointer output = ImageType::New();
256 for(unsigned int i=0; i<dim; i++) {
257 size[i] = lrint((input->GetLargestPossibleRegion().GetSize()[i]*
258 input->GetSpacing()[i])/(double)working_image->GetSpacing()[i]);
261 // The index of the input is not necessarily zero, so we have to
262 // take it into account (not done)
264 IndexType index = input->GetLargestPossibleRegion().GetIndex();
265 region.SetSize(size);
266 for(unsigned int i=0; i<dim; i++) {
268 std::cerr << "Index diff from zero : " << index << ". not done yet !" << std::endl;
272 // output->SetLargestPossibleRegion(region);
273 output->SetRegions(region);
274 output->SetSpacing(working_image->GetSpacing());
275 PointType origin = input->GetOrigin();
276 for(unsigned int i=0; i<dim; i++) {
277 origin[i] = index[i]*input->GetSpacing()[i] + input->GetOrigin()[i];
279 output->SetOrigin(origin);
280 // output->SetOrigin(input->GetOrigin());
283 output->FillBuffer(m_BackgroundValue);
284 typename PasteFilterType::Pointer padFilter = PasteFilterType::New();
285 // typename PasteFilterType::InputImageIndexType index;
286 for(unsigned int i=0; i<dim; i++) {
287 index[i] = -index[i]*input->GetSpacing()[i]/(double)working_image->GetSpacing()[i]
288 + lrint((working_image->GetOrigin()[i] - input->GetOrigin()[i])/working_image->GetSpacing()[i]);
290 padFilter->SetSourceImage(working_image);
291 padFilter->SetDestinationImage(output);
292 padFilter->SetDestinationIndex(index);
293 padFilter->SetSourceRegion(working_image->GetLargestPossibleRegion());
295 working_image = padFilter->GetOutput();
298 // Resize object like input
299 working_image = clitk::ResizeImageLike<ImageType>(working_image, input, GetBackgroundValue());
300 StopCurrentStep<ImageType>(working_image);
303 //--------------------------------------------------------------------
304 //--------------------------------------------------------------------
306 if (m_IntermediateSpacingFlag) {
307 StartNewStep("Resample object to intermediate spacing");
308 typedef clitk::ResampleImageWithOptionsFilter<ImageType> ResampleFilterType;
309 typename ResampleFilterType::Pointer resampleFilter = ResampleFilterType::New();
310 resampleFilter->SetInput(working_image);
311 resampleFilter->SetDefaultPixelValue(0);
312 resampleFilter->SetOutputIsoSpacing(m_IntermediateSpacing);
313 resampleFilter->SetGaussianFilteringEnabled(false);
314 // resampleFilter->SetVerboseOptions(true);
315 resampleFilter->Update();
316 working_image = resampleFilter->GetOutput();
317 StopCurrentStep<ImageType>(working_image);
320 // Keep object image (with resampline and pad)
321 object_resampled = working_image;
323 // Step 3: compute rel pos in object
324 StartNewStep("Relative Position Map");
325 typedef itk::RelativePositionPropImageFilter<ImageType, FloatImageType> RelPosFilterType;
326 typename RelPosFilterType::Pointer relPosFilter;
328 typename FloatImageType::Pointer m_FuzzyMap;
329 for(int i=0; i<GetNumberOfAngles(); i++) {
331 relPosFilter = RelPosFilterType::New();
332 relPosFilter->SetInput(working_image);
333 relPosFilter->SetAlpha1(m_Angle1[i]); // xy plane
334 relPosFilter->SetAlpha2(m_Angle2[i]);
335 relPosFilter->SetK1(M_PI/2.0); // Opening parameter, default = pi/2
336 relPosFilter->SetFast(true);
337 relPosFilter->SetRadius(1); // seems sufficient in this case
338 // relPosFilter->SetVerboseProgress(true);
339 relPosFilter->Update();
340 relPos = relPosFilter->GetOutput();
342 if (GetNumberOfAngles() != 1) {
343 // Creation of the first m_FuzzyMap
345 m_FuzzyMap = clitk::NewImageLike<FloatImageType>(relPos, true);
346 m_FuzzyMap->FillBuffer(0.0);
349 // Add to current fuzzy map
350 typedef itk::AddImageFilter<FloatImageType, FloatImageType, FloatImageType> AddImageFilter;
351 typename AddImageFilter::Pointer addFilter = AddImageFilter::New();
352 addFilter->SetInput1(m_FuzzyMap);
353 addFilter->SetInput2(relPos);
355 m_FuzzyMap = addFilter->GetOutput();
357 else m_FuzzyMap = relPos;
360 // Divide by the number of relpos
361 if (GetNumberOfAngles() != 1) {
362 typedef itk::DivideByConstantImageFilter<FloatImageType, float, FloatImageType> DivideFilter;
363 typename DivideFilter::Pointer divideFilter = DivideFilter::New();
364 divideFilter->SetInput(m_FuzzyMap);
365 divideFilter->SetConstant(GetNumberOfAngles());
366 divideFilter->Update();
367 m_FuzzyMap = divideFilter->GetOutput();
371 StopCurrentStep<FloatImageType>(relPos);
373 //--------------------------------------------------------------------
374 //--------------------------------------------------------------------
375 StartNewStep("Map Threshold");
377 typedef itk::BinaryThresholdImageFilter<FloatImageType, ImageType> BinaryThresholdImageFilterType;
378 typename BinaryThresholdImageFilterType::Pointer thresholdFilter = BinaryThresholdImageFilterType::New();
379 thresholdFilter->SetInput(relPos);
380 thresholdFilter->SetOutsideValue(m_BackgroundValue);
381 thresholdFilter->SetInsideValue(m_BackgroundValue+1);
382 thresholdFilter->SetLowerThreshold(m_FuzzyThreshold);
383 thresholdFilter->Update();
384 working_image = thresholdFilter->GetOutput();
385 StopCurrentStep<ImageType>(working_image);
387 //--------------------------------------------------------------------
388 //--------------------------------------------------------------------
389 StartNewStep("Post Processing: erosion with initial mask");
390 // Step 2 : erosion with initial mask to exclude pixels that were
391 // inside the resampled version and outside the original mask
392 typedef itk::BinaryBallStructuringElement<unsigned int, ImageDimension> StructuringElementType;
393 StructuringElementType kernel;
395 kernel.CreateStructuringElement();
396 typedef itk::BinaryErodeImageFilter<ImageType, ImageType, StructuringElementType> ErodeFilterType;
397 typename ErodeFilterType::Pointer erodeFilter = ErodeFilterType::New();
398 erodeFilter->SetInput(working_image);
399 erodeFilter->SetKernel(kernel);
400 erodeFilter->SetBackgroundValue(m_BackgroundValue);
401 erodeFilter->SetErodeValue(m_BackgroundValue+1);
402 erodeFilter->Update();
403 working_image = erodeFilter->GetOutput();
404 StopCurrentStep<ImageType>(working_image);
406 //--------------------------------------------------------------------
407 //--------------------------------------------------------------------
408 // Step 5: resample to initial spacing
409 if (m_IntermediateSpacingFlag) {
410 StartNewStep("Resample to come back to the same sampling than input");
411 typedef clitk::ResampleImageWithOptionsFilter<ImageType> ResampleFilterType;
412 typename ResampleFilterType::Pointer resampleFilter = ResampleFilterType::New();
413 resampleFilter->SetDefaultPixelValue(m_BackgroundValue);
414 resampleFilter->SetInput(working_image);
415 resampleFilter->SetOutputSpacing(input->GetSpacing());
416 resampleFilter->SetGaussianFilteringEnabled(false);
417 // resampleFilter->SetVerboseOptions(true);
418 resampleFilter->Update();
419 working_image = resampleFilter->GetOutput();
420 StopCurrentStep<ImageType>(working_image);
423 //--------------------------------------------------------------------
424 //--------------------------------------------------------------------
425 // Pre Step 6: pad if not the same size : it can occur when downsample and upsample
426 //if (!HaveSameSizeAndSpacing(working_image, input)) {
427 if (working_image->GetLargestPossibleRegion() != input->GetLargestPossibleRegion()) {
428 StartNewStep("Pad to get the same size than input");
429 typename ImageType::Pointer temp = ImageType::New();
430 temp->CopyInformation(input);
431 temp->SetRegions(input->GetLargestPossibleRegion()); // Do not forget !!
433 temp->FillBuffer(m_BackgroundValue);
434 typename PasteFilterType::Pointer padFilter2 = PasteFilterType::New();
435 padFilter2->SetSourceImage(working_image);
436 padFilter2->SetDestinationImage(temp);
437 padFilter2->SetDestinationIndex(input->GetLargestPossibleRegion().GetIndex());
438 padFilter2->SetSourceRegion(working_image->GetLargestPossibleRegion());
439 padFilter2->Update();
440 working_image = padFilter2->GetOutput();
441 StopCurrentStep<ImageType>(working_image);
444 //DD("[debug] Rel Pos : no padding after");
447 //--------------------------------------------------------------------
448 //--------------------------------------------------------------------
449 // Step 6: combine input+thresholded relpos
450 StartNewStep("Combine with initial input (boolean And)");
451 typedef clitk::BooleanOperatorLabelImageFilter<ImageType> BoolFilterType;
452 typename BoolFilterType::Pointer combineFilter = BoolFilterType::New();
453 combineFilter->SetBackgroundValue(m_BackgroundValue);
454 combineFilter->SetBackgroundValue1(m_BackgroundValue);
455 combineFilter->SetBackgroundValue2(m_BackgroundValue);
456 combineFilter->SetForegroundValue(m_BackgroundValue+1);
457 combineFilter->SetInput1(input);
458 combineFilter->SetInput2(working_image);
459 if (GetInverseOrientationFlag())
460 combineFilter->SetOperationType(BoolFilterType::AndNot);
462 if (GetCombineWithOrFlag())
463 combineFilter->SetOperationType(BoolFilterType::Or);
465 combineFilter->SetOperationType(BoolFilterType::And);
467 combineFilter->InPlaceOff(); // Do not modify initial input (!)
468 combineFilter->Update();
469 working_image = combineFilter->GetOutput();
471 // Remove (if needed the object from the support)
472 if (GetRemoveObjectFlag()) {
473 combineFilter = BoolFilterType::New();
474 combineFilter->SetInput1(working_image);
475 combineFilter->SetInput2(object);
476 combineFilter->SetOperationType(BoolFilterType::AndNot);
477 combineFilter->InPlaceOn();
478 combineFilter->Update();
479 working_image = combineFilter->GetOutput();
481 // In the other case, we must *add* the initial object to keep it
482 // but not more than the initial support
484 combineFilter = BoolFilterType::New();
485 combineFilter->SetInput1(working_image);
486 combineFilter->SetInput2(object);
487 combineFilter->SetOperationType(BoolFilterType::Or);
488 combineFilter->InPlaceOn();
489 combineFilter->Update();
490 working_image = combineFilter->GetOutput(); // not needed because InPlaceOn ?
491 combineFilter = BoolFilterType::New();
492 combineFilter->SetInput1(working_image);
493 combineFilter->SetInput2(input);
494 combineFilter->SetOperationType(BoolFilterType::And);
495 combineFilter->InPlaceOn();
496 combineFilter->Update();
497 working_image = combineFilter->GetOutput();
500 StopCurrentStep<ImageType>(working_image);
502 //--------------------------------------------------------------------
503 //--------------------------------------------------------------------
505 if (GetAutoCropFlag()) {
506 StartNewStep("Final AutoCrop");
507 typedef clitk::AutoCropFilter<ImageType> CropFilterType;
508 typename CropFilterType::Pointer cropFilter = CropFilterType::New();
509 cropFilter->SetInput(working_image);
510 cropFilter->ReleaseDataFlagOff();
511 cropFilter->Update();
512 working_image = cropFilter->GetOutput();
513 StopCurrentStep<ImageType>(working_image);
516 //--------------------------------------------------------------------
517 //--------------------------------------------------------------------
519 // Final Step -> set output
520 this->SetNthOutput(0, working_image);
521 // this->GraftOutput(working_image);
523 //--------------------------------------------------------------------