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://www.centreleonberard.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"
25 #include "clitkCropLikeImageFilter.h"
29 #include <itkStatisticsLabelMapFilter.h>
30 #include <itkLabelImageToStatisticsLabelMapFilter.h>
31 #include <itkRegionOfInterestImageFilter.h>
32 #include <itkBinaryThresholdImageFilter.h>
33 #include <itkBinaryErodeImageFilter.h>
34 #include <itkBinaryBallStructuringElement.h>
35 #include <itkAddImageFilter.h>
36 #if ITK_VERSION_MAJOR >= 4
37 #include <itkDivideImageFilter.h>
39 #include <itkDivideByConstantImageFilter.h>
43 #include "RelativePositionPropImageFilter.h"
45 //--------------------------------------------------------------------
46 template <class ImageType>
47 clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
48 AddRelativePositionConstraintToLabelImageFilter():
50 itk::ImageToImageFilter<ImageType, ImageType>()
52 this->SetNumberOfRequiredInputs(2);
53 SetFuzzyThreshold(0.6);
54 SetBackgroundValue(0);
55 SetObjectBackgroundValue(0);
56 ClearOrientationType();
57 IntermediateSpacingFlagOn();
58 SetIntermediateSpacing(10);
60 InverseOrientationFlagOff();
62 CombineWithOrFlagOff();
65 FuzzyMapOnlyFlagOff();
69 //--------------------------------------------------------------------
72 //--------------------------------------------------------------------
73 template <class ImageType>
75 clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
76 SetInput(const ImageType * image)
78 // Process object is not const-correct so the const casting is required.
79 this->SetNthInput(0, const_cast<ImageType *>(image));
81 //--------------------------------------------------------------------
84 //--------------------------------------------------------------------
85 template <class ImageType>
87 clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
88 SetInputObject(const ImageType * image)
90 // Process object is not const-correct so the const casting is required.
91 this->SetNthInput(1, const_cast<ImageType *>(image));
93 //--------------------------------------------------------------------
96 //--------------------------------------------------------------------
97 template <class ImageType>
99 clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
100 ClearOrientationType()
102 m_OrientationTypeString.clear();
103 m_OrientationType.clear();
107 //--------------------------------------------------------------------
110 //--------------------------------------------------------------------
111 template <class ImageType>
113 clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
116 return m_OrientationType.size();
118 //--------------------------------------------------------------------
121 //--------------------------------------------------------------------
122 template <class ImageType>
124 clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
125 AddOrientationTypeString(std::string t)
127 m_OrientationTypeString.push_back(t);
129 if (t == "LeftTo") { AddOrientationType(LeftTo); return; }
130 if (t == "RightTo") { AddOrientationType(RightTo); return; }
131 if (t == "AntTo") { AddOrientationType(AntTo); return; }
132 if (t == "PostTo") { AddOrientationType(PostTo); return; }
133 if (t == "SupTo") { AddOrientationType(SupTo); return; }
134 if (t == "InfTo") { AddOrientationType(InfTo); return; }
136 if (t == "NotLeftTo") { AddOrientationType(LeftTo); InverseOrientationFlagOn(); return; }
137 if (t == "NotRightTo") { AddOrientationType(RightTo); InverseOrientationFlagOn(); return; }
138 if (t == "NotAntTo") { AddOrientationType(AntTo); InverseOrientationFlagOn(); return; }
139 if (t == "NotPostTo") { AddOrientationType(PostTo); InverseOrientationFlagOn(); return; }
140 if (t == "NotSupTo") { AddOrientationType(SupTo); InverseOrientationFlagOn(); return; }
141 if (t == "NotInfTo") { AddOrientationType(InfTo); InverseOrientationFlagOn(); return; }
143 if (t == "Angle") return;
145 clitkExceptionMacro("Error, you must provide LeftTo,RightTo or AntTo,PostTo or SupTo,InfTo (or NotLeftTo, NotRightTo etc) but you give " << t);
147 //--------------------------------------------------------------------
150 //--------------------------------------------------------------------
151 template <class ImageType>
153 clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
154 GenerateOutputInformation()
156 ImagePointer input = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(0));
157 ImagePointer outputImage = this->GetOutput(0);
158 outputImage->SetRegions(outputImage->GetLargestPossibleRegion());
160 //--------------------------------------------------------------------
163 //--------------------------------------------------------------------
164 template <class ImageType>
166 clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
167 GenerateInputRequestedRegion()
170 itk::ImageToImageFilter<ImageType, ImageType>::GenerateInputRequestedRegion();
171 // Get input pointers and set requested region to common region
172 ImagePointer input1 = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(0));
173 ImagePointer input2 = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(1));
174 input1->SetRequestedRegion(input1->GetLargestPossibleRegion());
175 input2->SetRequestedRegion(input2->GetLargestPossibleRegion());
177 //--------------------------------------------------------------------
180 //--------------------------------------------------------------------
181 template <class ImageType>
183 clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
184 AddAnglesInRad(double a, double b)
186 m_OrientationTypeString.push_back("Angle");
187 m_OrientationType.push_back(Angle);
188 m_Angle1.push_back(a);
189 m_Angle2.push_back(b);
191 //--------------------------------------------------------------------
194 //--------------------------------------------------------------------
195 template <class ImageType>
197 clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
198 AddAnglesInDeg(double a, double b)
200 AddAnglesInRad(clitk::deg2rad(a), clitk::deg2rad(b));
202 //--------------------------------------------------------------------
205 //--------------------------------------------------------------------
206 template <class ImageType>
208 clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
209 AddOrientationType(OrientationTypeEnumeration orientation)
211 m_OrientationType.push_back(orientation);
212 switch (orientation) {
214 m_Angle1.push_back(clitk::deg2rad(0));
215 m_Angle2.push_back(clitk::deg2rad(0));
218 m_Angle1.push_back(clitk::deg2rad(180));
219 m_Angle2.push_back(clitk::deg2rad(0));
222 m_Angle1.push_back(clitk::deg2rad(90));
223 m_Angle2.push_back(clitk::deg2rad(0));
226 m_Angle1.push_back(clitk::deg2rad(-90));
227 m_Angle2.push_back(clitk::deg2rad(0));
230 m_Angle1.push_back(clitk::deg2rad(0));
231 m_Angle2.push_back(clitk::deg2rad(90));
234 m_Angle1.push_back(clitk::deg2rad(0));
235 m_Angle2.push_back(clitk::deg2rad(-90));
248 //--------------------------------------------------------------------
251 //--------------------------------------------------------------------
252 template <class ImageType>
254 clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
257 DD((int)this->GetBackgroundValue());
258 DD((int)this->GetObjectBackgroundValue());
259 DDV(this->GetOrientationTypeString(), (uint)this->GetNumberOfAngles());
260 DD(this->GetIntermediateSpacingFlag());
261 DD(this->GetIntermediateSpacing());
262 DD(this->GetFuzzyThreshold());
263 DD(this->GetAutoCropFlag());
264 DD(this->GetInverseOrientationFlag());
265 DD(this->GetRemoveObjectFlag());
266 DD(this->GetCombineWithOrFlag());
268 //--------------------------------------------------------------------
271 //--------------------------------------------------------------------
272 template <class ImageType>
274 clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
277 if (GetNumberOfAngles() <1) {
278 clitkExceptionMacro("Add at least one orientation type");
281 if (GetVerboseOptionFlag()) {
282 for(int i=0; i<GetNumberOfAngles(); i++) {
283 std::cout << "Orientation \t:" << GetOrientationTypeString(i) << std::endl;
285 std::cout << "Interm Spacing \t:" << GetIntermediateSpacingFlag() << " " << GetIntermediateSpacing() << "mm" << std::endl;
286 std::cout << "Fuzzy threshold\t:" << GetFuzzyThreshold() << std::endl;
287 std::cout << "AutoCrop \t:" << GetAutoCropFlag() << std::endl;
288 std::cout << "InverseOrient \t:" << GetInverseOrientationFlag() << std::endl;
289 std::cout << "RemoveObject \t:" << GetRemoveObjectFlag() << std::endl;
290 std::cout << "CombineWithOr \t:" << GetCombineWithOrFlag() << std::endl;
294 input = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(0));
295 object = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(1));
296 static const unsigned int dim = ImageType::ImageDimension;
298 // Step 2: object pad to input image -> we want to compute the
299 // relative position for each point belonging to the input image
300 // domain, so we have to extend (pad) the object image to fit the
302 working_image = object;
303 if (!clitk::HaveSameSizeAndSpacing<ImageType, ImageType>(input, working_image)) {
304 StartNewStep("Pad (resize) object to input size");
306 if (0) { // OLD VERSION (TO REMOVE)
307 StartNewStep("Pad object to image size");
308 typename ImageType::Pointer output = ImageType::New();
310 for(unsigned int i=0; i<dim; i++) {
311 size[i] = lrint((input->GetLargestPossibleRegion().GetSize()[i]*
312 input->GetSpacing()[i])/(double)working_image->GetSpacing()[i]);
315 // The index of the input is not necessarily zero, so we have to
316 // take it into account (not done)
318 IndexType index = input->GetLargestPossibleRegion().GetIndex();
319 region.SetSize(size);
320 for(unsigned int i=0; i<dim; i++) {
322 std::cerr << "Index diff from zero : " << index << ". not done yet !" << std::endl;
326 // output->SetLargestPossibleRegion(region);
327 output->SetRegions(region);
328 output->SetSpacing(working_image->GetSpacing());
329 PointType origin = input->GetOrigin();
330 for(unsigned int i=0; i<dim; i++) {
331 origin[i] = index[i]*input->GetSpacing()[i] + input->GetOrigin()[i];
333 output->SetOrigin(origin);
334 // output->SetOrigin(input->GetOrigin());
337 output->FillBuffer(m_BackgroundValue);
338 typename PasteFilterType::Pointer padFilter = PasteFilterType::New();
339 // typename PasteFilterType::InputImageIndexType index;
340 for(unsigned int i=0; i<dim; i++) {
341 index[i] = -index[i]*input->GetSpacing()[i]/(double)working_image->GetSpacing()[i]
342 + lrint((working_image->GetOrigin()[i] - input->GetOrigin()[i])/working_image->GetSpacing()[i]);
344 padFilter->SetSourceImage(working_image);
345 padFilter->SetDestinationImage(output);
346 padFilter->SetDestinationIndex(index);
347 padFilter->SetSourceRegion(working_image->GetLargestPossibleRegion());
349 working_image = padFilter->GetOutput();
352 // Resize object like input
353 working_image = clitk::ResizeImageLike<ImageType>(working_image, input, GetBackgroundValue());
354 StopCurrentStep<ImageType>(working_image);
357 //--------------------------------------------------------------------
358 //--------------------------------------------------------------------
360 if (m_IntermediateSpacingFlag) {
361 StartNewStep("Resample object to intermediate spacing");
362 typedef clitk::ResampleImageWithOptionsFilter<ImageType> ResampleFilterType;
363 typename ResampleFilterType::Pointer resampleFilter = ResampleFilterType::New();
364 resampleFilter->SetInput(working_image);
365 resampleFilter->SetDefaultPixelValue(0);
366 resampleFilter->SetOutputIsoSpacing(m_IntermediateSpacing);
367 resampleFilter->SetGaussianFilteringEnabled(false);
368 // resampleFilter->SetVerboseOptions(true);
369 resampleFilter->Update();
370 working_image = resampleFilter->GetOutput();
371 StopCurrentStep<ImageType>(working_image);
374 // Keep object image (with resampline and pad)
375 object_resampled = working_image;
377 // Step 3: compute rel pos in object
378 StartNewStep("Relative Position Map");
379 typedef itk::RelativePositionPropImageFilter<ImageType, FloatImageType> RelPosFilterType;
380 typename RelPosFilterType::Pointer relPosFilter;
382 for(int i=0; i<GetNumberOfAngles(); i++) {
384 relPosFilter = RelPosFilterType::New();
385 relPosFilter->SetFast(GetFastFlag());
386 relPosFilter->SetRadius(GetRadius());
387 relPosFilter->SetInput(working_image);
388 relPosFilter->SetAlpha1(m_Angle1[i]); // xy plane
389 relPosFilter->SetAlpha2(m_Angle2[i]);
390 relPosFilter->SetK1(M_PI/2.0); // Opening parameter, default = pi/2
392 // relPosFilter->SetFast(true);
393 // relPosFilter->SetRadius(1); // seems sufficient in this case
395 // relPosFilter->SetVerboseProgress(true);
396 relPosFilter->Update();
397 relPos = relPosFilter->GetOutput();
399 if (GetNumberOfAngles() != 1) {
400 // Creation of the first m_FuzzyMap
402 m_FuzzyMap = clitk::NewImageLike<FloatImageType>(relPos, true);
403 m_FuzzyMap->FillBuffer(0.0);
406 // Add to current fuzzy map
407 typedef itk::AddImageFilter<FloatImageType, FloatImageType, FloatImageType> AddImageFilter;
408 typename AddImageFilter::Pointer addFilter = AddImageFilter::New();
409 addFilter->SetInput1(m_FuzzyMap);
410 addFilter->SetInput2(relPos);
412 m_FuzzyMap = addFilter->GetOutput();
414 else m_FuzzyMap = relPos;
417 // Divide by the number of relpos
418 if (GetNumberOfAngles() != 1) {
419 #if ITK_VERSION_MAJOR >= 4
420 typedef itk::DivideImageFilter<FloatImageType, FloatImageType, FloatImageType> DivideFilter;
421 typename DivideFilter::Pointer divideFilter = DivideFilter::New();
422 divideFilter->SetConstant2(GetNumberOfAngles());
424 typedef itk::DivideByConstantImageFilter<FloatImageType, float, FloatImageType> DivideFilter;
425 typename DivideFilter::Pointer divideFilter = DivideFilter::New();
426 divideFilter->SetConstant(GetNumberOfAngles());
428 divideFilter->SetInput(m_FuzzyMap);
429 divideFilter->Update();
430 m_FuzzyMap = divideFilter->GetOutput();
434 StopCurrentStep<FloatImageType>(relPos);
435 if (GetFuzzyMapOnlyFlag()) return;
437 //--------------------------------------------------------------------
438 //--------------------------------------------------------------------
439 StartNewStep("Map Threshold");
441 typedef itk::BinaryThresholdImageFilter<FloatImageType, ImageType> BinaryThresholdImageFilterType;
442 typename BinaryThresholdImageFilterType::Pointer thresholdFilter = BinaryThresholdImageFilterType::New();
443 thresholdFilter->SetInput(relPos);
444 thresholdFilter->SetOutsideValue(m_BackgroundValue);
445 thresholdFilter->SetInsideValue(m_BackgroundValue+1);
446 thresholdFilter->SetLowerThreshold(m_FuzzyThreshold);
447 thresholdFilter->Update();
448 working_image = thresholdFilter->GetOutput();
449 StopCurrentStep<ImageType>(working_image);
451 //--------------------------------------------------------------------
452 //--------------------------------------------------------------------
453 StartNewStep("Post Processing: erosion with initial mask");
454 // Step 2 : erosion with initial mask to exclude pixels that were
455 // inside the resampled version and outside the original mask
456 typedef itk::BinaryBallStructuringElement<unsigned int, ImageDimension> StructuringElementType;
457 StructuringElementType kernel;
459 kernel.CreateStructuringElement();
460 typedef itk::BinaryErodeImageFilter<ImageType, ImageType, StructuringElementType> ErodeFilterType;
461 typename ErodeFilterType::Pointer erodeFilter = ErodeFilterType::New();
462 erodeFilter->SetInput(working_image);
463 erodeFilter->SetKernel(kernel);
464 erodeFilter->SetBackgroundValue(m_BackgroundValue);
465 erodeFilter->SetErodeValue(m_BackgroundValue+1);
466 erodeFilter->Update();
467 working_image = erodeFilter->GetOutput();
468 StopCurrentStep<ImageType>(working_image);
470 //--------------------------------------------------------------------
471 //--------------------------------------------------------------------
472 // Step 5: resample to initial spacing
473 if (m_IntermediateSpacingFlag) {
474 StartNewStep("Resample to come back to the same sampling than input");
475 typedef clitk::ResampleImageWithOptionsFilter<ImageType> ResampleFilterType;
476 typename ResampleFilterType::Pointer resampleFilter = ResampleFilterType::New();
477 resampleFilter->SetDefaultPixelValue(m_BackgroundValue);
478 resampleFilter->SetInput(working_image);
479 resampleFilter->SetOutputSpacing(input->GetSpacing());
480 resampleFilter->SetGaussianFilteringEnabled(false);
481 // resampleFilter->SetVerboseOptions(true);
482 resampleFilter->Update();
483 working_image = resampleFilter->GetOutput();
484 StopCurrentStep<ImageType>(working_image);
487 //--------------------------------------------------------------------
488 //--------------------------------------------------------------------
489 // Pre Step 6: pad if not the same size : it can occur when downsample and upsample
490 //if (!HaveSameSizeAndSpacing(working_image, input)) {
491 if (working_image->GetLargestPossibleRegion() != input->GetLargestPossibleRegion()) {
492 StartNewStep("Pad to get the same size than input");
493 typename ImageType::Pointer temp = ImageType::New();
494 temp->CopyInformation(input);
495 temp->SetRegions(input->GetLargestPossibleRegion()); // Do not forget !!
497 temp->FillBuffer(m_BackgroundValue);
498 typename PasteFilterType::Pointer padFilter2 = PasteFilterType::New();
499 padFilter2->SetSourceImage(working_image);
500 padFilter2->SetDestinationImage(temp);
501 padFilter2->SetDestinationIndex(input->GetLargestPossibleRegion().GetIndex());
502 padFilter2->SetSourceRegion(working_image->GetLargestPossibleRegion());
503 padFilter2->Update();
504 working_image = padFilter2->GetOutput();
505 StopCurrentStep<ImageType>(working_image);
508 //DD("[debug] Rel Pos : no padding after");
511 //--------------------------------------------------------------------
512 //--------------------------------------------------------------------
513 // Step 6: combine input+thresholded relpos
514 StartNewStep("Combine with initial input (boolean And)");
515 typedef clitk::BooleanOperatorLabelImageFilter<ImageType> BoolFilterType;
516 typename BoolFilterType::Pointer combineFilter = BoolFilterType::New();
517 combineFilter->SetBackgroundValue(m_BackgroundValue);
518 combineFilter->SetBackgroundValue1(m_BackgroundValue);
519 combineFilter->SetBackgroundValue2(m_BackgroundValue);
520 combineFilter->SetForegroundValue(m_BackgroundValue+1);
521 combineFilter->SetInput1(input);
522 combineFilter->SetInput2(working_image);
523 if (GetInverseOrientationFlag())
524 combineFilter->SetOperationType(BoolFilterType::AndNot);
526 if (GetCombineWithOrFlag())
527 combineFilter->SetOperationType(BoolFilterType::Or);
529 combineFilter->SetOperationType(BoolFilterType::And);
531 combineFilter->InPlaceOff(); // Do not modify initial input (!)
532 combineFilter->Update();
533 working_image = combineFilter->GetOutput();
535 // Remove (if needed the object from the support)
536 if (GetRemoveObjectFlag()) {
537 combineFilter = BoolFilterType::New();
538 combineFilter->SetInput1(working_image);
539 combineFilter->SetInput2(object);
540 combineFilter->SetOperationType(BoolFilterType::AndNot);
541 combineFilter->InPlaceOn();
542 combineFilter->Update();
543 working_image = combineFilter->GetOutput();
545 // In the other case, we must *add* the initial object to keep it
546 // but not more than the initial support
548 combineFilter = BoolFilterType::New();
549 combineFilter->SetInput1(working_image);
550 combineFilter->SetInput2(object);
551 combineFilter->SetOperationType(BoolFilterType::Or);
552 combineFilter->InPlaceOn();
553 combineFilter->Update();
554 working_image = combineFilter->GetOutput(); // not needed because InPlaceOn ?
555 combineFilter = BoolFilterType::New();
556 combineFilter->SetInput1(working_image);
557 combineFilter->SetInput2(input);
558 combineFilter->SetOperationType(BoolFilterType::And);
559 combineFilter->InPlaceOn();
560 combineFilter->Update();
561 working_image = combineFilter->GetOutput();
564 StopCurrentStep<ImageType>(working_image);
566 //--------------------------------------------------------------------
567 //--------------------------------------------------------------------
569 if (GetAutoCropFlag()) {
570 StartNewStep("Final AutoCrop");
571 typedef clitk::AutoCropFilter<ImageType> CropFilterType;
572 typename CropFilterType::Pointer cropFilter = CropFilterType::New();
573 cropFilter->SetInput(working_image);
574 cropFilter->ReleaseDataFlagOff();
575 cropFilter->Update();
576 working_image = cropFilter->GetOutput();
577 StopCurrentStep<ImageType>(working_image);
580 //--------------------------------------------------------------------
581 //--------------------------------------------------------------------
583 // Final Step -> set output
584 this->SetNthOutput(0, working_image);
585 // this->GraftOutput(working_image);
587 //--------------------------------------------------------------------