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"
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 #if ITK_VERSION_MAJOR >= 4
36 #include <itkDivideImageFilter.h>
38 #include <itkDivideByConstantImageFilter.h>
42 #include "RelativePositionPropImageFilter.h"
44 //--------------------------------------------------------------------
45 template <class ImageType>
46 clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
47 AddRelativePositionConstraintToLabelImageFilter():
49 itk::ImageToImageFilter<ImageType, ImageType>()
51 this->SetNumberOfRequiredInputs(2);
52 SetFuzzyThreshold(0.6);
53 SetBackgroundValue(0);
54 SetObjectBackgroundValue(0);
55 ClearOrientationType();
56 IntermediateSpacingFlagOn();
57 SetIntermediateSpacing(10);
59 InverseOrientationFlagOff();
61 CombineWithOrFlagOff();
65 //--------------------------------------------------------------------
68 //--------------------------------------------------------------------
69 template <class ImageType>
71 clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
72 SetInput(const ImageType * image)
74 // Process object is not const-correct so the const casting is required.
75 this->SetNthInput(0, const_cast<ImageType *>(image));
77 //--------------------------------------------------------------------
80 //--------------------------------------------------------------------
81 template <class ImageType>
83 clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
84 SetInputObject(const ImageType * image)
86 // Process object is not const-correct so the const casting is required.
87 this->SetNthInput(1, const_cast<ImageType *>(image));
89 //--------------------------------------------------------------------
92 //--------------------------------------------------------------------
93 template <class ImageType>
95 clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
96 ClearOrientationType()
98 m_OrientationTypeString.clear();
99 m_OrientationType.clear();
103 //--------------------------------------------------------------------
106 //--------------------------------------------------------------------
107 template <class ImageType>
109 clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
112 return m_OrientationType.size();
114 //--------------------------------------------------------------------
117 //--------------------------------------------------------------------
118 template <class ImageType>
120 clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
121 AddOrientationTypeString(std::string t)
123 m_OrientationTypeString.push_back(t);
125 if (t == "LeftTo") { AddOrientationType(LeftTo); return; }
126 if (t == "RightTo") { AddOrientationType(RightTo); return; }
127 if (t == "AntTo") { AddOrientationType(AntTo); return; }
128 if (t == "PostTo") { AddOrientationType(PostTo); return; }
129 if (t == "SupTo") { AddOrientationType(SupTo); return; }
130 if (t == "InfTo") { AddOrientationType(InfTo); return; }
132 if (t == "NotLeftTo") { AddOrientationType(LeftTo); InverseOrientationFlagOn(); return; }
133 if (t == "NotRightTo") { AddOrientationType(RightTo); InverseOrientationFlagOn(); return; }
134 if (t == "NotAntTo") { AddOrientationType(AntTo); InverseOrientationFlagOn(); return; }
135 if (t == "NotPostTo") { AddOrientationType(PostTo); InverseOrientationFlagOn(); return; }
136 if (t == "NotSupTo") { AddOrientationType(SupTo); InverseOrientationFlagOn(); return; }
137 if (t == "NotInfTo") { AddOrientationType(InfTo); InverseOrientationFlagOn(); return; }
139 clitkExceptionMacro("Error, you must provide LeftTo,RightTo or AntTo,PostTo or SupTo,InfTo (or NotLeftTo, NotRightTo etc) but you give " << t);
141 //--------------------------------------------------------------------
144 //--------------------------------------------------------------------
145 template <class ImageType>
147 clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
148 GenerateOutputInformation()
150 ImagePointer input = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(0));
151 ImagePointer outputImage = this->GetOutput(0);
152 outputImage->SetRegions(outputImage->GetLargestPossibleRegion());
154 //--------------------------------------------------------------------
157 //--------------------------------------------------------------------
158 template <class ImageType>
160 clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
161 GenerateInputRequestedRegion()
164 itk::ImageToImageFilter<ImageType, ImageType>::GenerateInputRequestedRegion();
165 // Get input pointers and set requested region to common region
166 ImagePointer input1 = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(0));
167 ImagePointer input2 = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(1));
168 input1->SetRequestedRegion(input1->GetLargestPossibleRegion());
169 input2->SetRequestedRegion(input2->GetLargestPossibleRegion());
171 //--------------------------------------------------------------------
174 //--------------------------------------------------------------------
175 template <class ImageType>
177 clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
178 AddAngles(double a, double b)
180 AddOrientationTypeString("Angle");
181 m_Angle1.push_back(a);
182 m_Angle2.push_back(b);
184 //--------------------------------------------------------------------
187 //--------------------------------------------------------------------
188 template <class ImageType>
190 clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
191 AddOrientationType(OrientationTypeEnumeration orientation)
193 m_OrientationType.push_back(orientation);
194 switch (orientation) {
196 m_Angle1.push_back(clitk::deg2rad(0));
197 m_Angle2.push_back(clitk::deg2rad(0));
200 m_Angle1.push_back(clitk::deg2rad(180));
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(-90));
209 m_Angle2.push_back(clitk::deg2rad(0));
212 m_Angle1.push_back(clitk::deg2rad(0));
213 m_Angle2.push_back(clitk::deg2rad(90));
216 m_Angle1.push_back(clitk::deg2rad(0));
217 m_Angle2.push_back(clitk::deg2rad(-90));
230 //--------------------------------------------------------------------
233 //--------------------------------------------------------------------
234 template <class ImageType>
236 clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
239 DD((int)this->GetBackgroundValue());
240 DD((int)this->GetObjectBackgroundValue());
241 DDV(this->GetOrientationTypeString(), (uint)this->GetNumberOfAngles());
242 DD(this->GetIntermediateSpacingFlag());
243 DD(this->GetIntermediateSpacing());
244 DD(this->GetFuzzyThreshold());
245 DD(this->GetAutoCropFlag());
246 DD(this->GetInverseOrientationFlag());
247 DD(this->GetRemoveObjectFlag());
248 DD(this->GetCombineWithOrFlag());
250 //--------------------------------------------------------------------
253 //--------------------------------------------------------------------
254 template <class ImageType>
256 clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
259 if (GetNumberOfAngles() <1) {
260 clitkExceptionMacro("Add at least one orientation type");
263 if (GetVerboseOptionFlag()) {
264 for(int i=0; i<GetNumberOfAngles(); i++) {
265 std::cout << "Orientation \t:" << GetOrientationTypeString(i) << std::endl;
267 std::cout << "Interm Spacing \t:" << GetIntermediateSpacingFlag() << " " << GetIntermediateSpacing() << "mm" << std::endl;
268 std::cout << "Fuzzy threshold\t:" << GetFuzzyThreshold() << std::endl;
269 std::cout << "AutoCrop \t:" << GetAutoCropFlag() << std::endl;
270 std::cout << "InverseOrient \t:" << GetInverseOrientationFlag() << std::endl;
271 std::cout << "RemoveObject \t:" << GetRemoveObjectFlag() << std::endl;
272 std::cout << "CombineWithOr \t:" << GetCombineWithOrFlag() << std::endl;
276 input = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(0));
277 object = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(1));
278 static const unsigned int dim = ImageType::ImageDimension;
280 // Step 2: object pad to input image -> we want to compute the
281 // relative position for each point belonging to the input image
282 // domain, so we have to extend (pad) the object image to fit the
284 working_image = object;
285 if (!clitk::HaveSameSizeAndSpacing<ImageType, ImageType>(input, working_image)) {
286 StartNewStep("Pad (resize) object to input size");
288 if (0) { // OLD VERSION (TO REMOVE)
289 StartNewStep("Pad object to image size");
290 typename ImageType::Pointer output = ImageType::New();
292 for(unsigned int i=0; i<dim; i++) {
293 size[i] = lrint((input->GetLargestPossibleRegion().GetSize()[i]*
294 input->GetSpacing()[i])/(double)working_image->GetSpacing()[i]);
297 // The index of the input is not necessarily zero, so we have to
298 // take it into account (not done)
300 IndexType index = input->GetLargestPossibleRegion().GetIndex();
301 region.SetSize(size);
302 for(unsigned int i=0; i<dim; i++) {
304 std::cerr << "Index diff from zero : " << index << ". not done yet !" << std::endl;
308 // output->SetLargestPossibleRegion(region);
309 output->SetRegions(region);
310 output->SetSpacing(working_image->GetSpacing());
311 PointType origin = input->GetOrigin();
312 for(unsigned int i=0; i<dim; i++) {
313 origin[i] = index[i]*input->GetSpacing()[i] + input->GetOrigin()[i];
315 output->SetOrigin(origin);
316 // output->SetOrigin(input->GetOrigin());
319 output->FillBuffer(m_BackgroundValue);
320 typename PasteFilterType::Pointer padFilter = PasteFilterType::New();
321 // typename PasteFilterType::InputImageIndexType index;
322 for(unsigned int i=0; i<dim; i++) {
323 index[i] = -index[i]*input->GetSpacing()[i]/(double)working_image->GetSpacing()[i]
324 + lrint((working_image->GetOrigin()[i] - input->GetOrigin()[i])/working_image->GetSpacing()[i]);
326 padFilter->SetSourceImage(working_image);
327 padFilter->SetDestinationImage(output);
328 padFilter->SetDestinationIndex(index);
329 padFilter->SetSourceRegion(working_image->GetLargestPossibleRegion());
331 working_image = padFilter->GetOutput();
334 // Resize object like input
335 working_image = clitk::ResizeImageLike<ImageType>(working_image, input, GetBackgroundValue());
336 StopCurrentStep<ImageType>(working_image);
339 //--------------------------------------------------------------------
340 //--------------------------------------------------------------------
342 if (m_IntermediateSpacingFlag) {
343 StartNewStep("Resample object to intermediate spacing");
344 typedef clitk::ResampleImageWithOptionsFilter<ImageType> ResampleFilterType;
345 typename ResampleFilterType::Pointer resampleFilter = ResampleFilterType::New();
346 resampleFilter->SetInput(working_image);
347 resampleFilter->SetDefaultPixelValue(0);
348 resampleFilter->SetOutputIsoSpacing(m_IntermediateSpacing);
349 resampleFilter->SetGaussianFilteringEnabled(false);
350 // resampleFilter->SetVerboseOptions(true);
351 resampleFilter->Update();
352 working_image = resampleFilter->GetOutput();
353 StopCurrentStep<ImageType>(working_image);
356 // Keep object image (with resampline and pad)
357 object_resampled = working_image;
359 // Step 3: compute rel pos in object
360 StartNewStep("Relative Position Map");
361 typedef itk::RelativePositionPropImageFilter<ImageType, FloatImageType> RelPosFilterType;
362 typename RelPosFilterType::Pointer relPosFilter;
364 typename FloatImageType::Pointer m_FuzzyMap;
365 for(int i=0; i<GetNumberOfAngles(); i++) {
367 relPosFilter = RelPosFilterType::New();
368 relPosFilter->SetInput(working_image);
369 relPosFilter->SetAlpha1(m_Angle1[i]); // xy plane
370 relPosFilter->SetAlpha2(m_Angle2[i]);
371 relPosFilter->SetK1(M_PI/2.0); // Opening parameter, default = pi/2
372 relPosFilter->SetFast(true);
373 relPosFilter->SetRadius(1); // seems sufficient in this case
374 // relPosFilter->SetVerboseProgress(true);
375 relPosFilter->Update();
376 relPos = relPosFilter->GetOutput();
378 if (GetNumberOfAngles() != 1) {
379 // Creation of the first m_FuzzyMap
381 m_FuzzyMap = clitk::NewImageLike<FloatImageType>(relPos, true);
382 m_FuzzyMap->FillBuffer(0.0);
385 // Add to current fuzzy map
386 typedef itk::AddImageFilter<FloatImageType, FloatImageType, FloatImageType> AddImageFilter;
387 typename AddImageFilter::Pointer addFilter = AddImageFilter::New();
388 addFilter->SetInput1(m_FuzzyMap);
389 addFilter->SetInput2(relPos);
391 m_FuzzyMap = addFilter->GetOutput();
393 else m_FuzzyMap = relPos;
396 // Divide by the number of relpos
397 if (GetNumberOfAngles() != 1) {
398 #if ITK_VERSION_MAJOR >= 4
399 typedef itk::DivideImageFilter<FloatImageType, FloatImageType, FloatImageType> DivideFilter;
400 typename DivideFilter::Pointer divideFilter = DivideFilter::New();
401 divideFilter->SetConstant2(GetNumberOfAngles());
403 typedef itk::DivideByConstantImageFilter<FloatImageType, float, FloatImageType> DivideFilter;
404 typename DivideFilter::Pointer divideFilter = DivideFilter::New();
405 divideFilter->SetConstant(GetNumberOfAngles());
407 divideFilter->SetInput(m_FuzzyMap);
408 divideFilter->Update();
409 m_FuzzyMap = divideFilter->GetOutput();
413 StopCurrentStep<FloatImageType>(relPos);
415 //--------------------------------------------------------------------
416 //--------------------------------------------------------------------
417 StartNewStep("Map Threshold");
419 typedef itk::BinaryThresholdImageFilter<FloatImageType, ImageType> BinaryThresholdImageFilterType;
420 typename BinaryThresholdImageFilterType::Pointer thresholdFilter = BinaryThresholdImageFilterType::New();
421 thresholdFilter->SetInput(relPos);
422 thresholdFilter->SetOutsideValue(m_BackgroundValue);
423 thresholdFilter->SetInsideValue(m_BackgroundValue+1);
424 thresholdFilter->SetLowerThreshold(m_FuzzyThreshold);
425 thresholdFilter->Update();
426 working_image = thresholdFilter->GetOutput();
427 StopCurrentStep<ImageType>(working_image);
429 //--------------------------------------------------------------------
430 //--------------------------------------------------------------------
431 StartNewStep("Post Processing: erosion with initial mask");
432 // Step 2 : erosion with initial mask to exclude pixels that were
433 // inside the resampled version and outside the original mask
434 typedef itk::BinaryBallStructuringElement<unsigned int, ImageDimension> StructuringElementType;
435 StructuringElementType kernel;
437 kernel.CreateStructuringElement();
438 typedef itk::BinaryErodeImageFilter<ImageType, ImageType, StructuringElementType> ErodeFilterType;
439 typename ErodeFilterType::Pointer erodeFilter = ErodeFilterType::New();
440 erodeFilter->SetInput(working_image);
441 erodeFilter->SetKernel(kernel);
442 erodeFilter->SetBackgroundValue(m_BackgroundValue);
443 erodeFilter->SetErodeValue(m_BackgroundValue+1);
444 erodeFilter->Update();
445 working_image = erodeFilter->GetOutput();
446 StopCurrentStep<ImageType>(working_image);
448 //--------------------------------------------------------------------
449 //--------------------------------------------------------------------
450 // Step 5: resample to initial spacing
451 if (m_IntermediateSpacingFlag) {
452 StartNewStep("Resample to come back to the same sampling than input");
453 typedef clitk::ResampleImageWithOptionsFilter<ImageType> ResampleFilterType;
454 typename ResampleFilterType::Pointer resampleFilter = ResampleFilterType::New();
455 resampleFilter->SetDefaultPixelValue(m_BackgroundValue);
456 resampleFilter->SetInput(working_image);
457 resampleFilter->SetOutputSpacing(input->GetSpacing());
458 resampleFilter->SetGaussianFilteringEnabled(false);
459 // resampleFilter->SetVerboseOptions(true);
460 resampleFilter->Update();
461 working_image = resampleFilter->GetOutput();
462 StopCurrentStep<ImageType>(working_image);
465 //--------------------------------------------------------------------
466 //--------------------------------------------------------------------
467 // Pre Step 6: pad if not the same size : it can occur when downsample and upsample
468 //if (!HaveSameSizeAndSpacing(working_image, input)) {
469 if (working_image->GetLargestPossibleRegion() != input->GetLargestPossibleRegion()) {
470 StartNewStep("Pad to get the same size than input");
471 typename ImageType::Pointer temp = ImageType::New();
472 temp->CopyInformation(input);
473 temp->SetRegions(input->GetLargestPossibleRegion()); // Do not forget !!
475 temp->FillBuffer(m_BackgroundValue);
476 typename PasteFilterType::Pointer padFilter2 = PasteFilterType::New();
477 padFilter2->SetSourceImage(working_image);
478 padFilter2->SetDestinationImage(temp);
479 padFilter2->SetDestinationIndex(input->GetLargestPossibleRegion().GetIndex());
480 padFilter2->SetSourceRegion(working_image->GetLargestPossibleRegion());
481 padFilter2->Update();
482 working_image = padFilter2->GetOutput();
483 StopCurrentStep<ImageType>(working_image);
486 //DD("[debug] Rel Pos : no padding after");
489 //--------------------------------------------------------------------
490 //--------------------------------------------------------------------
491 // Step 6: combine input+thresholded relpos
492 StartNewStep("Combine with initial input (boolean And)");
493 typedef clitk::BooleanOperatorLabelImageFilter<ImageType> BoolFilterType;
494 typename BoolFilterType::Pointer combineFilter = BoolFilterType::New();
495 combineFilter->SetBackgroundValue(m_BackgroundValue);
496 combineFilter->SetBackgroundValue1(m_BackgroundValue);
497 combineFilter->SetBackgroundValue2(m_BackgroundValue);
498 combineFilter->SetForegroundValue(m_BackgroundValue+1);
499 combineFilter->SetInput1(input);
500 combineFilter->SetInput2(working_image);
501 if (GetInverseOrientationFlag())
502 combineFilter->SetOperationType(BoolFilterType::AndNot);
504 if (GetCombineWithOrFlag())
505 combineFilter->SetOperationType(BoolFilterType::Or);
507 combineFilter->SetOperationType(BoolFilterType::And);
509 combineFilter->InPlaceOff(); // Do not modify initial input (!)
510 combineFilter->Update();
511 working_image = combineFilter->GetOutput();
513 // Remove (if needed the object from the support)
514 if (GetRemoveObjectFlag()) {
515 combineFilter = BoolFilterType::New();
516 combineFilter->SetInput1(working_image);
517 combineFilter->SetInput2(object);
518 combineFilter->SetOperationType(BoolFilterType::AndNot);
519 combineFilter->InPlaceOn();
520 combineFilter->Update();
521 working_image = combineFilter->GetOutput();
523 // In the other case, we must *add* the initial object to keep it
524 // but not more than the initial support
526 combineFilter = BoolFilterType::New();
527 combineFilter->SetInput1(working_image);
528 combineFilter->SetInput2(object);
529 combineFilter->SetOperationType(BoolFilterType::Or);
530 combineFilter->InPlaceOn();
531 combineFilter->Update();
532 working_image = combineFilter->GetOutput(); // not needed because InPlaceOn ?
533 combineFilter = BoolFilterType::New();
534 combineFilter->SetInput1(working_image);
535 combineFilter->SetInput2(input);
536 combineFilter->SetOperationType(BoolFilterType::And);
537 combineFilter->InPlaceOn();
538 combineFilter->Update();
539 working_image = combineFilter->GetOutput();
542 StopCurrentStep<ImageType>(working_image);
544 //--------------------------------------------------------------------
545 //--------------------------------------------------------------------
547 if (GetAutoCropFlag()) {
548 StartNewStep("Final AutoCrop");
549 typedef clitk::AutoCropFilter<ImageType> CropFilterType;
550 typename CropFilterType::Pointer cropFilter = CropFilterType::New();
551 cropFilter->SetInput(working_image);
552 cropFilter->ReleaseDataFlagOff();
553 cropFilter->Update();
554 working_image = cropFilter->GetOutput();
555 StopCurrentStep<ImageType>(working_image);
558 //--------------------------------------------------------------------
559 //--------------------------------------------------------------------
561 // Final Step -> set output
562 this->SetNthOutput(0, working_image);
563 // this->GraftOutput(working_image);
565 //--------------------------------------------------------------------