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();
67 //--------------------------------------------------------------------
70 //--------------------------------------------------------------------
71 template <class ImageType>
73 clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
74 SetInput(const ImageType * image)
76 // Process object is not const-correct so the const casting is required.
77 this->SetNthInput(0, const_cast<ImageType *>(image));
79 //--------------------------------------------------------------------
82 //--------------------------------------------------------------------
83 template <class ImageType>
85 clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
86 SetInputObject(const ImageType * image)
88 // Process object is not const-correct so the const casting is required.
89 this->SetNthInput(1, const_cast<ImageType *>(image));
91 //--------------------------------------------------------------------
94 //--------------------------------------------------------------------
95 template <class ImageType>
97 clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
98 ClearOrientationType()
100 m_OrientationTypeString.clear();
101 m_OrientationType.clear();
105 //--------------------------------------------------------------------
108 //--------------------------------------------------------------------
109 template <class ImageType>
111 clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
114 return m_OrientationType.size();
116 //--------------------------------------------------------------------
119 //--------------------------------------------------------------------
120 template <class ImageType>
122 clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
123 AddOrientationTypeString(std::string t)
125 m_OrientationTypeString.push_back(t);
127 if (t == "LeftTo") { AddOrientationType(LeftTo); return; }
128 if (t == "RightTo") { AddOrientationType(RightTo); return; }
129 if (t == "AntTo") { AddOrientationType(AntTo); return; }
130 if (t == "PostTo") { AddOrientationType(PostTo); return; }
131 if (t == "SupTo") { AddOrientationType(SupTo); return; }
132 if (t == "InfTo") { AddOrientationType(InfTo); return; }
134 if (t == "NotLeftTo") { AddOrientationType(LeftTo); InverseOrientationFlagOn(); return; }
135 if (t == "NotRightTo") { AddOrientationType(RightTo); InverseOrientationFlagOn(); return; }
136 if (t == "NotAntTo") { AddOrientationType(AntTo); InverseOrientationFlagOn(); return; }
137 if (t == "NotPostTo") { AddOrientationType(PostTo); InverseOrientationFlagOn(); return; }
138 if (t == "NotSupTo") { AddOrientationType(SupTo); InverseOrientationFlagOn(); return; }
139 if (t == "NotInfTo") { AddOrientationType(InfTo); InverseOrientationFlagOn(); return; }
141 if (t == "Angle") return;
143 clitkExceptionMacro("Error, you must provide LeftTo,RightTo or AntTo,PostTo or SupTo,InfTo (or NotLeftTo, NotRightTo etc) but you give " << t);
145 //--------------------------------------------------------------------
148 //--------------------------------------------------------------------
149 template <class ImageType>
151 clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
152 GenerateOutputInformation()
154 ImagePointer input = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(0));
155 ImagePointer outputImage = this->GetOutput(0);
156 outputImage->SetRegions(outputImage->GetLargestPossibleRegion());
158 //--------------------------------------------------------------------
161 //--------------------------------------------------------------------
162 template <class ImageType>
164 clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
165 GenerateInputRequestedRegion()
168 itk::ImageToImageFilter<ImageType, ImageType>::GenerateInputRequestedRegion();
169 // Get input pointers and set requested region to common region
170 ImagePointer input1 = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(0));
171 ImagePointer input2 = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(1));
172 input1->SetRequestedRegion(input1->GetLargestPossibleRegion());
173 input2->SetRequestedRegion(input2->GetLargestPossibleRegion());
175 //--------------------------------------------------------------------
178 //--------------------------------------------------------------------
179 template <class ImageType>
181 clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
182 AddAngles(double a, double b)
184 m_OrientationTypeString.push_back("Angle");
185 m_OrientationType.push_back(Angle);
186 m_Angle1.push_back(a);
187 m_Angle2.push_back(b);
189 //--------------------------------------------------------------------
192 //--------------------------------------------------------------------
193 template <class ImageType>
195 clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
196 AddOrientationType(OrientationTypeEnumeration orientation)
198 m_OrientationType.push_back(orientation);
199 switch (orientation) {
201 m_Angle1.push_back(clitk::deg2rad(0));
202 m_Angle2.push_back(clitk::deg2rad(0));
205 m_Angle1.push_back(clitk::deg2rad(180));
206 m_Angle2.push_back(clitk::deg2rad(0));
209 m_Angle1.push_back(clitk::deg2rad(90));
210 m_Angle2.push_back(clitk::deg2rad(0));
213 m_Angle1.push_back(clitk::deg2rad(-90));
214 m_Angle2.push_back(clitk::deg2rad(0));
217 m_Angle1.push_back(clitk::deg2rad(0));
218 m_Angle2.push_back(clitk::deg2rad(90));
221 m_Angle1.push_back(clitk::deg2rad(0));
222 m_Angle2.push_back(clitk::deg2rad(-90));
235 //--------------------------------------------------------------------
238 //--------------------------------------------------------------------
239 template <class ImageType>
241 clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
244 DD((int)this->GetBackgroundValue());
245 DD((int)this->GetObjectBackgroundValue());
246 DDV(this->GetOrientationTypeString(), (uint)this->GetNumberOfAngles());
247 DD(this->GetIntermediateSpacingFlag());
248 DD(this->GetIntermediateSpacing());
249 DD(this->GetFuzzyThreshold());
250 DD(this->GetAutoCropFlag());
251 DD(this->GetInverseOrientationFlag());
252 DD(this->GetRemoveObjectFlag());
253 DD(this->GetCombineWithOrFlag());
255 //--------------------------------------------------------------------
258 //--------------------------------------------------------------------
259 template <class ImageType>
261 clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
264 if (GetNumberOfAngles() <1) {
265 clitkExceptionMacro("Add at least one orientation type");
268 if (GetVerboseOptionFlag()) {
269 for(int i=0; i<GetNumberOfAngles(); i++) {
270 std::cout << "Orientation \t:" << GetOrientationTypeString(i) << std::endl;
272 std::cout << "Interm Spacing \t:" << GetIntermediateSpacingFlag() << " " << GetIntermediateSpacing() << "mm" << std::endl;
273 std::cout << "Fuzzy threshold\t:" << GetFuzzyThreshold() << std::endl;
274 std::cout << "AutoCrop \t:" << GetAutoCropFlag() << std::endl;
275 std::cout << "InverseOrient \t:" << GetInverseOrientationFlag() << std::endl;
276 std::cout << "RemoveObject \t:" << GetRemoveObjectFlag() << std::endl;
277 std::cout << "CombineWithOr \t:" << GetCombineWithOrFlag() << std::endl;
281 input = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(0));
282 object = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(1));
283 static const unsigned int dim = ImageType::ImageDimension;
285 // Step 2: object pad to input image -> we want to compute the
286 // relative position for each point belonging to the input image
287 // domain, so we have to extend (pad) the object image to fit the
289 working_image = object;
290 if (!clitk::HaveSameSizeAndSpacing<ImageType, ImageType>(input, working_image)) {
291 StartNewStep("Pad (resize) object to input size");
293 if (0) { // OLD VERSION (TO REMOVE)
294 StartNewStep("Pad object to image size");
295 typename ImageType::Pointer output = ImageType::New();
297 for(unsigned int i=0; i<dim; i++) {
298 size[i] = lrint((input->GetLargestPossibleRegion().GetSize()[i]*
299 input->GetSpacing()[i])/(double)working_image->GetSpacing()[i]);
302 // The index of the input is not necessarily zero, so we have to
303 // take it into account (not done)
305 IndexType index = input->GetLargestPossibleRegion().GetIndex();
306 region.SetSize(size);
307 for(unsigned int i=0; i<dim; i++) {
309 std::cerr << "Index diff from zero : " << index << ". not done yet !" << std::endl;
313 // output->SetLargestPossibleRegion(region);
314 output->SetRegions(region);
315 output->SetSpacing(working_image->GetSpacing());
316 PointType origin = input->GetOrigin();
317 for(unsigned int i=0; i<dim; i++) {
318 origin[i] = index[i]*input->GetSpacing()[i] + input->GetOrigin()[i];
320 output->SetOrigin(origin);
321 // output->SetOrigin(input->GetOrigin());
324 output->FillBuffer(m_BackgroundValue);
325 typename PasteFilterType::Pointer padFilter = PasteFilterType::New();
326 // typename PasteFilterType::InputImageIndexType index;
327 for(unsigned int i=0; i<dim; i++) {
328 index[i] = -index[i]*input->GetSpacing()[i]/(double)working_image->GetSpacing()[i]
329 + lrint((working_image->GetOrigin()[i] - input->GetOrigin()[i])/working_image->GetSpacing()[i]);
331 padFilter->SetSourceImage(working_image);
332 padFilter->SetDestinationImage(output);
333 padFilter->SetDestinationIndex(index);
334 padFilter->SetSourceRegion(working_image->GetLargestPossibleRegion());
336 working_image = padFilter->GetOutput();
339 // Resize object like input
340 working_image = clitk::ResizeImageLike<ImageType>(working_image, input, GetBackgroundValue());
341 StopCurrentStep<ImageType>(working_image);
344 //--------------------------------------------------------------------
345 //--------------------------------------------------------------------
347 if (m_IntermediateSpacingFlag) {
348 StartNewStep("Resample object to intermediate spacing");
349 typedef clitk::ResampleImageWithOptionsFilter<ImageType> ResampleFilterType;
350 typename ResampleFilterType::Pointer resampleFilter = ResampleFilterType::New();
351 resampleFilter->SetInput(working_image);
352 resampleFilter->SetDefaultPixelValue(0);
353 resampleFilter->SetOutputIsoSpacing(m_IntermediateSpacing);
354 resampleFilter->SetGaussianFilteringEnabled(false);
355 // resampleFilter->SetVerboseOptions(true);
356 resampleFilter->Update();
357 working_image = resampleFilter->GetOutput();
358 StopCurrentStep<ImageType>(working_image);
361 // Keep object image (with resampline and pad)
362 object_resampled = working_image;
364 // Step 3: compute rel pos in object
365 StartNewStep("Relative Position Map");
366 typedef itk::RelativePositionPropImageFilter<ImageType, FloatImageType> RelPosFilterType;
367 typename RelPosFilterType::Pointer relPosFilter;
369 for(int i=0; i<GetNumberOfAngles(); i++) {
371 relPosFilter = RelPosFilterType::New();
372 relPosFilter->SetInput(working_image);
373 relPosFilter->SetAlpha1(m_Angle1[i]); // xy plane
374 relPosFilter->SetAlpha2(m_Angle2[i]);
375 relPosFilter->SetK1(M_PI/2.0); // Opening parameter, default = pi/2
376 relPosFilter->SetFast(true);
377 relPosFilter->SetRadius(1); // seems sufficient in this case
378 // relPosFilter->SetVerboseProgress(true);
379 relPosFilter->Update();
380 relPos = relPosFilter->GetOutput();
382 if (GetNumberOfAngles() != 1) {
383 // Creation of the first m_FuzzyMap
385 m_FuzzyMap = clitk::NewImageLike<FloatImageType>(relPos, true);
386 m_FuzzyMap->FillBuffer(0.0);
389 // Add to current fuzzy map
390 typedef itk::AddImageFilter<FloatImageType, FloatImageType, FloatImageType> AddImageFilter;
391 typename AddImageFilter::Pointer addFilter = AddImageFilter::New();
392 addFilter->SetInput1(m_FuzzyMap);
393 addFilter->SetInput2(relPos);
395 m_FuzzyMap = addFilter->GetOutput();
397 else m_FuzzyMap = relPos;
400 // Divide by the number of relpos
401 if (GetNumberOfAngles() != 1) {
402 #if ITK_VERSION_MAJOR >= 4
403 typedef itk::DivideImageFilter<FloatImageType, FloatImageType, FloatImageType> DivideFilter;
404 typename DivideFilter::Pointer divideFilter = DivideFilter::New();
405 divideFilter->SetConstant2(GetNumberOfAngles());
407 typedef itk::DivideByConstantImageFilter<FloatImageType, float, FloatImageType> DivideFilter;
408 typename DivideFilter::Pointer divideFilter = DivideFilter::New();
409 divideFilter->SetConstant(GetNumberOfAngles());
411 divideFilter->SetInput(m_FuzzyMap);
412 divideFilter->Update();
413 m_FuzzyMap = divideFilter->GetOutput();
417 StopCurrentStep<FloatImageType>(relPos);
418 if (GetFuzzyMapOnlyFlag()) return;
420 //--------------------------------------------------------------------
421 //--------------------------------------------------------------------
422 StartNewStep("Map Threshold");
424 typedef itk::BinaryThresholdImageFilter<FloatImageType, ImageType> BinaryThresholdImageFilterType;
425 typename BinaryThresholdImageFilterType::Pointer thresholdFilter = BinaryThresholdImageFilterType::New();
426 thresholdFilter->SetInput(relPos);
427 thresholdFilter->SetOutsideValue(m_BackgroundValue);
428 thresholdFilter->SetInsideValue(m_BackgroundValue+1);
429 thresholdFilter->SetLowerThreshold(m_FuzzyThreshold);
430 thresholdFilter->Update();
431 working_image = thresholdFilter->GetOutput();
432 StopCurrentStep<ImageType>(working_image);
434 //--------------------------------------------------------------------
435 //--------------------------------------------------------------------
436 StartNewStep("Post Processing: erosion with initial mask");
437 // Step 2 : erosion with initial mask to exclude pixels that were
438 // inside the resampled version and outside the original mask
439 typedef itk::BinaryBallStructuringElement<unsigned int, ImageDimension> StructuringElementType;
440 StructuringElementType kernel;
442 kernel.CreateStructuringElement();
443 typedef itk::BinaryErodeImageFilter<ImageType, ImageType, StructuringElementType> ErodeFilterType;
444 typename ErodeFilterType::Pointer erodeFilter = ErodeFilterType::New();
445 erodeFilter->SetInput(working_image);
446 erodeFilter->SetKernel(kernel);
447 erodeFilter->SetBackgroundValue(m_BackgroundValue);
448 erodeFilter->SetErodeValue(m_BackgroundValue+1);
449 erodeFilter->Update();
450 working_image = erodeFilter->GetOutput();
451 StopCurrentStep<ImageType>(working_image);
453 //--------------------------------------------------------------------
454 //--------------------------------------------------------------------
455 // Step 5: resample to initial spacing
456 if (m_IntermediateSpacingFlag) {
457 StartNewStep("Resample to come back to the same sampling than input");
458 typedef clitk::ResampleImageWithOptionsFilter<ImageType> ResampleFilterType;
459 typename ResampleFilterType::Pointer resampleFilter = ResampleFilterType::New();
460 resampleFilter->SetDefaultPixelValue(m_BackgroundValue);
461 resampleFilter->SetInput(working_image);
462 resampleFilter->SetOutputSpacing(input->GetSpacing());
463 resampleFilter->SetGaussianFilteringEnabled(false);
464 // resampleFilter->SetVerboseOptions(true);
465 resampleFilter->Update();
466 working_image = resampleFilter->GetOutput();
467 StopCurrentStep<ImageType>(working_image);
470 //--------------------------------------------------------------------
471 //--------------------------------------------------------------------
472 // Pre Step 6: pad if not the same size : it can occur when downsample and upsample
473 //if (!HaveSameSizeAndSpacing(working_image, input)) {
474 if (working_image->GetLargestPossibleRegion() != input->GetLargestPossibleRegion()) {
475 StartNewStep("Pad to get the same size than input");
476 typename ImageType::Pointer temp = ImageType::New();
477 temp->CopyInformation(input);
478 temp->SetRegions(input->GetLargestPossibleRegion()); // Do not forget !!
480 temp->FillBuffer(m_BackgroundValue);
481 typename PasteFilterType::Pointer padFilter2 = PasteFilterType::New();
482 padFilter2->SetSourceImage(working_image);
483 padFilter2->SetDestinationImage(temp);
484 padFilter2->SetDestinationIndex(input->GetLargestPossibleRegion().GetIndex());
485 padFilter2->SetSourceRegion(working_image->GetLargestPossibleRegion());
486 padFilter2->Update();
487 working_image = padFilter2->GetOutput();
488 StopCurrentStep<ImageType>(working_image);
491 //DD("[debug] Rel Pos : no padding after");
494 //--------------------------------------------------------------------
495 //--------------------------------------------------------------------
496 // Step 6: combine input+thresholded relpos
497 StartNewStep("Combine with initial input (boolean And)");
498 typedef clitk::BooleanOperatorLabelImageFilter<ImageType> BoolFilterType;
499 typename BoolFilterType::Pointer combineFilter = BoolFilterType::New();
500 combineFilter->SetBackgroundValue(m_BackgroundValue);
501 combineFilter->SetBackgroundValue1(m_BackgroundValue);
502 combineFilter->SetBackgroundValue2(m_BackgroundValue);
503 combineFilter->SetForegroundValue(m_BackgroundValue+1);
504 combineFilter->SetInput1(input);
505 combineFilter->SetInput2(working_image);
506 if (GetInverseOrientationFlag())
507 combineFilter->SetOperationType(BoolFilterType::AndNot);
509 if (GetCombineWithOrFlag())
510 combineFilter->SetOperationType(BoolFilterType::Or);
512 combineFilter->SetOperationType(BoolFilterType::And);
514 combineFilter->InPlaceOff(); // Do not modify initial input (!)
515 combineFilter->Update();
516 working_image = combineFilter->GetOutput();
518 // Remove (if needed the object from the support)
519 if (GetRemoveObjectFlag()) {
520 combineFilter = BoolFilterType::New();
521 combineFilter->SetInput1(working_image);
522 combineFilter->SetInput2(object);
523 combineFilter->SetOperationType(BoolFilterType::AndNot);
524 combineFilter->InPlaceOn();
525 combineFilter->Update();
526 working_image = combineFilter->GetOutput();
528 // In the other case, we must *add* the initial object to keep it
529 // but not more than the initial support
531 combineFilter = BoolFilterType::New();
532 combineFilter->SetInput1(working_image);
533 combineFilter->SetInput2(object);
534 combineFilter->SetOperationType(BoolFilterType::Or);
535 combineFilter->InPlaceOn();
536 combineFilter->Update();
537 working_image = combineFilter->GetOutput(); // not needed because InPlaceOn ?
538 combineFilter = BoolFilterType::New();
539 combineFilter->SetInput1(working_image);
540 combineFilter->SetInput2(input);
541 combineFilter->SetOperationType(BoolFilterType::And);
542 combineFilter->InPlaceOn();
543 combineFilter->Update();
544 working_image = combineFilter->GetOutput();
547 StopCurrentStep<ImageType>(working_image);
549 //--------------------------------------------------------------------
550 //--------------------------------------------------------------------
552 if (GetAutoCropFlag()) {
553 StartNewStep("Final AutoCrop");
554 typedef clitk::AutoCropFilter<ImageType> CropFilterType;
555 typename CropFilterType::Pointer cropFilter = CropFilterType::New();
556 cropFilter->SetInput(working_image);
557 cropFilter->ReleaseDataFlagOff();
558 cropFilter->Update();
559 working_image = cropFilter->GetOutput();
560 StopCurrentStep<ImageType>(working_image);
563 //--------------------------------------------------------------------
564 //--------------------------------------------------------------------
566 // Final Step -> set output
567 this->SetNthOutput(0, working_image);
568 // this->GraftOutput(working_image);
570 //--------------------------------------------------------------------