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 "clitkSetBackgroundImageFilter.h"
21 #include "clitkSliceBySliceRelativePositionFilter.h"
22 #include "clitkCropLikeImageFilter.h"
23 #include "clitkMemoryUsage.h"
26 #include <itkConnectedComponentImageFilter.h>
27 #include <itkRelabelComponentImageFilter.h>
28 #include <itkBinaryThresholdImageFilter.h>
29 #include <itkPasteImageFilter.h>
30 #include <itkStatisticsLabelMapFilter.h>
31 #include <itkBinaryBallStructuringElement.h>
32 #include <itkBinaryDilateImageFilter.h>
33 #include <itkConstantPadImageFilter.h>
34 #include <itkImageSliceIteratorWithIndex.h>
35 #include <itkBinaryMorphologicalOpeningImageFilter.h>
36 #include <itkImageDuplicator.h>
37 #include <itkSignedMaurerDistanceMapImageFilter.h>
41 //--------------------------------------------------------------------
42 template<class ImageType, class TMaskImageType>
43 typename ImageType::Pointer
44 SetBackground(const ImageType * input,
45 const TMaskImageType * mask,
46 typename TMaskImageType::PixelType maskBG,
47 typename ImageType::PixelType outValue,
49 typedef SetBackgroundImageFilter<ImageType, TMaskImageType, ImageType>
50 SetBackgroundImageFilterType;
51 typename SetBackgroundImageFilterType::Pointer setBackgroundFilter
52 = SetBackgroundImageFilterType::New();
53 // if (inPlace) setBackgroundFilter->ReleaseDataFlagOn(); // No seg fault
54 setBackgroundFilter->SetInPlace(inPlace); // This is important to keep memory low
55 setBackgroundFilter->SetInput(input);
56 setBackgroundFilter->SetInput2(mask);
57 setBackgroundFilter->SetMaskValue(maskBG);
58 setBackgroundFilter->SetOutsideValue(outValue);
59 setBackgroundFilter->Update();
60 return setBackgroundFilter->GetOutput();
62 //--------------------------------------------------------------------
65 //--------------------------------------------------------------------
66 template<class ImageType>
67 int GetNumberOfConnectedComponentLabels(const ImageType * input,
68 typename ImageType::PixelType BG,
69 bool isFullyConnected) {
70 // Connected Component label
71 typedef itk::ConnectedComponentImageFilter<ImageType, ImageType> ConnectFilterType;
72 typename ConnectFilterType::Pointer connectFilter = ConnectFilterType::New();
73 connectFilter->SetInput(input);
74 connectFilter->SetBackgroundValue(BG);
75 connectFilter->SetFullyConnected(isFullyConnected);
76 connectFilter->Update();
79 return connectFilter->GetObjectCount();
81 //--------------------------------------------------------------------
83 //--------------------------------------------------------------------
85 Warning : in this cas, we consider outputType like inputType, not
86 InternalImageType. Be sure it fits.
88 template<class ImageType>
89 typename ImageType::Pointer
90 Labelize(const ImageType * input,
91 typename ImageType::PixelType BG,
92 bool isFullyConnected,
93 int minimalComponentSize) {
94 // InternalImageType for storing large number of component
95 typedef itk::Image<int, ImageType::ImageDimension> InternalImageType;
97 // Connected Component label
98 typedef itk::ConnectedComponentImageFilter<ImageType, InternalImageType> ConnectFilterType;
99 typename ConnectFilterType::Pointer connectFilter = ConnectFilterType::New();
100 // connectFilter->ReleaseDataFlagOn();
101 connectFilter->SetInput(input);
102 connectFilter->SetBackgroundValue(BG);
103 connectFilter->SetFullyConnected(isFullyConnected);
105 // Sort by size and remove too small area.
106 typedef itk::RelabelComponentImageFilter<InternalImageType, ImageType> RelabelFilterType;
107 typename RelabelFilterType::Pointer relabelFilter = RelabelFilterType::New();
108 // relabelFilter->ReleaseDataFlagOn(); // if yes, fail when ExplosionControlledThresholdConnectedImageFilter ???
109 relabelFilter->SetInput(connectFilter->GetOutput());
110 relabelFilter->SetMinimumObjectSize(minimalComponentSize);
111 relabelFilter->Update();
114 typename ImageType::Pointer output = relabelFilter->GetOutput();
117 //--------------------------------------------------------------------
120 //--------------------------------------------------------------------
122 Warning : in this cas, we consider outputType like inputType, not
123 InternalImageType. Be sure it fits.
125 template<class ImageType>
126 typename ImageType::Pointer
127 LabelizeAndCountNumberOfObjects(const ImageType * input,
128 typename ImageType::PixelType BG,
129 bool isFullyConnected,
130 int minimalComponentSize,
132 // InternalImageType for storing large number of component
133 typedef itk::Image<int, ImageType::ImageDimension> InternalImageType;
135 // Connected Component label
136 typedef itk::ConnectedComponentImageFilter<ImageType, InternalImageType> ConnectFilterType;
137 typename ConnectFilterType::Pointer connectFilter = ConnectFilterType::New();
138 // connectFilter->ReleaseDataFlagOn();
139 connectFilter->SetInput(input);
140 connectFilter->SetBackgroundValue(BG);
141 connectFilter->SetFullyConnected(isFullyConnected);
143 // Sort by size and remove too small area.
144 typedef itk::RelabelComponentImageFilter<InternalImageType, ImageType> RelabelFilterType;
145 typename RelabelFilterType::Pointer relabelFilter = RelabelFilterType::New();
146 // relabelFilter->ReleaseDataFlagOn(); // if yes, fail when ExplosionControlledThresholdConnectedImageFilter ???
147 relabelFilter->SetInput(connectFilter->GetOutput());
148 relabelFilter->SetMinimumObjectSize(minimalComponentSize);
149 relabelFilter->Update();
151 nb = relabelFilter->GetNumberOfObjects();
152 // DD(relabelFilter->GetOriginalNumberOfObjects());
153 // DD(relabelFilter->GetSizeOfObjectsInPhysicalUnits()[0]);
156 typename ImageType::Pointer output = relabelFilter->GetOutput();
159 //--------------------------------------------------------------------
163 //--------------------------------------------------------------------
164 template<class ImageType>
165 typename ImageType::Pointer
166 RemoveLabels(const ImageType * input,
167 typename ImageType::PixelType BG,
168 std::vector<typename ImageType::PixelType> & labelsToRemove) {
169 assert(labelsToRemove.size() != 0);
170 typename ImageType::Pointer working_image;// = input;
171 for (unsigned int i=0; i <labelsToRemove.size(); i++) {
172 typedef SetBackgroundImageFilter<ImageType, ImageType> SetBackgroundImageFilterType;
173 typename SetBackgroundImageFilterType::Pointer setBackgroundFilter = SetBackgroundImageFilterType::New();
174 setBackgroundFilter->SetInput(input);
175 setBackgroundFilter->SetInput2(input);
176 setBackgroundFilter->SetMaskValue(labelsToRemove[i]);
177 setBackgroundFilter->SetOutsideValue(BG);
178 setBackgroundFilter->Update();
179 working_image = setBackgroundFilter->GetOutput();
181 return working_image;
183 //--------------------------------------------------------------------
186 //--------------------------------------------------------------------
187 template<class ImageType>
188 typename ImageType::Pointer
189 KeepLabels(const ImageType * input,
190 typename ImageType::PixelType BG,
191 typename ImageType::PixelType FG,
192 typename ImageType::PixelType firstKeep,
193 typename ImageType::PixelType lastKeep,
195 typedef itk::BinaryThresholdImageFilter<ImageType, ImageType> BinarizeFilterType;
196 typename BinarizeFilterType::Pointer binarizeFilter = BinarizeFilterType::New();
197 binarizeFilter->SetInput(input);
198 binarizeFilter->SetLowerThreshold(firstKeep);
199 if (useLastKeep) binarizeFilter->SetUpperThreshold(lastKeep);
200 binarizeFilter->SetInsideValue(FG);
201 binarizeFilter->SetOutsideValue(BG);
202 binarizeFilter->Update();
203 return binarizeFilter->GetOutput();
205 //--------------------------------------------------------------------
208 //--------------------------------------------------------------------
209 template<class ImageType>
210 typename ImageType::Pointer
211 LabelizeAndSelectLabels(const ImageType * input,
212 typename ImageType::PixelType BG,
213 typename ImageType::PixelType FG,
214 bool isFullyConnected,
215 int minimalComponentSize,
216 LabelizeParameters<typename ImageType::PixelType> * param)
218 typename ImageType::Pointer working_image;
219 working_image = Labelize<ImageType>(input, BG, isFullyConnected, minimalComponentSize);
220 if (param->GetLabelsToRemove().size() != 0)
221 working_image = RemoveLabels<ImageType>(working_image, BG, param->GetLabelsToRemove());
222 working_image = KeepLabels<ImageType>(working_image,
224 param->GetFirstKeep(),
225 param->GetLastKeep(),
226 param->GetUseLastKeep());
227 return working_image;
229 //--------------------------------------------------------------------
232 //--------------------------------------------------------------------
233 template<class MaskImageType>
234 typename MaskImageType::Pointer
235 SliceBySliceRelativePosition(const MaskImageType * input,
236 const MaskImageType * object,
239 std::string orientation,
240 bool uniqueConnectedComponent,
243 bool singleObjectCCL)
245 typedef clitk::SliceBySliceRelativePositionFilter<MaskImageType> SliceRelPosFilterType;
246 typename SliceRelPosFilterType::Pointer sliceRelPosFilter = SliceRelPosFilterType::New();
247 sliceRelPosFilter->VerboseStepFlagOff();
248 sliceRelPosFilter->WriteStepFlagOff();
249 sliceRelPosFilter->SetInput(input);
250 sliceRelPosFilter->SetInputObject(object);
251 sliceRelPosFilter->SetDirection(direction);
252 sliceRelPosFilter->SetFuzzyThreshold(threshold);
253 sliceRelPosFilter->AddOrientationTypeString(orientation);
254 sliceRelPosFilter->SetIntermediateSpacingFlag((spacing != -1));
255 sliceRelPosFilter->SetIntermediateSpacing(spacing);
256 sliceRelPosFilter->SetUniqueConnectedComponentBySliceFlag(uniqueConnectedComponent);
257 sliceRelPosFilter->ObjectCCLSelectionFlagOff();
258 sliceRelPosFilter->SetUseTheLargestObjectCCLFlag(singleObjectCCL);
259 // sliceRelPosFilter->SetInverseOrientationFlag(inverseflag);
260 sliceRelPosFilter->SetAutoCropFlag(autocropFlag);
261 sliceRelPosFilter->IgnoreEmptySliceObjectFlagOn();
262 sliceRelPosFilter->Update();
263 return sliceRelPosFilter->GetOutput();
265 //--------------------------------------------------------------------
268 //--------------------------------------------------------------------
269 template<class MaskImageType>
270 typename MaskImageType::Pointer
271 SliceBySliceRelativePosition(const MaskImageType * input,
272 const MaskImageType * object,
277 bool uniqueConnectedComponent,
280 bool singleObjectCCL)
282 typedef clitk::SliceBySliceRelativePositionFilter<MaskImageType> SliceRelPosFilterType;
283 typename SliceRelPosFilterType::Pointer sliceRelPosFilter = SliceRelPosFilterType::New();
284 sliceRelPosFilter->VerboseStepFlagOff();
285 sliceRelPosFilter->WriteStepFlagOff();
286 sliceRelPosFilter->SetInput(input);
287 sliceRelPosFilter->SetInputObject(object);
288 sliceRelPosFilter->SetDirection(direction);
289 sliceRelPosFilter->SetFuzzyThreshold(threshold);
290 // sliceRelPosFilter->AddOrientationTypeString(orientation);
291 sliceRelPosFilter->AddAnglesInRad(angle, 0.0);
292 sliceRelPosFilter->SetIntermediateSpacingFlag((spacing != -1));
293 sliceRelPosFilter->SetIntermediateSpacing(spacing);
294 sliceRelPosFilter->SetUniqueConnectedComponentBySliceFlag(uniqueConnectedComponent);
295 sliceRelPosFilter->ObjectCCLSelectionFlagOff();
296 sliceRelPosFilter->SetUseTheLargestObjectCCLFlag(singleObjectCCL);
297 sliceRelPosFilter->SetInverseOrientationFlag(inverseflag);
298 sliceRelPosFilter->SetAutoCropFlag(autocropFlag);
299 sliceRelPosFilter->IgnoreEmptySliceObjectFlagOn();
300 sliceRelPosFilter->Update();
301 return sliceRelPosFilter->GetOutput();
303 //--------------------------------------------------------------------
306 //--------------------------------------------------------------------
307 template<class ImageType>
309 FindExtremaPointInAGivenDirection(const ImageType * input,
310 typename ImageType::PixelType bg,
311 int direction, bool opposite,
312 typename ImageType::PointType & point)
314 typename ImageType::PointType dummy;
315 return FindExtremaPointInAGivenDirection(input, bg, direction,
316 opposite, dummy, 0, point);
318 //--------------------------------------------------------------------
321 //--------------------------------------------------------------------
322 template<class ImageType>
324 FindExtremaPointInAGivenDirection(const ImageType * input,
325 typename ImageType::PixelType bg,
326 int direction, bool opposite,
327 typename ImageType::PointType refpoint,
329 typename ImageType::PointType & point)
332 loop over input pixels, store the index in the fg that is max
333 according to the given direction.
335 typedef itk::ImageRegionConstIteratorWithIndex<ImageType> IteratorType;
336 IteratorType iter(input, input->GetLargestPossibleRegion());
338 typename ImageType::IndexType max = input->GetLargestPossibleRegion().GetIndex();
339 if (opposite) max = max+input->GetLargestPossibleRegion().GetSize();
341 while (!iter.IsAtEnd()) {
342 if (iter.Get() != bg) {
343 bool test = iter.GetIndex()[direction] > max[direction];
344 if (opposite) test = !test;
346 typename ImageType::PointType p;
347 input->TransformIndexToPhysicalPoint(iter.GetIndex(), p);
348 if ((distanceMax==0) || (p.EuclideanDistanceTo(refpoint) < distanceMax)) {
349 max = iter.GetIndex();
356 if (!found) return false;
357 input->TransformIndexToPhysicalPoint(max, point); // half of the pixel
360 //--------------------------------------------------------------------
363 //--------------------------------------------------------------------
364 template<class ImageType>
365 typename ImageType::Pointer
366 CropImageRemoveGreaterThan(const ImageType * image,
367 int dim, double min, bool autoCrop,
368 typename ImageType::PixelType BG)
370 return CropImageAlongOneAxis<ImageType>(image, dim,
371 image->GetOrigin()[dim],
375 //--------------------------------------------------------------------
378 //--------------------------------------------------------------------
379 template<class ImageType>
380 typename ImageType::Pointer
381 CropImageRemoveLowerThan(const ImageType * image,
382 int dim, double max, bool autoCrop,
383 typename ImageType::PixelType BG)
385 typename ImageType::PointType p;
387 image->TransformIndexToPhysicalPoint(image->GetLargestPossibleRegion().GetIndex()+
388 image->GetLargestPossibleRegion().GetSize(), p);
390 return CropImageAlongOneAxis<ImageType>(image, dim, max, p[dim], autoCrop, BG);
392 //--------------------------------------------------------------------
395 //--------------------------------------------------------------------
396 template<class ImageType>
397 typename ImageType::Pointer
398 CropImageAlongOneAxis(const ImageType * image,
399 int dim, double min, double max,
400 bool autoCrop, typename ImageType::PixelType BG)
402 // Compute region size
403 typename ImageType::RegionType region;
404 typename ImageType::SizeType size = image->GetLargestPossibleRegion().GetSize();
407 typename ImageType::PointType p = image->GetOrigin(); // not at pixel center !
408 if (min > p[dim]) p[dim] = min; // Check if not outside the image
409 typename ImageType::IndexType start;
410 image->TransformPhysicalPointToIndex(p, start);
412 // Size of the region
413 // -1 because last point is size -1
414 double m = image->GetOrigin()[dim] + (size[dim]-1)*image->GetSpacing()[dim];
415 if (max > m) p[dim] = m; // Check if not outside the image
418 typename ImageType::IndexType end;
419 image->TransformPhysicalPointToIndex(p, end);
420 size[dim] = abs(end[dim]-start[dim])+1;// +1 because we want to include the point.
423 region.SetIndex(start);
424 region.SetSize(size);
427 typedef itk::RegionOfInterestImageFilter<ImageType, ImageType> CropFilterType;
428 typename CropFilterType::Pointer cropFilter = CropFilterType::New();
429 cropFilter->SetInput(image);
430 cropFilter->SetRegionOfInterest(region);
431 cropFilter->Update();
432 typename ImageType::Pointer result = cropFilter->GetOutput();
436 result = AutoCrop<ImageType>(result, BG);
440 //--------------------------------------------------------------------
443 //--------------------------------------------------------------------
444 template<class ImageType>
446 ComputeCentroids(const ImageType * image,
447 typename ImageType::PixelType BG,
448 std::vector<typename ImageType::PointType> & centroids)
450 typedef long LabelType;
451 static const unsigned int Dim = ImageType::ImageDimension;
452 typedef itk::ShapeLabelObject< LabelType, Dim > LabelObjectType;
453 typedef itk::LabelMap< LabelObjectType > LabelMapType;
454 typedef itk::LabelImageToLabelMapFilter<ImageType, LabelMapType> ImageToMapFilterType;
455 typename ImageToMapFilterType::Pointer imageToLabelFilter = ImageToMapFilterType::New();
456 typedef itk::ShapeLabelMapFilter<LabelMapType, ImageType> ShapeFilterType;
457 typename ShapeFilterType::Pointer statFilter = ShapeFilterType::New();
458 imageToLabelFilter->SetBackgroundValue(BG);
459 imageToLabelFilter->SetInput(image);
460 statFilter->SetInput(imageToLabelFilter->GetOutput());
461 statFilter->Update();
462 typename LabelMapType::Pointer labelMap = statFilter->GetOutput();
465 typename ImageType::PointType dummy;
466 centroids.push_back(dummy); // label 0 -> no centroid, use dummy point for BG
467 //DS FIXME (not useful ! to change ..)
468 for(uint i=0; i<labelMap->GetNumberOfLabelObjects(); i++) {
469 int label = labelMap->GetLabels()[i];
470 centroids.push_back(labelMap->GetLabelObject(label)->GetCentroid());
473 //--------------------------------------------------------------------
476 //--------------------------------------------------------------------
477 template<class ImageType, class LabelType>
478 typename itk::LabelMap< itk::ShapeLabelObject<LabelType, ImageType::ImageDimension> >::Pointer
479 ComputeLabelMap(const ImageType * image,
480 typename ImageType::PixelType BG,
481 bool computePerimeterFlag)
483 static const unsigned int Dim = ImageType::ImageDimension;
484 typedef itk::ShapeLabelObject< LabelType, Dim > LabelObjectType;
485 typedef itk::LabelMap< LabelObjectType > LabelMapType;
486 typedef itk::LabelImageToLabelMapFilter<ImageType, LabelMapType> ImageToMapFilterType;
487 typename ImageToMapFilterType::Pointer imageToLabelFilter = ImageToMapFilterType::New();
488 typedef itk::ShapeLabelMapFilter<LabelMapType, ImageType> ShapeFilterType;
489 typename ShapeFilterType::Pointer statFilter = ShapeFilterType::New();
490 imageToLabelFilter->SetBackgroundValue(BG);
491 imageToLabelFilter->SetInput(image);
492 statFilter->SetInput(imageToLabelFilter->GetOutput());
493 statFilter->SetComputePerimeter(computePerimeterFlag);
494 statFilter->Update();
495 return statFilter->GetOutput();
497 //--------------------------------------------------------------------
500 //--------------------------------------------------------------------
501 template<class ImageType>
503 ComputeCentroids2(const ImageType * image,
504 typename ImageType::PixelType BG,
505 std::vector<typename ImageType::PointType> & centroids)
507 typedef long LabelType;
508 static const unsigned int Dim = ImageType::ImageDimension;
509 typedef itk::ShapeLabelObject< LabelType, Dim > LabelObjectType;
510 typedef itk::LabelMap< LabelObjectType > LabelMapType;
511 typedef itk::LabelImageToLabelMapFilter<ImageType, LabelMapType> ImageToMapFilterType;
512 typename ImageToMapFilterType::Pointer imageToLabelFilter = ImageToMapFilterType::New();
513 typedef itk::ShapeLabelMapFilter<LabelMapType, ImageType> ShapeFilterType;
514 typename ShapeFilterType::Pointer statFilter = ShapeFilterType::New();
515 imageToLabelFilter->SetBackgroundValue(BG);
516 imageToLabelFilter->SetInput(image);
517 statFilter->SetInput(imageToLabelFilter->GetOutput());
518 statFilter->Update();
519 typename LabelMapType::Pointer labelMap = statFilter->GetOutput();
522 typename ImageType::PointType dummy;
523 centroids.push_back(dummy); // label 0 -> no centroid, use dummy point
524 for(uint i=1; i<labelMap->GetNumberOfLabelObjects()+1; i++) {
525 centroids.push_back(labelMap->GetLabelObject(i)->GetCentroid());
528 for(uint i=1; i<labelMap->GetNumberOfLabelObjects()+1; i++) {
529 DD(labelMap->GetLabelObject(i)->GetBinaryPrincipalAxes());
530 DD(labelMap->GetLabelObject(i)->GetBinaryFlatness());
531 DD(labelMap->GetLabelObject(i)->GetRoundness ());
533 // search for the point on the boundary alog PA
538 //--------------------------------------------------------------------
541 //--------------------------------------------------------------------
542 template<class ImageType>
544 PointsUtils<ImageType>::Convert2DTo3D(const PointType2D & p2D,
545 const ImageType * image,
550 index3D[0] = index3D[1] = 0;
551 index3D[2] = image->GetLargestPossibleRegion().GetIndex()[2]+slice;
552 image->TransformIndexToPhysicalPoint(index3D, p3D);
555 // p3D[2] = p[2];//(image->GetLargestPossibleRegion().GetIndex()[2]+slice)*image->GetSpacing()[2]
556 // + image->GetOrigin()[2];
558 //--------------------------------------------------------------------
561 //--------------------------------------------------------------------
562 template<class ImageType>
564 PointsUtils<ImageType>::Convert2DMapTo3DList(const MapPoint2DType & map,
565 const ImageType * image,
566 VectorPoint3DType & list)
568 typename MapPoint2DType::const_iterator iter = map.begin();
569 while (iter != map.end()) {
571 Convert2DTo3D(iter->second, image, iter->first, p);
576 //--------------------------------------------------------------------
579 //--------------------------------------------------------------------
580 template<class ImageType>
582 PointsUtils<ImageType>::Convert2DListTo3DList(const VectorPoint2DType & p2D,
584 const ImageType * image,
585 VectorPoint3DType & list)
587 for(uint i=0; i<p2D.size(); i++) {
589 Convert2DTo3D(p2D[i], image, slice, p);
593 //--------------------------------------------------------------------
596 //--------------------------------------------------------------------
597 template<class ImageType>
599 WriteListOfLandmarks(std::vector<typename ImageType::PointType> points,
600 std::string filename)
603 openFileForWriting(os, filename);
604 os << "LANDMARKS1" << std::endl;
605 for(uint i=0; i<points.size(); i++) {
606 const typename ImageType::PointType & p = points[i];
607 // Write it in the file
608 os << i << " " << p[0] << " " << p[1] << " " << p[2] << " 0 0 " << std::endl;
612 //--------------------------------------------------------------------
615 //--------------------------------------------------------------------
616 template<class ImageType>
617 typename ImageType::Pointer
618 Dilate(const ImageType * image, double radiusInMM,
619 typename ImageType::PixelType BG,
620 typename ImageType::PixelType FG,
623 typename ImageType::SizeType r;
624 for(uint i=0; i<ImageType::ImageDimension; i++)
625 r[i] = (uint)lrint(radiusInMM/image->GetSpacing()[i]);
626 return Dilate<ImageType>(image, r, BG, FG, extendSupport);
628 //--------------------------------------------------------------------
631 //--------------------------------------------------------------------
632 template<class ImageType>
633 typename ImageType::Pointer
634 Dilate(const ImageType * image, typename ImageType::PointType radiusInMM,
635 typename ImageType::PixelType BG,
636 typename ImageType::PixelType FG,
639 typename ImageType::SizeType r;
640 for(uint i=0; i<ImageType::ImageDimension; i++)
641 r[i] = (uint)lrint(radiusInMM[i]/image->GetSpacing()[i]);
642 return Dilate<ImageType>(image, r, BG, FG, extendSupport);
644 //--------------------------------------------------------------------
647 //--------------------------------------------------------------------
648 template<class ImageType>
649 typename ImageType::Pointer
650 Dilate(const ImageType * image, typename ImageType::SizeType radius,
651 typename ImageType::PixelType BG,
652 typename ImageType::PixelType FG,
655 // Create kernel for dilatation
656 typedef itk::BinaryBallStructuringElement<typename ImageType::PixelType,
657 ImageType::ImageDimension> KernelType;
658 KernelType structuringElement;
659 structuringElement.SetRadius(radius);
660 structuringElement.CreateStructuringElement();
662 typename ImageType::Pointer output;
664 typedef itk::ConstantPadImageFilter<ImageType, ImageType> PadFilterType;
665 typename PadFilterType::Pointer padFilter = PadFilterType::New();
666 padFilter->SetInput(image);
667 typename ImageType::SizeType lower;
668 typename ImageType::SizeType upper;
669 for(uint i=0; i<3; i++) {
670 lower[i] = upper[i] = 2*(radius[i]+1);
672 padFilter->SetPadLowerBound(lower);
673 padFilter->SetPadUpperBound(upper);
675 output = padFilter->GetOutput();
679 typedef itk::BinaryDilateImageFilter<ImageType, ImageType , KernelType> DilateFilterType;
680 typename DilateFilterType::Pointer dilateFilter = DilateFilterType::New();
681 dilateFilter->SetBackgroundValue(BG);
682 dilateFilter->SetForegroundValue(FG);
683 dilateFilter->SetBoundaryToForeground(false);
684 dilateFilter->SetKernel(structuringElement);
685 if (extendSupport) dilateFilter->SetInput(output);
686 else dilateFilter->SetInput(image);
687 dilateFilter->Update();
688 return dilateFilter->GetOutput();
690 //--------------------------------------------------------------------
693 //--------------------------------------------------------------------
694 template<class ImageType>
695 typename ImageType::Pointer
696 Opening(const ImageType * image, typename ImageType::SizeType radius,
697 typename ImageType::PixelType BG,
698 typename ImageType::PixelType FG)
701 typedef itk::BinaryBallStructuringElement<typename ImageType::PixelType,
702 ImageType::ImageDimension> KernelType;
703 KernelType structuringElement;
704 structuringElement.SetRadius(radius);
705 structuringElement.CreateStructuringElement();
708 typedef itk::BinaryMorphologicalOpeningImageFilter<ImageType, ImageType , KernelType> OpeningFilterType;
709 typename OpeningFilterType::Pointer open = OpeningFilterType::New();
710 open->SetInput(image);
711 open->SetBackgroundValue(BG);
712 open->SetForegroundValue(FG);
713 open->SetKernel(structuringElement);
715 return open->GetOutput();
717 //--------------------------------------------------------------------
721 //--------------------------------------------------------------------
722 template<class ValueType, class VectorType>
723 void ConvertOption(std::string optionName, uint given,
724 ValueType * values, VectorType & p,
725 uint dim, bool required)
727 if (required && (given == 0)) {
728 clitkExceptionMacro("The option --" << optionName << " must be set and have 1 or "
729 << dim << " values.");
732 for(uint i=0; i<dim; i++) p[i] = values[0];
736 for(uint i=0; i<dim; i++) p[i] = values[i];
739 if (given == 0) return;
740 clitkExceptionMacro("The option --" << optionName << " must have 1 or "
741 << dim << " values.");
743 //--------------------------------------------------------------------
746 //--------------------------------------------------------------------
748 http://www.gamedev.net/community/forums/topic.asp?topic_id=542870
749 Assuming the points are (Ax,Ay) (Bx,By) and (Cx,Cy), you need to compute:
750 (Bx - Ax) * (Cy - Ay) - (By - Ay) * (Cx - Ax)
751 This will equal zero if the point C is on the line formed by
752 points A and B, and will have a different sign depending on the
753 side. Which side this is depends on the orientation of your (x,y)
754 coordinates, but you can plug test values for A,B and C into this
755 formula to determine whether negative values are to the left or to
757 => to accelerate, start with formula, when change sign -> stop and fill
759 offsetToKeep = is used to determine which side of the line we
760 keep. The point along the mainDirection but 'offsetToKeep' mm away
764 template<class ImageType>
766 SliceBySliceSetBackgroundFromLineSeparation(ImageType * input,
767 std::vector<typename ImageType::PointType> & lA,
768 std::vector<typename ImageType::PointType> & lB,
769 typename ImageType::PixelType BG,
774 assert((mainDirection==0) || (mainDirection==1));
775 typename ImageType::PointType offset;
776 offset[0] = offset[1] = offset[2] = 0.0;
777 offset[mainDirection] = offsetToKeep;
778 SliceBySliceSetBackgroundFromLineSeparation<ImageType>(input, lA, lB, BG, offset, keepIfEqual);
780 template<class ImageType>
782 SliceBySliceSetBackgroundFromLineSeparation(ImageType * input,
783 std::vector<typename ImageType::PointType> & lA,
784 std::vector<typename ImageType::PointType> & lB,
785 typename ImageType::PixelType BG,
786 typename ImageType::PointType offsetToKeep,
789 typedef itk::ImageSliceIteratorWithIndex<ImageType> SliceIteratorType;
790 SliceIteratorType siter = SliceIteratorType(input, input->GetLargestPossibleRegion());
791 siter.SetFirstDirection(0);
792 siter.SetSecondDirection(1);
795 typename ImageType::PointType A;
796 typename ImageType::PointType B;
797 typename ImageType::PointType C;
798 assert(lA.size() == lB.size());
799 while ((i<lA.size()) && (!siter.IsAtEnd())) {
800 // Check that the current slice correspond to the current point
801 input->TransformIndexToPhysicalPoint(siter.GetIndex(), C);
802 if ((fabs(C[2] - lA[i][2]))>0.01) { // is !equal with a tolerance of 0.01 mm
805 // Define A,B,C points
809 // Check that the line is not a point (A=B)
810 bool p = (A[0] == B[0]) && (A[1] == B[1]);
813 //C[mainDirection] += offsetToKeep; // I know I must keep this point
814 C[0] += offsetToKeep[0];
815 C[1] += offsetToKeep[1];
816 //C[2] += offsetToKeep[2];
817 double s = (B[0] - A[0]) * (C[1] - A[1]) - (B[1] - A[1]) * (C[0] - A[0]);
818 bool isPositive = s<0;
819 while (!siter.IsAtEndOfSlice()) {
820 while (!siter.IsAtEndOfLine()) {
821 // Very slow, I know ... but image should be very small
822 input->TransformIndexToPhysicalPoint(siter.GetIndex(), C);
823 double s = (B[0] - A[0]) * (C[1] - A[1]) - (B[1] - A[1]) * (C[0] - A[0]);
825 if (!keepIfEqual) siter.Set(BG); // on the line, we decide to remove
828 if (s > 0) siter.Set(BG);
831 if (s < 0) siter.Set(BG);
840 } // End of current slice
844 //--------------------------------------------------------------------
847 //--------------------------------------------------------------------
848 template<class ImageType>
850 AndNot(ImageType * input,
851 const ImageType * object,
852 typename ImageType::PixelType BG)
854 typename ImageType::Pointer o;
856 if (!clitk::HaveSameSizeAndSpacing<ImageType, ImageType>(input, object)) {
857 o = clitk::ResizeImageLike<ImageType>(object, input, BG);
861 typedef clitk::BooleanOperatorLabelImageFilter<ImageType> BoolFilterType;
862 typename BoolFilterType::Pointer boolFilter = BoolFilterType::New();
863 boolFilter->InPlaceOn();
864 boolFilter->SetInput1(input);
865 if (resized) boolFilter->SetInput2(o);
866 else boolFilter->SetInput2(object);
867 boolFilter->SetBackgroundValue1(BG);
868 boolFilter->SetBackgroundValue2(BG);
869 boolFilter->SetOperationType(BoolFilterType::AndNot);
870 boolFilter->Update();
872 //--------------------------------------------------------------------
875 //--------------------------------------------------------------------
876 template<class ImageType>
878 And(ImageType * input,
879 const ImageType * object,
880 typename ImageType::PixelType BG)
882 typename ImageType::Pointer o;
884 if (!clitk::HaveSameSizeAndSpacing<ImageType, ImageType>(input, object)) {
885 o = clitk::ResizeImageLike<ImageType>(object, input, BG);
889 typedef clitk::BooleanOperatorLabelImageFilter<ImageType> BoolFilterType;
890 typename BoolFilterType::Pointer boolFilter = BoolFilterType::New();
891 boolFilter->InPlaceOn();
892 boolFilter->SetInput1(input);
893 if (resized) boolFilter->SetInput2(o);
894 else boolFilter->SetInput2(object);
895 boolFilter->SetBackgroundValue1(BG);
896 boolFilter->SetBackgroundValue2(BG);
897 boolFilter->SetOperationType(BoolFilterType::And);
898 boolFilter->Update();
900 //--------------------------------------------------------------------
903 //--------------------------------------------------------------------
904 template<class ImageType>
906 Or(ImageType * input,
907 const ImageType * object,
908 typename ImageType::PixelType BG)
910 typename ImageType::Pointer o;
912 if (!clitk::HaveSameSizeAndSpacing<ImageType, ImageType>(input, object)) {
913 o = clitk::ResizeImageLike<ImageType>(object, input, BG);
917 typedef clitk::BooleanOperatorLabelImageFilter<ImageType> BoolFilterType;
918 typename BoolFilterType::Pointer boolFilter = BoolFilterType::New();
919 boolFilter->InPlaceOn();
920 boolFilter->SetInput1(input);
921 if (resized) boolFilter->SetInput2(o);
922 else boolFilter->SetInput2(object);
923 boolFilter->SetBackgroundValue1(BG);
924 boolFilter->SetBackgroundValue2(BG);
925 boolFilter->SetOperationType(BoolFilterType::Or);
926 boolFilter->Update();
928 //--------------------------------------------------------------------
931 //--------------------------------------------------------------------
932 template<class ImageType>
933 typename ImageType::Pointer
934 Binarize(const ImageType * input,
935 typename ImageType::PixelType lower,
936 typename ImageType::PixelType upper,
937 typename ImageType::PixelType BG,
938 typename ImageType::PixelType FG)
940 typedef itk::BinaryThresholdImageFilter<ImageType, ImageType> BinaryThresholdFilterType;
941 typename BinaryThresholdFilterType::Pointer binarizeFilter = BinaryThresholdFilterType::New();
942 binarizeFilter->SetInput(input);
943 binarizeFilter->InPlaceOff();
944 binarizeFilter->SetLowerThreshold(lower);
945 binarizeFilter->SetUpperThreshold(upper);
946 binarizeFilter->SetInsideValue(FG);
947 binarizeFilter->SetOutsideValue(BG);
948 binarizeFilter->Update();
949 return binarizeFilter->GetOutput();
951 //--------------------------------------------------------------------
954 //--------------------------------------------------------------------
955 template<class ImageType>
957 GetMinMaxPointPosition(const ImageType * input,
958 typename ImageType::PointType & min,
959 typename ImageType::PointType & max)
961 typename ImageType::IndexType index = input->GetLargestPossibleRegion().GetIndex();
962 input->TransformIndexToPhysicalPoint(index, min);
963 index = index+input->GetLargestPossibleRegion().GetSize();
964 input->TransformIndexToPhysicalPoint(index, max);
966 //--------------------------------------------------------------------
969 //--------------------------------------------------------------------
970 template<class ImageType>
971 typename ImageType::PointType
972 FindExtremaPointInAGivenLine(const ImageType * input,
975 typename ImageType::PointType p,
976 typename ImageType::PixelType BG,
979 // Which direction ? Increasing or decreasing.
983 // Transform to pixel index
984 typename ImageType::IndexType index;
985 input->TransformPhysicalPointToIndex(p, index);
987 // Loop while inside the mask;
988 while (input->GetPixel(index) != BG) {
989 index[dimension] += d;
992 // Transform back to Physical Units
993 typename ImageType::PointType result;
994 input->TransformIndexToPhysicalPoint(index, result);
996 // Check that is is not too far away
997 double distance = p.EuclideanDistanceTo(result);
998 if (distance > distanceMax) {
999 result = p; // Get back to initial value
1004 //--------------------------------------------------------------------
1007 //--------------------------------------------------------------------
1008 template<class PointType>
1010 IsOnTheSameLineSide(PointType C, PointType A, PointType B, PointType like)
1012 // Look at the position of point 'like' according to the AB line
1013 double s = (B[0] - A[0]) * (like[1] - A[1]) - (B[1] - A[1]) * (like[0] - A[0]);
1014 bool negative = s<0;
1016 // Look the C position
1017 s = (B[0] - A[0]) * (C[1] - A[1]) - (B[1] - A[1]) * (C[0] - A[0]);
1019 if (negative && (s<=0)) return true;
1020 if (!negative && (s>=0)) return true;
1023 //--------------------------------------------------------------------
1026 //--------------------------------------------------------------------
1027 /* Consider an input object, for each slice, find the extrema
1028 position according to a given direction and build a line segment
1029 passing throught this point in a given direction. Output is a
1030 vector of line (from point A to B), for each slice;
1032 template<class ImageType>
1034 SliceBySliceBuildLineSegmentAccordingToExtremaPosition(const ImageType * input,
1035 typename ImageType::PixelType BG,
1037 int extremaDirection,
1038 bool extremaOppositeFlag,
1041 std::vector<typename ImageType::PointType> & A,
1042 std::vector<typename ImageType::PointType> & B)
1045 typedef typename itk::Image<typename ImageType::PixelType, ImageType::ImageDimension-1> SliceType;
1047 // Build the list of slices
1048 std::vector<typename SliceType::Pointer> slices;
1049 clitk::ExtractSlices<ImageType>(input, sliceDimension, slices);
1051 // Build the list of 2D points
1052 std::map<int, typename SliceType::PointType> position2D;
1053 for(uint i=0; i<slices.size(); i++) {
1054 typename SliceType::PointType p;
1056 clitk::FindExtremaPointInAGivenDirection<SliceType>(slices[i], BG,
1057 extremaDirection, extremaOppositeFlag, p);
1063 // Convert 2D points in slice into 3D points
1064 clitk::PointsUtils<ImageType>::Convert2DMapTo3DList(position2D, input, A);
1066 // Create additional point just right to the previous ones, on the
1067 // given lineDirection, in order to create a horizontal/vertical line.
1068 for(uint i=0; i<A.size(); i++) {
1069 typename ImageType::PointType p = A[i];
1070 p[lineDirection] += 10;
1073 A[i][extremaDirection] += margin;
1074 B[i][extremaDirection] += margin;
1078 //--------------------------------------------------------------------
1081 //--------------------------------------------------------------------
1082 template<class ImageType>
1083 typename ImageType::Pointer
1084 SliceBySliceKeepMainCCL(const ImageType * input,
1085 typename ImageType::PixelType BG,
1086 typename ImageType::PixelType FG) {
1089 const int d = ImageType::ImageDimension-1;
1090 typedef typename itk::Image<typename ImageType::PixelType, d> SliceType;
1091 std::vector<typename SliceType::Pointer> slices;
1092 clitk::ExtractSlices<ImageType>(input, d, slices);
1094 // Labelize and keep the main one
1095 std::vector<typename SliceType::Pointer> o;
1096 for(uint i=0; i<slices.size(); i++) {
1097 o.push_back(clitk::Labelize<SliceType>(slices[i], BG, false, 1));
1098 o[i] = clitk::KeepLabels<SliceType>(o[i], BG, FG, 1, 1, true);
1102 typename ImageType::Pointer output;
1103 output = clitk::JoinSlices<ImageType>(o, input, d);
1106 //--------------------------------------------------------------------
1109 //--------------------------------------------------------------------
1110 template<class ImageType>
1111 typename ImageType::Pointer
1112 Clone(const ImageType * input) {
1113 typedef itk::ImageDuplicator<ImageType> DuplicatorType;
1114 typename DuplicatorType::Pointer duplicator = DuplicatorType::New();
1115 duplicator->SetInputImage(input);
1116 duplicator->Update();
1117 return duplicator->GetOutput();
1119 //--------------------------------------------------------------------
1122 //--------------------------------------------------------------------
1123 /* Consider an input object, start at A, for each slice (dim1):
1124 - compute the intersection between the AB line and the current slice
1125 - remove what is at lower or greater according to dim2 of this point
1128 template<class ImageType>
1129 typename ImageType::Pointer
1130 SliceBySliceSetBackgroundFromSingleLine(const ImageType * input,
1131 typename ImageType::PixelType BG,
1132 typename ImageType::PointType & A,
1133 typename ImageType::PointType & B,
1134 int dim1, int dim2, bool removeLowerPartFlag)
1138 typedef typename itk::Image<typename ImageType::PixelType, ImageType::ImageDimension-1> SliceType;
1139 typedef typename SliceType::Pointer SlicePointer;
1140 std::vector<SlicePointer> slices;
1141 clitk::ExtractSlices<ImageType>(input, dim1, slices);
1143 // Start at slice that contains A, and stop at B
1144 typename ImageType::IndexType Ap;
1145 typename ImageType::IndexType Bp;
1146 input->TransformPhysicalPointToIndex(A, Ap);
1147 input->TransformPhysicalPointToIndex(B, Bp);
1149 // Determine slice largest region
1150 typename SliceType::RegionType region = slices[0]->GetLargestPossibleRegion();
1151 typename SliceType::SizeType size = region.GetSize();
1152 typename SliceType::IndexType index = region.GetIndex();
1155 double a = (Bp[dim2]-Ap[dim2])/(Bp[dim1]-Ap[dim1]);
1156 double b = Ap[dim2];
1158 // Loop from slice A to slice B
1159 for(uint i=0; i<(Bp[dim1]-Ap[dim1]); i++) {
1160 // Compute intersection between line AB and current slice for the dim2
1162 // Change region (lower than dim2)
1163 if (removeLowerPartFlag) {
1164 size[dim2] = p-Ap[dim2];
1167 size[dim2] = slices[0]->GetLargestPossibleRegion().GetSize()[dim2]-p;
1170 region.SetSize(size);
1171 region.SetIndex(index);
1172 // Fill region with BG (simple region iterator)
1173 FillRegionWithValue<SliceType>(slices[i+Ap[dim1]], BG, region);
1175 typedef itk::ImageRegionIterator<SliceType> IteratorType;
1176 IteratorType iter(slices[i+Ap[dim1]], region);
1178 while (!iter.IsAtEnd()) {
1187 typename ImageType::Pointer output;
1188 output = clitk::JoinSlices<ImageType>(slices, input, dim1);
1191 //--------------------------------------------------------------------
1193 //--------------------------------------------------------------------
1194 /* Consider an input object, slice by slice, use the point A and set
1195 pixel to BG according to their position relatively to A
1197 template<class ImageType>
1198 typename ImageType::Pointer
1199 SliceBySliceSetBackgroundFromPoints(const ImageType * input,
1200 typename ImageType::PixelType BG,
1202 std::vector<typename ImageType::PointType> & A,
1203 bool removeGreaterThanXFlag,
1204 bool removeGreaterThanYFlag)
1208 typedef typename itk::Image<typename ImageType::PixelType, ImageType::ImageDimension-1> SliceType;
1209 typedef typename SliceType::Pointer SlicePointer;
1210 std::vector<SlicePointer> slices;
1211 clitk::ExtractSlices<ImageType>(input, sliceDim, slices);
1213 // Start at slice that contains A
1214 typename ImageType::IndexType Ap;
1216 // Determine slice largest region
1217 typename SliceType::RegionType region = slices[0]->GetLargestPossibleRegion();
1218 typename SliceType::SizeType size = region.GetSize();
1219 typename SliceType::IndexType index = region.GetIndex();
1221 // Loop from slice A to slice B
1222 for(uint i=0; i<A.size(); i++) {
1223 input->TransformPhysicalPointToIndex(A[i], Ap);
1224 uint sliceIndex = Ap[2] - input->GetLargestPossibleRegion().GetIndex()[2];
1225 if ((sliceIndex < 0) || (sliceIndex >= slices.size())) {
1226 continue; // do not consider this slice
1229 // Compute region for BG
1230 if (removeGreaterThanXFlag) {
1232 size[0] = region.GetSize()[0]-(index[0]-region.GetIndex()[0]);
1235 index[0] = region.GetIndex()[0];
1236 size[0] = Ap[0] - index[0];
1239 if (removeGreaterThanYFlag) {
1241 size[1] = region.GetSize()[1]-(index[1]-region.GetIndex()[1]);
1244 index[1] = region.GetIndex()[1];
1245 size[1] = Ap[1] - index[1];
1249 region.SetSize(size);
1250 region.SetIndex(index);
1252 // Fill region with BG (simple region iterator)
1253 FillRegionWithValue<SliceType>(slices[sliceIndex], BG, region);
1258 typename ImageType::Pointer output;
1259 output = clitk::JoinSlices<ImageType>(slices, input, sliceDim);
1262 //--------------------------------------------------------------------
1265 //--------------------------------------------------------------------
1266 template<class ImageType>
1268 FillRegionWithValue(ImageType * input, typename ImageType::PixelType value, typename ImageType::RegionType & region)
1270 typedef itk::ImageRegionIterator<ImageType> IteratorType;
1271 IteratorType iter(input, region);
1273 while (!iter.IsAtEnd()) {
1278 //--------------------------------------------------------------------
1281 //--------------------------------------------------------------------
1282 template<class ImageType>
1284 GetMinMaxBoundary(ImageType * input, typename ImageType::PointType & min,
1285 typename ImageType::PointType & max)
1287 typedef typename ImageType::PointType PointType;
1288 typedef typename ImageType::IndexType IndexType;
1289 IndexType min_i, max_i;
1290 min_i = input->GetLargestPossibleRegion().GetIndex();
1291 for(uint i=0; i<ImageType::ImageDimension; i++)
1292 max_i[i] = input->GetLargestPossibleRegion().GetSize()[i] + min_i[i];
1293 input->TransformIndexToPhysicalPoint(min_i, min);
1294 input->TransformIndexToPhysicalPoint(max_i, max);
1296 //--------------------------------------------------------------------
1299 //--------------------------------------------------------------------
1300 template<class ImageType>
1301 typename itk::Image<float, ImageType::ImageDimension>::Pointer
1302 DistanceMap(const ImageType * input, typename ImageType::PixelType BG)//,
1303 // typename itk::Image<float, ImageType::ImageDimension>::Pointer dmap)
1305 typedef itk::Image<float,ImageType::ImageDimension> FloatImageType;
1306 typedef itk::SignedMaurerDistanceMapImageFilter<ImageType, FloatImageType> DistanceMapFilterType;
1307 typename DistanceMapFilterType::Pointer filter = DistanceMapFilterType::New();
1308 filter->SetInput(input);
1309 filter->SetUseImageSpacing(true);
1310 filter->SquaredDistanceOff();
1311 filter->SetBackgroundValue(BG);
1313 return filter->GetOutput();
1315 //--------------------------------------------------------------------
1318 //--------------------------------------------------------------------
1319 template<class ImageType>
1321 SliceBySliceBuildLineSegmentAccordingToMinimalDistanceBetweenStructures(const ImageType * S1,
1322 const ImageType * S2,
1323 typename ImageType::PixelType BG,
1325 std::vector<typename ImageType::PointType> & A,
1326 std::vector<typename ImageType::PointType> & B)
1329 typedef typename itk::Image<typename ImageType::PixelType, 2> SliceType;
1330 typedef typename SliceType::Pointer SlicePointer;
1331 std::vector<SlicePointer> slices_s1;
1332 std::vector<SlicePointer> slices_s2;
1333 clitk::ExtractSlices<ImageType>(S1, sliceDimension, slices_s1);
1334 clitk::ExtractSlices<ImageType>(S2, sliceDimension, slices_s2);
1336 assert(slices_s1.size() == slices_s2.size());
1339 typedef itk::Image<float,2> FloatImageType;
1340 typedef itk::SignedMaurerDistanceMapImageFilter<SliceType, FloatImageType> DistanceMapFilterType;
1341 std::vector<typename FloatImageType::Pointer> dmaps1;
1342 std::vector<typename FloatImageType::Pointer> dmaps2;
1343 typename FloatImageType::Pointer dmap;
1346 for(uint i=0; i<slices_s1.size(); i++) {
1347 // Compute dmap for S1 *TO PUT IN FONCTION*
1348 dmap = clitk::DistanceMap<SliceType>(slices_s1[i], BG);
1349 dmaps1.push_back(dmap);
1350 //writeImage<FloatImageType>(dmap, "dmap1.mha");
1351 // Compute dmap for S2
1352 dmap = clitk::DistanceMap<SliceType>(slices_s2[i], BG);
1353 dmaps2.push_back(dmap);
1354 //writeImage<FloatImageType>(dmap, "dmap2.mha");
1356 // Look in S2 for the point the closest to S1
1357 typename SliceType::PointType p = ComputeClosestPoint<SliceType>(slices_s1[i], dmaps2[i], BG);
1358 typename ImageType::PointType p3D;
1359 clitk::PointsUtils<ImageType>::Convert2DTo3D(p, S1, i, p3D);
1362 // Look in S2 for the point the closest to S1
1363 p = ComputeClosestPoint<SliceType>(slices_s2[i], dmaps1[i], BG);
1364 clitk::PointsUtils<ImageType>::Convert2DTo3D(p, S2, i, p3D);
1371 typedef itk::Image<float,3> FT;
1372 FT::Pointer f = FT::New();
1373 typename FT::Pointer d1 = clitk::JoinSlices<FT>(dmaps1, S1, 2);
1374 typename FT::Pointer d2 = clitk::JoinSlices<FT>(dmaps2, S2, 2);
1375 writeImage<FT>(d1, "d1.mha");
1376 writeImage<FT>(d2, "d2.mha");
1379 //--------------------------------------------------------------------
1382 //--------------------------------------------------------------------
1383 template<class ImageType>
1384 typename ImageType::PointType
1385 ComputeClosestPoint(const ImageType * input,
1386 const itk::Image<float, ImageType::ImageDimension> * dmap,
1387 typename ImageType::PixelType & BG)
1389 // Loop dmap + S2, if FG, get min
1390 typedef itk::Image<float,ImageType::ImageDimension> FloatImageType;
1391 typedef itk::ImageRegionConstIteratorWithIndex<ImageType> ImageIteratorType;
1392 typedef itk::ImageRegionConstIterator<FloatImageType> DMapIteratorType;
1393 ImageIteratorType iter1(input, input->GetLargestPossibleRegion());
1394 DMapIteratorType iter2(dmap, dmap->GetLargestPossibleRegion());
1398 double dmin = 100000.0;
1399 typename ImageType::IndexType indexmin;
1401 while (!iter1.IsAtEnd()) {
1402 if (iter1.Get() != BG) {
1403 double d = iter2.Get();
1405 indexmin = iter1.GetIndex();
1414 typename ImageType::PointType p;
1415 input->TransformIndexToPhysicalPoint(indexmin, p);
1418 //--------------------------------------------------------------------
1421 //--------------------------------------------------------------------
1422 template<class ImageType>
1423 typename ImageType::Pointer
1424 RemoveNegativeIndexFromRegion(ImageType * input) {
1425 typedef itk::ChangeInformationImageFilter< ImageType > InfoFilterType;
1426 typename InfoFilterType::Pointer indexChangeFilter = InfoFilterType::New();
1427 indexChangeFilter->ChangeRegionOn();
1428 // The next line is commented because not exist in itk 3
1429 // typename InfoFilterType::OutputImageOffsetValueType indexShift[3];
1431 typename ImageType::IndexType index = input->GetLargestPossibleRegion().GetIndex();
1432 for(uint i=0;i<ImageType::ImageDimension; i++)
1433 indexShift[i] = (index[i]<0 ? -index[i]:0);
1434 typename ImageType::PointType origin;
1435 for(uint i=0;i<ImageType::ImageDimension; i++)
1436 origin[i] = input->GetOrigin()[i] - indexShift[i]*input->GetSpacing()[i];
1437 indexChangeFilter->SetOutputOffset( indexShift );
1438 indexChangeFilter->SetInput(input);
1439 indexChangeFilter->SetOutputOrigin(origin);
1440 indexChangeFilter->ChangeOriginOn();
1441 indexChangeFilter->Update();
1442 return indexChangeFilter->GetOutput();
1444 //--------------------------------------------------------------------
1447 } // end of namespace