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 #include <itkDivideImageFilter.h>
39 #include "RelativePositionPropImageFilter.h"
41 //--------------------------------------------------------------------
42 template <class ImageType>
43 clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
44 AddRelativePositionConstraintToLabelImageFilter():
46 itk::ImageToImageFilter<ImageType, ImageType>()
48 this->SetNumberOfRequiredInputs(2);
49 SetFuzzyThreshold(0.6);
50 SetBackgroundValue(0);
51 SetObjectBackgroundValue(0);
52 ClearOrientationType();
53 IntermediateSpacingFlagOn();
54 SetIntermediateSpacing(10);
56 InverseOrientationFlagOff();
58 CombineWithOrFlagOff();
61 FuzzyMapOnlyFlagOff();
64 SetK1(std::acos(-1.0)/2);
66 //--------------------------------------------------------------------
69 //--------------------------------------------------------------------
70 template <class ImageType>
72 clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
73 SetInput(const ImageType * image)
75 // Process object is not const-correct so the const casting is required.
76 this->SetNthInput(0, const_cast<ImageType *>(image));
78 //--------------------------------------------------------------------
81 //--------------------------------------------------------------------
82 template <class ImageType>
84 clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
85 SetInputObject(const ImageType * image)
87 // Process object is not const-correct so the const casting is required.
88 this->SetNthInput(1, const_cast<ImageType *>(image));
90 //--------------------------------------------------------------------
93 //--------------------------------------------------------------------
94 template <class ImageType>
96 clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
97 ClearOrientationType()
99 m_OrientationTypeString.clear();
100 m_OrientationType.clear();
104 //--------------------------------------------------------------------
107 //--------------------------------------------------------------------
108 template <class ImageType>
110 clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
113 return m_OrientationType.size();
115 //--------------------------------------------------------------------
118 //--------------------------------------------------------------------
119 template <class ImageType>
121 clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
122 AddOrientationTypeString(std::string t)
124 m_OrientationTypeString.push_back(t);
126 if (t == "LeftTo") { AddOrientationType(LeftTo); return; }
127 if (t == "RightTo") { AddOrientationType(RightTo); return; }
128 if (t == "AntTo") { AddOrientationType(AntTo); return; }
129 if (t == "PostTo") { AddOrientationType(PostTo); return; }
130 if (t == "SupTo") { AddOrientationType(SupTo); return; }
131 if (t == "InfTo") { AddOrientationType(InfTo); return; }
133 if (t == "NotLeftTo") { AddOrientationType(LeftTo); InverseOrientationFlagOn(); return; }
134 if (t == "NotRightTo") { AddOrientationType(RightTo); InverseOrientationFlagOn(); return; }
135 if (t == "NotAntTo") { AddOrientationType(AntTo); InverseOrientationFlagOn(); return; }
136 if (t == "NotPostTo") { AddOrientationType(PostTo); InverseOrientationFlagOn(); return; }
137 if (t == "NotSupTo") { AddOrientationType(SupTo); InverseOrientationFlagOn(); return; }
138 if (t == "NotInfTo") { AddOrientationType(InfTo); InverseOrientationFlagOn(); return; }
140 if (t == "Angle") return;
142 clitkExceptionMacro("Error, you must provide LeftTo,RightTo or AntTo,PostTo or SupTo,InfTo (or NotLeftTo, NotRightTo etc) but you give " << t);
144 //--------------------------------------------------------------------
147 //--------------------------------------------------------------------
148 template <class ImageType>
150 clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
151 GenerateOutputInformation()
153 ImagePointer input = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(0));
154 ImagePointer outputImage = this->GetOutput(0);
155 outputImage->SetRegions(outputImage->GetLargestPossibleRegion());
157 //--------------------------------------------------------------------
160 //--------------------------------------------------------------------
161 template <class ImageType>
163 clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
164 GenerateInputRequestedRegion()
167 itk::ImageToImageFilter<ImageType, ImageType>::GenerateInputRequestedRegion();
168 // Get input pointers and set requested region to common region
169 ImagePointer input1 = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(0));
170 ImagePointer input2 = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(1));
171 input1->SetRequestedRegion(input1->GetLargestPossibleRegion());
172 input2->SetRequestedRegion(input2->GetLargestPossibleRegion());
174 //--------------------------------------------------------------------
177 //--------------------------------------------------------------------
178 template <class ImageType>
180 clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
181 AddAnglesInRad(double a, double b)
183 m_OrientationTypeString.push_back("Angle");
184 m_OrientationType.push_back(Angle);
185 m_Angle1.push_back(a);
186 m_Angle2.push_back(b);
188 //--------------------------------------------------------------------
191 //--------------------------------------------------------------------
192 template <class ImageType>
194 clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
195 AddAnglesInDeg(double a, double b)
197 AddAnglesInRad(clitk::deg2rad(a), clitk::deg2rad(b));
199 //--------------------------------------------------------------------
202 //--------------------------------------------------------------------
203 template <class ImageType>
205 clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
206 AddOrientationType(OrientationTypeEnumeration orientation)
208 m_OrientationType.push_back(orientation);
209 switch (orientation) {
211 m_Angle1.push_back(clitk::deg2rad(0));
212 m_Angle2.push_back(clitk::deg2rad(0));
215 m_Angle1.push_back(clitk::deg2rad(180));
216 m_Angle2.push_back(clitk::deg2rad(0));
219 m_Angle1.push_back(clitk::deg2rad(90));
220 m_Angle2.push_back(clitk::deg2rad(0));
223 m_Angle1.push_back(clitk::deg2rad(-90));
224 m_Angle2.push_back(clitk::deg2rad(0));
227 m_Angle1.push_back(clitk::deg2rad(0));
228 m_Angle2.push_back(clitk::deg2rad(90));
231 m_Angle1.push_back(clitk::deg2rad(0));
232 m_Angle2.push_back(clitk::deg2rad(-90));
245 //--------------------------------------------------------------------
248 //--------------------------------------------------------------------
249 template <class ImageType>
251 clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
254 DD((int)this->GetBackgroundValue());
255 DD((int)this->GetObjectBackgroundValue());
256 DDV(this->GetOrientationTypeString(), (uint)this->GetNumberOfAngles());
257 DD(this->GetIntermediateSpacingFlag());
258 DD(this->GetIntermediateSpacing());
259 DD(this->GetFuzzyThreshold());
260 DD(this->GetAutoCropFlag());
261 DD(this->GetInverseOrientationFlag());
262 DD(this->GetRemoveObjectFlag());
263 DD(this->GetCombineWithOrFlag());
265 //--------------------------------------------------------------------
268 //--------------------------------------------------------------------
269 template <class ImageType>
271 clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
274 if (GetNumberOfAngles() <1) {
275 clitkExceptionMacro("Add at least one orientation type");
278 if (GetVerboseOptionFlag()) {
279 for(int i=0; i<GetNumberOfAngles(); i++) {
280 std::cout << "Orientation \t:" << GetOrientationTypeString(i) << std::endl;
282 std::cout << "Interm Spacing \t:" << GetIntermediateSpacingFlag() << " " << GetIntermediateSpacing() << "mm" << std::endl;
283 std::cout << "Fuzzy threshold\t:" << GetFuzzyThreshold() << std::endl;
284 std::cout << "AutoCrop \t:" << GetAutoCropFlag() << std::endl;
285 std::cout << "InverseOrient \t:" << GetInverseOrientationFlag() << std::endl;
286 std::cout << "RemoveObject \t:" << GetRemoveObjectFlag() << std::endl;
287 std::cout << "CombineWithOr \t:" << GetCombineWithOrFlag() << std::endl;
291 input = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(0));
292 object = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(1));
293 static const unsigned int dim = ImageType::ImageDimension;
295 // Step 2: object pad to input image -> we want to compute the
296 // relative position for each point belonging to the input image
297 // domain, so we have to extend (pad) the object image to fit the
299 working_image = object;
300 if (!clitk::HaveSameSizeAndSpacing<ImageType, ImageType>(input, working_image)) {
301 StartNewStep("Pad (resize) object to input size");
303 if (0) { // OLD VERSION (TO REMOVE)
304 StartNewStep("Pad object to image size");
305 typename ImageType::Pointer output = ImageType::New();
307 for(unsigned int i=0; i<dim; i++) {
308 size[i] = lrint((input->GetLargestPossibleRegion().GetSize()[i]*
309 input->GetSpacing()[i])/(double)working_image->GetSpacing()[i]);
312 // The index of the input is not necessarily zero, so we have to
313 // take it into account (not done)
315 IndexType index = input->GetLargestPossibleRegion().GetIndex();
316 region.SetSize(size);
317 for(unsigned int i=0; i<dim; i++) {
319 std::cerr << "Index diff from zero : " << index << ". not done yet !" << std::endl;
323 // output->SetLargestPossibleRegion(region);
324 output->SetRegions(region);
325 output->SetSpacing(working_image->GetSpacing());
326 PointType origin = input->GetOrigin();
327 for(unsigned int i=0; i<dim; i++) {
328 origin[i] = index[i]*input->GetSpacing()[i] + input->GetOrigin()[i];
330 output->SetOrigin(origin);
331 // output->SetOrigin(input->GetOrigin());
334 output->FillBuffer(m_BackgroundValue);
335 typename PasteFilterType::Pointer padFilter = PasteFilterType::New();
336 // typename PasteFilterType::InputImageIndexType index;
337 for(unsigned int i=0; i<dim; i++) {
338 index[i] = -index[i]*input->GetSpacing()[i]/(double)working_image->GetSpacing()[i]
339 + lrint((working_image->GetOrigin()[i] - input->GetOrigin()[i])/working_image->GetSpacing()[i]);
341 padFilter->SetSourceImage(working_image);
342 padFilter->SetDestinationImage(output);
343 padFilter->SetDestinationIndex(index);
344 padFilter->SetSourceRegion(working_image->GetLargestPossibleRegion());
346 working_image = padFilter->GetOutput();
349 // Resize object like input
350 working_image = clitk::ResizeImageLike<ImageType>(working_image, input, GetBackgroundValue());
351 StopCurrentStep<ImageType>(working_image);
354 //--------------------------------------------------------------------
355 //--------------------------------------------------------------------
357 if (m_IntermediateSpacingFlag) {
358 StartNewStep("Resample object to intermediate spacing (" + toString(m_IntermediateSpacing) + ")");
359 typedef clitk::ResampleImageWithOptionsFilter<ImageType> ResampleFilterType;
360 typename ResampleFilterType::Pointer resampleFilter = ResampleFilterType::New();
361 resampleFilter->SetInput(working_image);
362 resampleFilter->SetDefaultPixelValue(0);
363 resampleFilter->SetOutputIsoSpacing(m_IntermediateSpacing);
364 resampleFilter->SetGaussianFilteringEnabled(false);
365 //resampleFilter->SetVerboseOptions(true);
366 resampleFilter->Update();
367 working_image = resampleFilter->GetOutput();
368 StopCurrentStep<ImageType>(working_image);
371 // Keep object image (with resampline and pad)
372 object_resampled = working_image;
374 // Step 3: compute rel pos in object
375 StartNewStep("Relative Position Map");
376 typedef itk::RelativePositionPropImageFilter<ImageType, FloatImageType> RelPosFilterType;
377 typename RelPosFilterType::Pointer relPosFilter;
379 for(int i=0; i<GetNumberOfAngles(); i++) {
381 relPosFilter = RelPosFilterType::New();
382 relPosFilter->SetFast(GetFastFlag());
383 relPosFilter->SetRadius(GetRadius());
384 relPosFilter->SetInput(working_image);
385 relPosFilter->SetAlpha1(m_Angle1[i]); // xy plane
386 relPosFilter->SetAlpha2(m_Angle2[i]);
387 relPosFilter->SetK1(GetK1());// M_PI/2.0); // Opening parameter, default = pi/2
388 // relPosFilter->SetVerboseProgress(true);
390 relPosFilter->Update();
391 relPos = relPosFilter->GetOutput();
393 if (GetNumberOfAngles() != 1) {
394 // Creation of the first m_FuzzyMap
396 m_FuzzyMap = clitk::NewImageLike<FloatImageType>(relPos, true);
397 m_FuzzyMap->FillBuffer(0.0);
400 // Add to current fuzzy map
401 typedef itk::AddImageFilter<FloatImageType, FloatImageType, FloatImageType> AddImageFilter;
402 typename AddImageFilter::Pointer addFilter = AddImageFilter::New();
403 addFilter->SetInput1(m_FuzzyMap);
404 addFilter->SetInput2(relPos);
406 m_FuzzyMap = addFilter->GetOutput();
408 else m_FuzzyMap = relPos;
411 // Divide by the number of relpos
412 if (GetNumberOfAngles() != 1) {
413 typedef itk::DivideImageFilter<FloatImageType, FloatImageType, FloatImageType> DivideFilter;
414 typename DivideFilter::Pointer divideFilter = DivideFilter::New();
415 divideFilter->SetConstant2(GetNumberOfAngles());
416 divideFilter->SetInput(m_FuzzyMap);
417 divideFilter->Update();
418 m_FuzzyMap = divideFilter->GetOutput();
422 StopCurrentStep<FloatImageType>(relPos);
423 if (GetFuzzyMapOnlyFlag()) {
424 // If the spacing is used, recompute fuzzy map with the input size/spacing
425 if (m_IntermediateSpacingFlag) {
426 StartNewStep("Resample FuzzyMap to come back to the same sampling than input");
427 typedef clitk::ResampleImageWithOptionsFilter<FloatImageType> ResampleFilterType;
428 typename ResampleFilterType::Pointer resampleFilter = ResampleFilterType::New();
429 resampleFilter->SetDefaultPixelValue(0.0); //Default fuzzymap value FIXME
430 resampleFilter->SetInput(relPos);
431 resampleFilter->SetOutputSpacing(input->GetSpacing());
432 resampleFilter->SetGaussianFilteringEnabled(false);
433 resampleFilter->Update();
434 relPos = m_FuzzyMap = resampleFilter->GetOutput();
435 StopCurrentStep<FloatImageType>(relPos);
437 // Need to put exactly the same size
438 if (relPos->GetLargestPossibleRegion() != input->GetLargestPossibleRegion()) {
439 StartNewStep("Pad to get the same size than input");
440 typename FloatImageType::Pointer temp = FloatImageType::New();
441 temp->CopyInformation(input);
442 temp->SetRegions(input->GetLargestPossibleRegion()); // Do not forget !!
444 temp->FillBuffer(0.0); // Default fuzzymap value FIXME
445 typename PasteFloatFilterType::Pointer padFilter2 = PasteFloatFilterType::New();
446 padFilter2->SetSourceImage(relPos);
447 padFilter2->SetDestinationImage(temp);
448 padFilter2->SetDestinationIndex(input->GetLargestPossibleRegion().GetIndex());
449 padFilter2->SetSourceRegion(relPos->GetLargestPossibleRegion());
450 padFilter2->Update();
451 relPos = padFilter2->GetOutput();
452 StopCurrentStep<FloatImageType>(relPos);
457 StopCurrentStep<FloatImageType>(relPos);
462 //--------------------------------------------------------------------
463 //--------------------------------------------------------------------
464 StartNewStep("Map Threshold");
466 typedef itk::BinaryThresholdImageFilter<FloatImageType, ImageType> BinaryThresholdImageFilterType;
467 typename BinaryThresholdImageFilterType::Pointer thresholdFilter = BinaryThresholdImageFilterType::New();
468 thresholdFilter->SetInput(relPos);
469 thresholdFilter->SetOutsideValue(m_BackgroundValue);
470 thresholdFilter->SetInsideValue(m_BackgroundValue+1);
471 thresholdFilter->SetLowerThreshold(m_FuzzyThreshold);
472 thresholdFilter->Update();
473 working_image = thresholdFilter->GetOutput();
474 StopCurrentStep<ImageType>(working_image);
476 //--------------------------------------------------------------------
477 //--------------------------------------------------------------------
478 StartNewStep("Post Processing: erosion with initial mask");
479 // Step 2 : erosion with initial mask to exclude pixels that were
480 // inside the resampled version and outside the original mask
481 typedef itk::BinaryBallStructuringElement<unsigned int, ImageDimension> StructuringElementType;
482 StructuringElementType kernel;
484 kernel.CreateStructuringElement();
485 typedef itk::BinaryErodeImageFilter<ImageType, ImageType, StructuringElementType> ErodeFilterType;
486 typename ErodeFilterType::Pointer erodeFilter = ErodeFilterType::New();
487 erodeFilter->SetInput(working_image);
488 erodeFilter->SetKernel(kernel);
489 erodeFilter->SetBackgroundValue(m_BackgroundValue);
490 erodeFilter->SetErodeValue(m_BackgroundValue+1);
491 erodeFilter->Update();
492 working_image = erodeFilter->GetOutput();
493 StopCurrentStep<ImageType>(working_image);
495 //--------------------------------------------------------------------
496 //--------------------------------------------------------------------
497 // Step 5: resample to initial spacing
498 if (m_IntermediateSpacingFlag) {
499 StartNewStep("Resample to come back to the same sampling than input");
500 typedef clitk::ResampleImageWithOptionsFilter<ImageType> ResampleFilterType;
501 typename ResampleFilterType::Pointer resampleFilter = ResampleFilterType::New();
502 resampleFilter->SetDefaultPixelValue(m_BackgroundValue);
503 resampleFilter->SetInput(working_image);
504 resampleFilter->SetOutputSpacing(input->GetSpacing());
505 resampleFilter->SetGaussianFilteringEnabled(false);
506 // resampleFilter->SetVerboseOptions(true);
507 resampleFilter->Update();
508 working_image = resampleFilter->GetOutput();
509 StopCurrentStep<ImageType>(working_image);
512 //--------------------------------------------------------------------
513 //--------------------------------------------------------------------
514 // Pre Step 6: pad if not the same size : it can occur when downsample and upsample
515 //if (!HaveSameSizeAndSpacing(working_image, input)) {
516 if (working_image->GetLargestPossibleRegion() != input->GetLargestPossibleRegion()) {
517 StartNewStep("Pad to get the same size than input");
518 typename ImageType::Pointer temp = ImageType::New();
519 temp->CopyInformation(input);
520 temp->SetRegions(input->GetLargestPossibleRegion()); // Do not forget !!
522 temp->FillBuffer(m_BackgroundValue);
523 typename PasteFilterType::Pointer padFilter2 = PasteFilterType::New();
524 padFilter2->SetSourceImage(working_image);
525 padFilter2->SetDestinationImage(temp);
526 padFilter2->SetDestinationIndex(input->GetLargestPossibleRegion().GetIndex());
527 padFilter2->SetSourceRegion(working_image->GetLargestPossibleRegion());
528 padFilter2->Update();
529 working_image = padFilter2->GetOutput();
530 StopCurrentStep<ImageType>(working_image);
533 //DD("[debug] Rel Pos : no padding after");
536 //--------------------------------------------------------------------
537 //--------------------------------------------------------------------
538 // Step 6: combine input+thresholded relpos
539 StartNewStep("Combine with initial input (boolean And)");
540 typedef clitk::BooleanOperatorLabelImageFilter<ImageType> BoolFilterType;
541 typename BoolFilterType::Pointer combineFilter = BoolFilterType::New();
542 combineFilter->SetBackgroundValue(m_BackgroundValue);
543 combineFilter->SetBackgroundValue1(m_BackgroundValue);
544 combineFilter->SetBackgroundValue2(m_BackgroundValue);
545 combineFilter->SetForegroundValue(m_BackgroundValue+1);
546 combineFilter->SetInput1(input);
547 combineFilter->SetInput2(working_image);
548 if (GetInverseOrientationFlag())
549 combineFilter->SetOperationType(BoolFilterType::AndNot);
551 if (GetCombineWithOrFlag())
552 combineFilter->SetOperationType(BoolFilterType::Or);
554 combineFilter->SetOperationType(BoolFilterType::And);
556 combineFilter->InPlaceOff(); // Do not modify initial input (!)
557 combineFilter->Update();
558 working_image = combineFilter->GetOutput();
560 // Remove (if needed the object from the support)
561 if (GetRemoveObjectFlag()) {
562 combineFilter = BoolFilterType::New();
563 combineFilter->SetInput1(working_image);
564 combineFilter->SetInput2(object);
565 combineFilter->SetOperationType(BoolFilterType::AndNot);
566 combineFilter->InPlaceOn();
567 combineFilter->Update();
568 working_image = combineFilter->GetOutput();
570 // In the other case, we must *add* the initial object to keep it
571 // but not more than the initial support
573 combineFilter = BoolFilterType::New();
574 combineFilter->SetInput1(working_image);
575 combineFilter->SetInput2(object);
576 combineFilter->SetOperationType(BoolFilterType::Or);
577 combineFilter->InPlaceOn();
578 combineFilter->Update();
579 working_image = combineFilter->GetOutput(); // not needed because InPlaceOn ?
580 combineFilter = BoolFilterType::New();
581 combineFilter->SetInput1(working_image);
582 combineFilter->SetInput2(input);
583 combineFilter->SetOperationType(BoolFilterType::And);
584 combineFilter->InPlaceOn();
585 combineFilter->Update();
586 working_image = combineFilter->GetOutput();
589 StopCurrentStep<ImageType>(working_image);
591 //--------------------------------------------------------------------
592 //--------------------------------------------------------------------
594 if (GetAutoCropFlag()) {
595 StartNewStep("Final AutoCrop");
596 typedef clitk::AutoCropFilter<ImageType> CropFilterType;
597 typename CropFilterType::Pointer cropFilter = CropFilterType::New();
598 cropFilter->SetInput(working_image);
599 cropFilter->ReleaseDataFlagOff();
600 cropFilter->Update();
601 working_image = cropFilter->GetOutput();
602 StopCurrentStep<ImageType>(working_image);
605 //--------------------------------------------------------------------
606 //--------------------------------------------------------------------
608 // Final Step -> set output
609 this->SetNthOutput(0, working_image);
610 // this->GraftOutput(working_image);
612 //--------------------------------------------------------------------