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>
43 void ComputeBBFromImageRegion(const ImageType * image,
44 typename ImageType::RegionType region,
45 typename itk::BoundingBox<unsigned long,
46 ImageType::ImageDimension>::Pointer bb) {
47 typedef typename ImageType::IndexType IndexType;
50 for(unsigned int i=0; i<image->GetImageDimension(); i++) {
51 firstIndex[i] = region.GetIndex()[i];
52 lastIndex[i] = firstIndex[i]+region.GetSize()[i];
55 typedef itk::BoundingBox<unsigned long,
56 ImageType::ImageDimension> BBType;
57 typedef typename BBType::PointType PointType;
60 image->TransformIndexToPhysicalPoint(firstIndex, firstPoint);
61 image->TransformIndexToPhysicalPoint(lastIndex, lastPoint);
63 bb->SetMaximum(lastPoint);
64 bb->SetMinimum(firstPoint);
66 //--------------------------------------------------------------------
69 //--------------------------------------------------------------------
70 template<int Dimension>
71 void ComputeBBIntersection(typename itk::BoundingBox<unsigned long, Dimension>::Pointer bbo,
72 typename itk::BoundingBox<unsigned long, Dimension>::Pointer bbi1,
73 typename itk::BoundingBox<unsigned long, Dimension>::Pointer bbi2) {
75 typedef itk::BoundingBox<unsigned long, Dimension> BBType;
76 typedef typename BBType::PointType PointType;
80 for(unsigned int i=0; i<Dimension; i++) {
81 firstPoint[i] = std::max(bbi1->GetMinimum()[i],
82 bbi2->GetMinimum()[i]);
83 lastPoint[i] = std::min(bbi1->GetMaximum()[i],
84 bbi2->GetMaximum()[i]);
87 bbo->SetMaximum(lastPoint);
88 bbo->SetMinimum(firstPoint);
90 //--------------------------------------------------------------------
93 //--------------------------------------------------------------------
94 template<class ImageType>
95 void ComputeRegionFromBB(const ImageType * image,
96 const typename itk::BoundingBox<unsigned long,
97 ImageType::ImageDimension>::Pointer bb,
98 typename ImageType::RegionType & region) {
100 typedef typename ImageType::IndexType IndexType;
101 typedef typename ImageType::PointType PointType;
102 typedef typename ImageType::RegionType RegionType;
103 typedef typename ImageType::SizeType SizeType;
105 // Region starting point
106 IndexType regionStart;
107 PointType start = bb->GetMinimum();
108 image->TransformPhysicalPointToIndex(start, regionStart);
112 PointType maxs = bb->GetMaximum();
113 PointType mins = bb->GetMinimum();
114 for(unsigned int i=0; i<ImageType::ImageDimension; i++) {
115 regionSize[i] = lrint((maxs[i] - mins[i])/image->GetSpacing()[i]);
119 region.SetIndex(regionStart);
120 region.SetSize(regionSize);
122 //--------------------------------------------------------------------
124 //--------------------------------------------------------------------
125 template<class ImageType, class TMaskImageType>
126 typename ImageType::Pointer
127 SetBackground(const ImageType * input,
128 const TMaskImageType * mask,
129 typename TMaskImageType::PixelType maskBG,
130 typename ImageType::PixelType outValue,
132 typedef SetBackgroundImageFilter<ImageType, TMaskImageType, ImageType>
133 SetBackgroundImageFilterType;
134 typename SetBackgroundImageFilterType::Pointer setBackgroundFilter
135 = SetBackgroundImageFilterType::New();
136 // if (inPlace) setBackgroundFilter->ReleaseDataFlagOn(); // No seg fault
137 setBackgroundFilter->SetInPlace(inPlace); // This is important to keep memory low
138 setBackgroundFilter->SetInput(input);
139 setBackgroundFilter->SetInput2(mask);
140 setBackgroundFilter->SetMaskValue(maskBG);
141 setBackgroundFilter->SetOutsideValue(outValue);
142 setBackgroundFilter->Update();
143 return setBackgroundFilter->GetOutput();
145 //--------------------------------------------------------------------
148 //--------------------------------------------------------------------
149 template<class ImageType>
150 int GetNumberOfConnectedComponentLabels(const ImageType * input,
151 typename ImageType::PixelType BG,
152 bool isFullyConnected) {
153 // Connected Component label
154 typedef itk::ConnectedComponentImageFilter<ImageType, ImageType> ConnectFilterType;
155 typename ConnectFilterType::Pointer connectFilter = ConnectFilterType::New();
156 connectFilter->SetInput(input);
157 connectFilter->SetBackgroundValue(BG);
158 connectFilter->SetFullyConnected(isFullyConnected);
159 connectFilter->Update();
162 return connectFilter->GetObjectCount();
164 //--------------------------------------------------------------------
166 //--------------------------------------------------------------------
168 Warning : in this cas, we consider outputType like inputType, not
169 InternalImageType. Be sure it fits.
171 template<class ImageType>
172 typename ImageType::Pointer
173 Labelize(const ImageType * input,
174 typename ImageType::PixelType BG,
175 bool isFullyConnected,
176 int minimalComponentSize) {
177 // InternalImageType for storing large number of component
178 typedef itk::Image<int, ImageType::ImageDimension> InternalImageType;
180 // Connected Component label
181 typedef itk::ConnectedComponentImageFilter<ImageType, InternalImageType> ConnectFilterType;
182 typename ConnectFilterType::Pointer connectFilter = ConnectFilterType::New();
183 // connectFilter->ReleaseDataFlagOn();
184 connectFilter->SetInput(input);
185 connectFilter->SetBackgroundValue(BG);
186 connectFilter->SetFullyConnected(isFullyConnected);
188 // Sort by size and remove too small area.
189 typedef itk::RelabelComponentImageFilter<InternalImageType, ImageType> RelabelFilterType;
190 typename RelabelFilterType::Pointer relabelFilter = RelabelFilterType::New();
191 // relabelFilter->ReleaseDataFlagOn(); // if yes, fail when ExplosionControlledThresholdConnectedImageFilter ???
192 relabelFilter->SetInput(connectFilter->GetOutput());
193 relabelFilter->SetMinimumObjectSize(minimalComponentSize);
194 relabelFilter->Update();
197 typename ImageType::Pointer output = relabelFilter->GetOutput();
200 //--------------------------------------------------------------------
203 //--------------------------------------------------------------------
205 Warning : in this cas, we consider outputType like inputType, not
206 InternalImageType. Be sure it fits.
208 template<class ImageType>
209 typename ImageType::Pointer
210 LabelizeAndCountNumberOfObjects(const ImageType * input,
211 typename ImageType::PixelType BG,
212 bool isFullyConnected,
213 int minimalComponentSize,
215 // InternalImageType for storing large number of component
216 typedef itk::Image<int, ImageType::ImageDimension> InternalImageType;
218 // Connected Component label
219 typedef itk::ConnectedComponentImageFilter<ImageType, InternalImageType> ConnectFilterType;
220 typename ConnectFilterType::Pointer connectFilter = ConnectFilterType::New();
221 // connectFilter->ReleaseDataFlagOn();
222 connectFilter->SetInput(input);
223 connectFilter->SetBackgroundValue(BG);
224 connectFilter->SetFullyConnected(isFullyConnected);
226 // Sort by size and remove too small area.
227 typedef itk::RelabelComponentImageFilter<InternalImageType, ImageType> RelabelFilterType;
228 typename RelabelFilterType::Pointer relabelFilter = RelabelFilterType::New();
229 // relabelFilter->ReleaseDataFlagOn(); // if yes, fail when ExplosionControlledThresholdConnectedImageFilter ???
230 relabelFilter->SetInput(connectFilter->GetOutput());
231 relabelFilter->SetMinimumObjectSize(minimalComponentSize);
232 relabelFilter->Update();
234 nb = relabelFilter->GetNumberOfObjects();
235 // DD(relabelFilter->GetOriginalNumberOfObjects());
236 // DD(relabelFilter->GetSizeOfObjectsInPhysicalUnits()[0]);
239 typename ImageType::Pointer output = relabelFilter->GetOutput();
242 //--------------------------------------------------------------------
246 //--------------------------------------------------------------------
247 template<class ImageType>
248 typename ImageType::Pointer
249 RemoveLabels(const ImageType * input,
250 typename ImageType::PixelType BG,
251 std::vector<typename ImageType::PixelType> & labelsToRemove) {
252 assert(labelsToRemove.size() != 0);
253 typename ImageType::Pointer working_image;// = input;
254 for (unsigned int i=0; i <labelsToRemove.size(); i++) {
255 typedef SetBackgroundImageFilter<ImageType, ImageType> SetBackgroundImageFilterType;
256 typename SetBackgroundImageFilterType::Pointer setBackgroundFilter = SetBackgroundImageFilterType::New();
257 setBackgroundFilter->SetInput(input);
258 setBackgroundFilter->SetInput2(input);
259 setBackgroundFilter->SetMaskValue(labelsToRemove[i]);
260 setBackgroundFilter->SetOutsideValue(BG);
261 setBackgroundFilter->Update();
262 working_image = setBackgroundFilter->GetOutput();
264 return working_image;
266 //--------------------------------------------------------------------
269 //--------------------------------------------------------------------
270 template<class ImageType>
271 typename ImageType::Pointer
272 KeepLabels(const ImageType * input,
273 typename ImageType::PixelType BG,
274 typename ImageType::PixelType FG,
275 typename ImageType::PixelType firstKeep,
276 typename ImageType::PixelType lastKeep,
278 typedef itk::BinaryThresholdImageFilter<ImageType, ImageType> BinarizeFilterType;
279 typename BinarizeFilterType::Pointer binarizeFilter = BinarizeFilterType::New();
280 binarizeFilter->SetInput(input);
281 binarizeFilter->SetLowerThreshold(firstKeep);
282 if (useLastKeep) binarizeFilter->SetUpperThreshold(lastKeep);
283 binarizeFilter->SetInsideValue(FG);
284 binarizeFilter->SetOutsideValue(BG);
285 binarizeFilter->Update();
286 return binarizeFilter->GetOutput();
288 //--------------------------------------------------------------------
291 //--------------------------------------------------------------------
292 template<class ImageType>
293 typename ImageType::Pointer
294 LabelizeAndSelectLabels(const ImageType * input,
295 typename ImageType::PixelType BG,
296 typename ImageType::PixelType FG,
297 bool isFullyConnected,
298 int minimalComponentSize,
299 LabelizeParameters<typename ImageType::PixelType> * param)
301 typename ImageType::Pointer working_image;
302 working_image = Labelize<ImageType>(input, BG, isFullyConnected, minimalComponentSize);
303 if (param->GetLabelsToRemove().size() != 0)
304 working_image = RemoveLabels<ImageType>(working_image, BG, param->GetLabelsToRemove());
305 working_image = KeepLabels<ImageType>(working_image,
307 param->GetFirstKeep(),
308 param->GetLastKeep(),
309 param->GetUseLastKeep());
310 return working_image;
312 //--------------------------------------------------------------------
315 //--------------------------------------------------------------------
316 template<class ImageType>
317 typename ImageType::Pointer
318 ResizeImageLike(const ImageType * input,
319 const itk::ImageBase<ImageType::ImageDimension> * like,
320 typename ImageType::PixelType backgroundValue)
322 typedef CropLikeImageFilter<ImageType> CropFilterType;
323 typename CropFilterType::Pointer cropFilter = CropFilterType::New();
324 cropFilter->SetInput(input);
325 cropFilter->SetCropLikeImage(like);
326 cropFilter->SetBackgroundValue(backgroundValue);
327 cropFilter->Update();
328 return cropFilter->GetOutput();
330 //--------------------------------------------------------------------
333 //--------------------------------------------------------------------
334 template<class MaskImageType>
335 typename MaskImageType::Pointer
336 SliceBySliceRelativePosition(const MaskImageType * input,
337 const MaskImageType * object,
340 std::string orientation,
341 bool uniqueConnectedComponent,
344 bool singleObjectCCL)
346 typedef SliceBySliceRelativePositionFilter<MaskImageType> SliceRelPosFilterType;
347 typename SliceRelPosFilterType::Pointer sliceRelPosFilter = SliceRelPosFilterType::New();
348 sliceRelPosFilter->VerboseStepFlagOff();
349 sliceRelPosFilter->WriteStepFlagOff();
350 sliceRelPosFilter->SetInput(input);
351 sliceRelPosFilter->SetInputObject(object);
352 sliceRelPosFilter->SetDirection(direction);
353 sliceRelPosFilter->SetFuzzyThreshold(threshold);
354 sliceRelPosFilter->AddOrientationTypeString(orientation);
355 sliceRelPosFilter->SetIntermediateSpacingFlag((spacing != -1));
356 sliceRelPosFilter->SetIntermediateSpacing(spacing);
357 sliceRelPosFilter->SetUniqueConnectedComponentBySliceFlag(uniqueConnectedComponent);
358 sliceRelPosFilter->ObjectCCLSelectionFlagOff();
359 sliceRelPosFilter->SetUseTheLargestObjectCCLFlag(singleObjectCCL);
360 // sliceRelPosFilter->SetInverseOrientationFlag(inverseflag);
361 sliceRelPosFilter->SetAutoCropFlag(autocropFlag);
362 sliceRelPosFilter->IgnoreEmptySliceObjectFlagOn();
363 sliceRelPosFilter->Update();
364 return sliceRelPosFilter->GetOutput();
366 //--------------------------------------------------------------------
369 //--------------------------------------------------------------------
370 template<class ImageType>
372 FindExtremaPointInAGivenDirection(const ImageType * input,
373 typename ImageType::PixelType bg,
374 int direction, bool opposite,
375 typename ImageType::PointType & point)
377 typename ImageType::PointType dummy;
378 return FindExtremaPointInAGivenDirection(input, bg, direction,
379 opposite, dummy, 0, point);
381 //--------------------------------------------------------------------
384 //--------------------------------------------------------------------
385 template<class ImageType>
387 FindExtremaPointInAGivenDirection(const ImageType * input,
388 typename ImageType::PixelType bg,
389 int direction, bool opposite,
390 typename ImageType::PointType refpoint,
392 typename ImageType::PointType & point)
395 loop over input pixels, store the index in the fg that is max
396 according to the given direction.
398 typedef itk::ImageRegionConstIteratorWithIndex<ImageType> IteratorType;
399 IteratorType iter(input, input->GetLargestPossibleRegion());
401 typename ImageType::IndexType max = input->GetLargestPossibleRegion().GetIndex();
402 if (opposite) max = max+input->GetLargestPossibleRegion().GetSize();
404 while (!iter.IsAtEnd()) {
405 if (iter.Get() != bg) {
406 bool test = iter.GetIndex()[direction] > max[direction];
407 if (opposite) test = !test;
409 typename ImageType::PointType p;
410 input->TransformIndexToPhysicalPoint(iter.GetIndex(), p);
411 if ((distanceMax==0) || (p.EuclideanDistanceTo(refpoint) < distanceMax)) {
412 max = iter.GetIndex();
419 if (!found) return false;
420 input->TransformIndexToPhysicalPoint(max, point);
423 //--------------------------------------------------------------------
426 //--------------------------------------------------------------------
427 template<class ImageType>
428 typename ImageType::Pointer
429 CropImageRemoveGreaterThan(const ImageType * image,
430 int dim, double min, bool autoCrop,
431 typename ImageType::PixelType BG)
433 return CropImageAlongOneAxis<ImageType>(image, dim,
434 image->GetOrigin()[dim],
438 //--------------------------------------------------------------------
441 //--------------------------------------------------------------------
442 template<class ImageType>
443 typename ImageType::Pointer
444 CropImageRemoveLowerThan(const ImageType * image,
445 int dim, double max, bool autoCrop,
446 typename ImageType::PixelType BG)
448 typename ImageType::PointType p;
449 image->TransformIndexToPhysicalPoint(image->GetLargestPossibleRegion().GetIndex()+
450 image->GetLargestPossibleRegion().GetSize(), p);
451 return CropImageAlongOneAxis<ImageType>(image, dim, max, p[dim], autoCrop, BG);
453 //--------------------------------------------------------------------
456 //--------------------------------------------------------------------
457 template<class ImageType>
458 typename ImageType::Pointer
459 CropImageAlongOneAxis(const ImageType * image,
460 int dim, double min, double max,
461 bool autoCrop, typename ImageType::PixelType BG)
463 // Compute region size
464 typename ImageType::RegionType region;
465 typename ImageType::SizeType size = image->GetLargestPossibleRegion().GetSize();
466 typename ImageType::PointType p = image->GetOrigin();
467 if (min > p[dim]) p[dim] = min; // Check if not outside the image
468 typename ImageType::IndexType start;
469 image->TransformPhysicalPointToIndex(p, start);
470 double m = image->GetOrigin()[dim] + size[dim]*image->GetSpacing()[dim];
471 if (max > m) p[dim] = m; // Check if not outside the image
473 typename ImageType::IndexType end;
474 image->TransformPhysicalPointToIndex(p, end);
475 size[dim] = abs(end[dim]-start[dim]);
476 region.SetIndex(start);
477 region.SetSize(size);
480 typedef itk::RegionOfInterestImageFilter<ImageType, ImageType> CropFilterType;
481 typename CropFilterType::Pointer cropFilter = CropFilterType::New();
482 cropFilter->SetInput(image);
483 cropFilter->SetRegionOfInterest(region);
484 cropFilter->Update();
485 typename ImageType::Pointer result = cropFilter->GetOutput();
489 result = AutoCrop<ImageType>(result, BG);
493 //--------------------------------------------------------------------
496 //--------------------------------------------------------------------
497 template<class ImageType>
499 ComputeCentroids(const ImageType * image,
500 typename ImageType::PixelType BG,
501 std::vector<typename ImageType::PointType> & centroids)
503 typedef long LabelType;
504 static const unsigned int Dim = ImageType::ImageDimension;
505 typedef itk::ShapeLabelObject< LabelType, Dim > LabelObjectType;
506 typedef itk::LabelMap< LabelObjectType > LabelMapType;
507 typedef itk::LabelImageToLabelMapFilter<ImageType, LabelMapType> ImageToMapFilterType;
508 typename ImageToMapFilterType::Pointer imageToLabelFilter = ImageToMapFilterType::New();
509 typedef itk::ShapeLabelMapFilter<LabelMapType, ImageType> ShapeFilterType;
510 typename ShapeFilterType::Pointer statFilter = ShapeFilterType::New();
511 imageToLabelFilter->SetBackgroundValue(BG);
512 imageToLabelFilter->SetInput(image);
513 statFilter->SetInput(imageToLabelFilter->GetOutput());
514 statFilter->Update();
515 typename LabelMapType::Pointer labelMap = statFilter->GetOutput();
518 typename ImageType::PointType dummy;
519 centroids.push_back(dummy); // label 0 -> no centroid, use dummy point for BG
520 //DS FIXME (not useful ! to change ..)
521 for(uint i=0; i<labelMap->GetNumberOfLabelObjects(); i++) {
522 int label = labelMap->GetLabels()[i];
523 centroids.push_back(labelMap->GetLabelObject(label)->GetCentroid());
526 //--------------------------------------------------------------------
529 //--------------------------------------------------------------------
530 template<class ImageType, class LabelType>
531 typename itk::LabelMap< itk::ShapeLabelObject<LabelType, ImageType::ImageDimension> >::Pointer
532 ComputeLabelMap(const ImageType * image,
533 typename ImageType::PixelType BG,
534 bool computePerimeterFlag)
536 static const unsigned int Dim = ImageType::ImageDimension;
537 typedef itk::ShapeLabelObject< LabelType, Dim > LabelObjectType;
538 typedef itk::LabelMap< LabelObjectType > LabelMapType;
539 typedef itk::LabelImageToLabelMapFilter<ImageType, LabelMapType> ImageToMapFilterType;
540 typename ImageToMapFilterType::Pointer imageToLabelFilter = ImageToMapFilterType::New();
541 typedef itk::ShapeLabelMapFilter<LabelMapType, ImageType> ShapeFilterType;
542 typename ShapeFilterType::Pointer statFilter = ShapeFilterType::New();
543 imageToLabelFilter->SetBackgroundValue(BG);
544 imageToLabelFilter->SetInput(image);
545 statFilter->SetInput(imageToLabelFilter->GetOutput());
546 statFilter->SetComputePerimeter(computePerimeterFlag);
547 statFilter->Update();
548 return statFilter->GetOutput();
550 //--------------------------------------------------------------------
553 //--------------------------------------------------------------------
554 template<class ImageType>
556 ComputeCentroids2(const ImageType * image,
557 typename ImageType::PixelType BG,
558 std::vector<typename ImageType::PointType> & centroids)
560 typedef long LabelType;
561 static const unsigned int Dim = ImageType::ImageDimension;
562 typedef itk::ShapeLabelObject< LabelType, Dim > LabelObjectType;
563 typedef itk::LabelMap< LabelObjectType > LabelMapType;
564 typedef itk::LabelImageToLabelMapFilter<ImageType, LabelMapType> ImageToMapFilterType;
565 typename ImageToMapFilterType::Pointer imageToLabelFilter = ImageToMapFilterType::New();
566 typedef itk::ShapeLabelMapFilter<LabelMapType, ImageType> ShapeFilterType;
567 typename ShapeFilterType::Pointer statFilter = ShapeFilterType::New();
568 imageToLabelFilter->SetBackgroundValue(BG);
569 imageToLabelFilter->SetInput(image);
570 statFilter->SetInput(imageToLabelFilter->GetOutput());
571 statFilter->Update();
572 typename LabelMapType::Pointer labelMap = statFilter->GetOutput();
575 typename ImageType::PointType dummy;
576 centroids.push_back(dummy); // label 0 -> no centroid, use dummy point
577 for(uint i=1; i<labelMap->GetNumberOfLabelObjects()+1; i++) {
578 centroids.push_back(labelMap->GetLabelObject(i)->GetCentroid());
581 for(uint i=1; i<labelMap->GetNumberOfLabelObjects()+1; i++) {
582 DD(labelMap->GetLabelObject(i)->GetBinaryPrincipalAxes());
583 DD(labelMap->GetLabelObject(i)->GetBinaryFlatness());
584 DD(labelMap->GetLabelObject(i)->GetRoundness ());
586 // search for the point on the boundary alog PA
591 //--------------------------------------------------------------------
594 //--------------------------------------------------------------------
595 template<class ImageType>
597 PointsUtils<ImageType>::Convert2DTo3D(const PointType2D & p2D,
598 const ImageType * image,
603 index3D[0] = index3D[1] = 0;
604 index3D[2] = image->GetLargestPossibleRegion().GetIndex()[2]+slice;
605 image->TransformIndexToPhysicalPoint(index3D, p3D);
608 // p3D[2] = p[2];//(image->GetLargestPossibleRegion().GetIndex()[2]+slice)*image->GetSpacing()[2]
609 // + image->GetOrigin()[2];
611 //--------------------------------------------------------------------
614 //--------------------------------------------------------------------
615 template<class ImageType>
617 PointsUtils<ImageType>::Convert2DMapTo3DList(const MapPoint2DType & map,
618 const ImageType * image,
619 VectorPoint3DType & list)
621 typename MapPoint2DType::const_iterator iter = map.begin();
622 while (iter != map.end()) {
624 Convert2DTo3D(iter->second, image, iter->first, p);
629 //--------------------------------------------------------------------
632 //--------------------------------------------------------------------
633 template<class ImageType>
635 PointsUtils<ImageType>::Convert2DListTo3DList(const VectorPoint2DType & p2D,
637 const ImageType * image,
638 VectorPoint3DType & list)
640 for(uint i=0; i<p2D.size(); i++) {
642 Convert2DTo3D(p2D[i], image, slice, p);
646 //--------------------------------------------------------------------
649 //--------------------------------------------------------------------
650 template<class ImageType>
652 WriteListOfLandmarks(std::vector<typename ImageType::PointType> points,
653 std::string filename)
656 openFileForWriting(os, filename);
657 os << "LANDMARKS1" << std::endl;
658 for(uint i=0; i<points.size(); i++) {
659 const typename ImageType::PointType & p = points[i];
660 // Write it in the file
661 os << i << " " << p[0] << " " << p[1] << " " << p[2] << " 0 0 " << std::endl;
665 //--------------------------------------------------------------------
668 //--------------------------------------------------------------------
669 template<class ImageType>
670 typename ImageType::Pointer
671 Dilate(const ImageType * image, double radiusInMM,
672 typename ImageType::PixelType BG,
673 typename ImageType::PixelType FG,
676 typename ImageType::SizeType r;
677 for(uint i=0; i<ImageType::ImageDimension; i++)
678 r[i] = (uint)lrint(radiusInMM/image->GetSpacing()[i]);
679 return Dilate<ImageType>(image, r, BG, FG, extendSupport);
681 //--------------------------------------------------------------------
684 //--------------------------------------------------------------------
685 template<class ImageType>
686 typename ImageType::Pointer
687 Dilate(const ImageType * image, typename ImageType::PointType radiusInMM,
688 typename ImageType::PixelType BG,
689 typename ImageType::PixelType FG,
692 typename ImageType::SizeType r;
693 for(uint i=0; i<ImageType::ImageDimension; i++)
694 r[i] = (uint)lrint(radiusInMM[i]/image->GetSpacing()[i]);
695 return Dilate<ImageType>(image, r, BG, FG, extendSupport);
697 //--------------------------------------------------------------------
700 //--------------------------------------------------------------------
701 template<class ImageType>
702 typename ImageType::Pointer
703 Dilate(const ImageType * image, typename ImageType::SizeType radius,
704 typename ImageType::PixelType BG,
705 typename ImageType::PixelType FG,
708 // Create kernel for dilatation
709 typedef itk::BinaryBallStructuringElement<typename ImageType::PixelType,
710 ImageType::ImageDimension> KernelType;
711 KernelType structuringElement;
712 structuringElement.SetRadius(radius);
713 structuringElement.CreateStructuringElement();
715 typename ImageType::Pointer output;
717 typedef itk::ConstantPadImageFilter<ImageType, ImageType> PadFilterType;
718 typename PadFilterType::Pointer padFilter = PadFilterType::New();
719 padFilter->SetInput(image);
720 typename ImageType::SizeType lower;
721 typename ImageType::SizeType upper;
722 for(uint i=0; i<3; i++) {
723 lower[i] = upper[i] = 2*(radius[i]+1);
725 padFilter->SetPadLowerBound(lower);
726 padFilter->SetPadUpperBound(upper);
728 output = padFilter->GetOutput();
732 typedef itk::BinaryDilateImageFilter<ImageType, ImageType , KernelType> DilateFilterType;
733 typename DilateFilterType::Pointer dilateFilter = DilateFilterType::New();
734 dilateFilter->SetBackgroundValue(BG);
735 dilateFilter->SetForegroundValue(FG);
736 dilateFilter->SetBoundaryToForeground(false);
737 dilateFilter->SetKernel(structuringElement);
738 if (extendSupport) dilateFilter->SetInput(output);
739 else dilateFilter->SetInput(image);
740 dilateFilter->Update();
741 return dilateFilter->GetOutput();
743 //--------------------------------------------------------------------
746 //--------------------------------------------------------------------
747 template<class ImageType>
748 typename ImageType::Pointer
749 Opening(const ImageType * image, typename ImageType::SizeType radius,
750 typename ImageType::PixelType BG,
751 typename ImageType::PixelType FG)
754 typedef itk::BinaryBallStructuringElement<typename ImageType::PixelType,
755 ImageType::ImageDimension> KernelType;
756 KernelType structuringElement;
757 structuringElement.SetRadius(radius);
758 structuringElement.CreateStructuringElement();
761 typedef itk::BinaryMorphologicalOpeningImageFilter<ImageType, ImageType , KernelType> OpeningFilterType;
762 typename OpeningFilterType::Pointer open = OpeningFilterType::New();
763 open->SetInput(image);
764 open->SetBackgroundValue(BG);
765 open->SetForegroundValue(FG);
766 open->SetKernel(structuringElement);
768 return open->GetOutput();
770 //--------------------------------------------------------------------
774 //--------------------------------------------------------------------
775 template<class ValueType, class VectorType>
776 void ConvertOption(std::string optionName, uint given,
777 ValueType * values, VectorType & p,
778 uint dim, bool required)
780 if (required && (given == 0)) {
781 clitkExceptionMacro("The option --" << optionName << " must be set and have 1 or "
782 << dim << " values.");
785 for(uint i=0; i<dim; i++) p[i] = values[0];
789 for(uint i=0; i<dim; i++) p[i] = values[i];
792 if (given == 0) return;
793 clitkExceptionMacro("The option --" << optionName << " must have 1 or "
794 << dim << " values.");
796 //--------------------------------------------------------------------
799 //--------------------------------------------------------------------
801 http://www.gamedev.net/community/forums/topic.asp?topic_id=542870
802 Assuming the points are (Ax,Ay) (Bx,By) and (Cx,Cy), you need to compute:
803 (Bx - Ax) * (Cy - Ay) - (By - Ay) * (Cx - Ax)
804 This will equal zero if the point C is on the line formed by
805 points A and B, and will have a different sign depending on the
806 side. Which side this is depends on the orientation of your (x,y)
807 coordinates, but you can plug test values for A,B and C into this
808 formula to determine whether negative values are to the left or to
810 => to accelerate, start with formula, when change sign -> stop and fill
812 offsetToKeep = is used to determine which side of the line we
813 keep. The point along the mainDirection but 'offsetToKeep' mm away
817 template<class ImageType>
819 SliceBySliceSetBackgroundFromLineSeparation(ImageType * input,
820 std::vector<typename ImageType::PointType> & lA,
821 std::vector<typename ImageType::PointType> & lB,
822 typename ImageType::PixelType BG,
826 assert((mainDirection==0) || (mainDirection==1));
827 typedef itk::ImageSliceIteratorWithIndex<ImageType> SliceIteratorType;
828 SliceIteratorType siter = SliceIteratorType(input,
829 input->GetLargestPossibleRegion());
830 siter.SetFirstDirection(0);
831 siter.SetSecondDirection(1);
834 typename ImageType::PointType A;
835 typename ImageType::PointType B;
836 typename ImageType::PointType C;
837 assert(lA.size() == lB.size());
838 while ((i<lA.size()) && (!siter.IsAtEnd())) {
839 // Check that the current slice correspond to the current point
840 input->TransformIndexToPhysicalPoint(siter.GetIndex(), C);
841 if ((fabs(C[2] - lA[i][2]))>0.01) { // is !equal with a tolerance of 0.01 mm
844 // Define A,B,C points
849 // Check that the line is not a point (A=B)
850 bool p = (A[0] == B[0]) && (A[1] == B[1]);
853 C[mainDirection] += offsetToKeep; // I know I must keep this point
854 double s = (B[0] - A[0]) * (C[1] - A[1]) - (B[1] - A[1]) * (C[0] - A[0]);
855 bool isPositive = s<0;
856 while (!siter.IsAtEndOfSlice()) {
857 while (!siter.IsAtEndOfLine()) {
858 // Very slow, I know ... but image should be very small
859 input->TransformIndexToPhysicalPoint(siter.GetIndex(), C);
860 double s = (B[0] - A[0]) * (C[1] - A[1]) - (B[1] - A[1]) * (C[0] - A[0]);
861 if (s == 0) siter.Set(BG); // on the line, we decide to remove
863 if (s > 0) siter.Set(BG);
866 if (s < 0) siter.Set(BG);
875 } // End of current slice
879 //--------------------------------------------------------------------
882 //--------------------------------------------------------------------
883 template<class ImageType>
885 AndNot(ImageType * input,
886 const ImageType * object,
887 typename ImageType::PixelType BG)
889 typename ImageType::Pointer o;
891 if (!clitk::HaveSameSizeAndSpacing<ImageType, ImageType>(input, object)) {
892 o = clitk::ResizeImageLike<ImageType>(object, input, BG);
896 typedef clitk::BooleanOperatorLabelImageFilter<ImageType> BoolFilterType;
897 typename BoolFilterType::Pointer boolFilter = BoolFilterType::New();
898 boolFilter->InPlaceOn();
899 boolFilter->SetInput1(input);
900 if (resized) boolFilter->SetInput2(o);
901 else boolFilter->SetInput2(object);
902 boolFilter->SetBackgroundValue1(BG);
903 boolFilter->SetBackgroundValue2(BG);
904 boolFilter->SetOperationType(BoolFilterType::AndNot);
905 boolFilter->Update();
907 //--------------------------------------------------------------------
910 //--------------------------------------------------------------------
911 template<class ImageType>
913 And(ImageType * input,
914 const ImageType * object,
915 typename ImageType::PixelType BG)
917 typename ImageType::Pointer o;
919 if (!clitk::HaveSameSizeAndSpacing<ImageType, ImageType>(input, object)) {
920 o = clitk::ResizeImageLike<ImageType>(object, input, BG);
924 typedef clitk::BooleanOperatorLabelImageFilter<ImageType> BoolFilterType;
925 typename BoolFilterType::Pointer boolFilter = BoolFilterType::New();
926 boolFilter->InPlaceOn();
927 boolFilter->SetInput1(input);
928 if (resized) boolFilter->SetInput2(o);
929 else boolFilter->SetInput2(object);
930 boolFilter->SetBackgroundValue1(BG);
931 boolFilter->SetBackgroundValue2(BG);
932 boolFilter->SetOperationType(BoolFilterType::And);
933 boolFilter->Update();
935 //--------------------------------------------------------------------
938 //--------------------------------------------------------------------
939 template<class ImageType>
941 Or(ImageType * input,
942 const ImageType * object,
943 typename ImageType::PixelType BG)
945 typename ImageType::Pointer o;
947 if (!clitk::HaveSameSizeAndSpacing<ImageType, ImageType>(input, object)) {
948 o = clitk::ResizeImageLike<ImageType>(object, input, BG);
952 typedef clitk::BooleanOperatorLabelImageFilter<ImageType> BoolFilterType;
953 typename BoolFilterType::Pointer boolFilter = BoolFilterType::New();
954 boolFilter->InPlaceOn();
955 boolFilter->SetInput1(input);
956 if (resized) boolFilter->SetInput2(o);
957 else boolFilter->SetInput2(object);
958 boolFilter->SetBackgroundValue1(BG);
959 boolFilter->SetBackgroundValue2(BG);
960 boolFilter->SetOperationType(BoolFilterType::Or);
961 boolFilter->Update();
963 //--------------------------------------------------------------------
966 //--------------------------------------------------------------------
967 template<class ImageType>
968 typename ImageType::Pointer
969 Binarize(const ImageType * input,
970 typename ImageType::PixelType lower,
971 typename ImageType::PixelType upper,
972 typename ImageType::PixelType BG,
973 typename ImageType::PixelType FG)
975 typedef itk::BinaryThresholdImageFilter<ImageType, ImageType> BinaryThresholdFilterType;
976 typename BinaryThresholdFilterType::Pointer binarizeFilter = BinaryThresholdFilterType::New();
977 binarizeFilter->SetInput(input);
978 binarizeFilter->InPlaceOff();
979 binarizeFilter->SetLowerThreshold(lower);
980 binarizeFilter->SetUpperThreshold(upper);
981 binarizeFilter->SetInsideValue(FG);
982 binarizeFilter->SetOutsideValue(BG);
983 binarizeFilter->Update();
984 return binarizeFilter->GetOutput();
986 //--------------------------------------------------------------------
989 //--------------------------------------------------------------------
990 template<class ImageType>
992 GetMinMaxPointPosition(const ImageType * input,
993 typename ImageType::PointType & min,
994 typename ImageType::PointType & max)
996 typename ImageType::IndexType index = input->GetLargestPossibleRegion().GetIndex();
997 input->TransformIndexToPhysicalPoint(index, min);
998 index = index+input->GetLargestPossibleRegion().GetSize();
999 input->TransformIndexToPhysicalPoint(index, max);
1001 //--------------------------------------------------------------------
1004 //--------------------------------------------------------------------
1005 template<class ImageType>
1006 typename ImageType::PointType
1007 FindExtremaPointInAGivenLine(const ImageType * input,
1010 typename ImageType::PointType p,
1011 typename ImageType::PixelType BG,
1014 // Which direction ? Increasing or decreasing.
1018 // Transform to pixel index
1019 typename ImageType::IndexType index;
1020 input->TransformPhysicalPointToIndex(p, index);
1022 // Loop while inside the mask;
1023 while (input->GetPixel(index) != BG) {
1024 index[dimension] += d;
1027 // Transform back to Physical Units
1028 typename ImageType::PointType result;
1029 input->TransformIndexToPhysicalPoint(index, result);
1031 // Check that is is not too far away
1032 double distance = p.EuclideanDistanceTo(result);
1033 if (distance > distanceMax) {
1034 result = p; // Get back to initial value
1039 //--------------------------------------------------------------------
1042 //--------------------------------------------------------------------
1043 template<class PointType>
1045 IsOnTheSameLineSide(PointType C, PointType A, PointType B, PointType like)
1047 // Look at the position of point 'like' according to the AB line
1048 double s = (B[0] - A[0]) * (like[1] - A[1]) - (B[1] - A[1]) * (like[0] - A[0]);
1049 bool negative = s<0;
1051 // Look the C position
1052 s = (B[0] - A[0]) * (C[1] - A[1]) - (B[1] - A[1]) * (C[0] - A[0]);
1054 if (negative && (s<=0)) return true;
1055 if (!negative && (s>=0)) return true;
1058 //--------------------------------------------------------------------
1061 //--------------------------------------------------------------------
1062 /* Consider an input object, for each slice, find the extrema
1063 position according to a given direction and build a line segment
1064 passing throught this point in a given direction. Output is a
1065 vector of line (from point A to B), for each slice;
1067 template<class ImageType>
1069 SliceBySliceBuildLineSegmentAccordingToExtremaPosition(const ImageType * input,
1070 typename ImageType::PixelType BG,
1072 int extremaDirection,
1073 bool extremaOppositeFlag,
1076 std::vector<typename ImageType::PointType> & A,
1077 std::vector<typename ImageType::PointType> & B)
1080 typedef typename itk::Image<typename ImageType::PixelType, ImageType::ImageDimension-1> SliceType;
1082 // Build the list of slices
1083 std::vector<typename SliceType::Pointer> slices;
1084 clitk::ExtractSlices<ImageType>(input, sliceDimension, slices);
1086 // Build the list of 2D points
1087 std::map<int, typename SliceType::PointType> position2D;
1088 for(uint i=0; i<slices.size(); i++) {
1089 typename SliceType::PointType p;
1091 clitk::FindExtremaPointInAGivenDirection<SliceType>(slices[i], BG,
1092 extremaDirection, extremaOppositeFlag, p);
1098 // Convert 2D points in slice into 3D points
1099 clitk::PointsUtils<ImageType>::Convert2DMapTo3DList(position2D, input, A);
1101 // Create additional point just right to the previous ones, on the
1102 // given lineDirection, in order to create a horizontal/vertical line.
1103 for(uint i=0; i<A.size(); i++) {
1104 typename ImageType::PointType p = A[i];
1105 p[lineDirection] += 10;
1108 A[i][extremaDirection] += margin;
1109 B[i][extremaDirection] += margin;
1113 //--------------------------------------------------------------------
1116 //--------------------------------------------------------------------
1117 template<class ImageType>
1118 typename ImageType::Pointer
1119 SliceBySliceKeepMainCCL(const ImageType * input,
1120 typename ImageType::PixelType BG,
1121 typename ImageType::PixelType FG) {
1124 const int d = ImageType::ImageDimension-1;
1125 typedef typename itk::Image<typename ImageType::PixelType, d> SliceType;
1126 std::vector<typename SliceType::Pointer> slices;
1127 clitk::ExtractSlices<ImageType>(input, d, slices);
1129 // Labelize and keep the main one
1130 std::vector<typename SliceType::Pointer> o;
1131 for(uint i=0; i<slices.size(); i++) {
1132 o.push_back(clitk::Labelize<SliceType>(slices[i], BG, false, 1));
1133 o[i] = clitk::KeepLabels<SliceType>(o[i], BG, FG, 1, 1, true);
1137 typename ImageType::Pointer output;
1138 output = clitk::JoinSlices<ImageType>(o, input, d);
1141 //--------------------------------------------------------------------
1144 //--------------------------------------------------------------------
1145 template<class ImageType>
1146 typename ImageType::Pointer
1147 Clone(const ImageType * input) {
1148 typedef itk::ImageDuplicator<ImageType> DuplicatorType;
1149 typename DuplicatorType::Pointer duplicator = DuplicatorType::New();
1150 duplicator->SetInputImage(input);
1151 duplicator->Update();
1152 return duplicator->GetOutput();
1154 //--------------------------------------------------------------------
1157 //--------------------------------------------------------------------
1158 /* Consider an input object, start at A, for each slice (dim1):
1159 - compute the intersection between the AB line and the current slice
1160 - remove what is at lower or greater according to dim2 of this point
1163 template<class ImageType>
1164 typename ImageType::Pointer
1165 SliceBySliceSetBackgroundFromSingleLine(const ImageType * input,
1166 typename ImageType::PixelType BG,
1167 typename ImageType::PointType & A,
1168 typename ImageType::PointType & B,
1169 int dim1, int dim2, bool removeLowerPartFlag)
1173 typedef typename itk::Image<typename ImageType::PixelType, ImageType::ImageDimension-1> SliceType;
1174 typedef typename SliceType::Pointer SlicePointer;
1175 std::vector<SlicePointer> slices;
1176 clitk::ExtractSlices<ImageType>(input, dim1, slices);
1178 // Start at slice that contains A, and stop at B
1179 typename ImageType::IndexType Ap;
1180 typename ImageType::IndexType Bp;
1181 input->TransformPhysicalPointToIndex(A, Ap);
1182 input->TransformPhysicalPointToIndex(B, Bp);
1184 // Determine slice largest region
1185 typename SliceType::RegionType region = slices[0]->GetLargestPossibleRegion();
1186 typename SliceType::SizeType size = region.GetSize();
1187 typename SliceType::IndexType index = region.GetIndex();
1190 double a = (Bp[dim2]-Ap[dim2])/(Bp[dim1]-Ap[dim1]);
1191 double b = Ap[dim2];
1193 // Loop from slice A to slice B
1194 for(uint i=0; i<(Bp[dim1]-Ap[dim1]); i++) {
1195 // Compute intersection between line AB and current slice for the dim2
1197 // Change region (lower than dim2)
1198 if (removeLowerPartFlag) {
1199 size[dim2] = p-Ap[dim2];
1202 size[dim2] = slices[0]->GetLargestPossibleRegion().GetSize()[dim2]-p;
1205 region.SetSize(size);
1206 region.SetIndex(index);
1207 // Fill region with BG (simple region iterator)
1208 FillRegionWithValue<SliceType>(slices[i+Ap[dim1]], BG, region);
1210 typedef itk::ImageRegionIterator<SliceType> IteratorType;
1211 IteratorType iter(slices[i+Ap[dim1]], region);
1213 while (!iter.IsAtEnd()) {
1222 typename ImageType::Pointer output;
1223 output = clitk::JoinSlices<ImageType>(slices, input, dim1);
1226 //--------------------------------------------------------------------
1228 //--------------------------------------------------------------------
1229 /* Consider an input object, slice by slice, use the point A and set
1230 pixel to BG according to their position relatively to A
1232 template<class ImageType>
1233 typename ImageType::Pointer
1234 SliceBySliceSetBackgroundFromPoints(const ImageType * input,
1235 typename ImageType::PixelType BG,
1237 std::vector<typename ImageType::PointType> & A,
1238 bool removeGreaterThanXFlag,
1239 bool removeGreaterThanYFlag)
1243 typedef typename itk::Image<typename ImageType::PixelType, ImageType::ImageDimension-1> SliceType;
1244 typedef typename SliceType::Pointer SlicePointer;
1245 std::vector<SlicePointer> slices;
1246 clitk::ExtractSlices<ImageType>(input, sliceDim, slices);
1248 // Start at slice that contains A
1249 typename ImageType::IndexType Ap;
1251 // Determine slice largest region
1252 typename SliceType::RegionType region = slices[0]->GetLargestPossibleRegion();
1253 typename SliceType::SizeType size = region.GetSize();
1254 typename SliceType::IndexType index = region.GetIndex();
1256 // Loop from slice A to slice B
1257 for(uint i=0; i<A.size(); i++) {
1258 input->TransformPhysicalPointToIndex(A[i], Ap);
1259 uint sliceIndex = Ap[2] - input->GetLargestPossibleRegion().GetIndex()[2];
1260 if ((sliceIndex < 0) || (sliceIndex >= slices.size())) {
1261 continue; // do not consider this slice
1264 // Compute region for BG
1265 if (removeGreaterThanXFlag) {
1267 size[0] = region.GetSize()[0]-(index[0]-region.GetIndex()[0]);
1270 index[0] = region.GetIndex()[0];
1271 size[0] = Ap[0] - index[0];
1274 if (removeGreaterThanYFlag) {
1276 size[1] = region.GetSize()[1]-(index[1]-region.GetIndex()[1]);
1279 index[1] = region.GetIndex()[1];
1280 size[1] = Ap[1] - index[1];
1284 region.SetSize(size);
1285 region.SetIndex(index);
1287 // Fill region with BG (simple region iterator)
1288 FillRegionWithValue<SliceType>(slices[sliceIndex], BG, region);
1293 typename ImageType::Pointer output;
1294 output = clitk::JoinSlices<ImageType>(slices, input, sliceDim);
1297 //--------------------------------------------------------------------
1300 //--------------------------------------------------------------------
1301 template<class ImageType>
1303 FillRegionWithValue(ImageType * input, typename ImageType::PixelType value, typename ImageType::RegionType & region)
1305 typedef itk::ImageRegionIterator<ImageType> IteratorType;
1306 IteratorType iter(input, region);
1308 while (!iter.IsAtEnd()) {
1313 //--------------------------------------------------------------------
1316 //--------------------------------------------------------------------
1317 template<class ImageType>
1319 GetMinMaxBoundary(ImageType * input, typename ImageType::PointType & min,
1320 typename ImageType::PointType & max)
1322 typedef typename ImageType::PointType PointType;
1323 typedef typename ImageType::IndexType IndexType;
1324 IndexType min_i, max_i;
1325 min_i = input->GetLargestPossibleRegion().GetIndex();
1326 for(uint i=0; i<ImageType::ImageDimension; i++)
1327 max_i[i] = input->GetLargestPossibleRegion().GetSize()[i] + min_i[i];
1328 input->TransformIndexToPhysicalPoint(min_i, min);
1329 input->TransformIndexToPhysicalPoint(max_i, max);
1331 //--------------------------------------------------------------------
1334 //--------------------------------------------------------------------
1335 template<class ImageType>
1336 typename itk::Image<float, ImageType::ImageDimension>::Pointer
1337 DistanceMap(const ImageType * input, typename ImageType::PixelType BG)//,
1338 // typename itk::Image<float, ImageType::ImageDimension>::Pointer dmap)
1340 typedef itk::Image<float,ImageType::ImageDimension> FloatImageType;
1341 typedef itk::SignedMaurerDistanceMapImageFilter<ImageType, FloatImageType> DistanceMapFilterType;
1342 typename DistanceMapFilterType::Pointer filter = DistanceMapFilterType::New();
1343 filter->SetInput(input);
1344 filter->SetUseImageSpacing(true);
1345 filter->SquaredDistanceOff();
1346 filter->SetBackgroundValue(BG);
1348 return filter->GetOutput();
1350 //--------------------------------------------------------------------
1353 //--------------------------------------------------------------------
1354 template<class ImageType>
1356 SliceBySliceBuildLineSegmentAccordingToMinimalDistanceBetweenStructures(const ImageType * S1,
1357 const ImageType * S2,
1358 typename ImageType::PixelType BG,
1360 std::vector<typename ImageType::PointType> & A,
1361 std::vector<typename ImageType::PointType> & B)
1364 typedef typename itk::Image<typename ImageType::PixelType, 2> SliceType;
1365 typedef typename SliceType::Pointer SlicePointer;
1366 std::vector<SlicePointer> slices_s1;
1367 std::vector<SlicePointer> slices_s2;
1368 clitk::ExtractSlices<ImageType>(S1, sliceDimension, slices_s1);
1369 clitk::ExtractSlices<ImageType>(S2, sliceDimension, slices_s2);
1371 assert(slices_s1.size() == slices_s2.size());
1374 typedef itk::Image<float,2> FloatImageType;
1375 typedef itk::SignedMaurerDistanceMapImageFilter<SliceType, FloatImageType> DistanceMapFilterType;
1376 std::vector<typename FloatImageType::Pointer> dmaps1;
1377 std::vector<typename FloatImageType::Pointer> dmaps2;
1378 typename FloatImageType::Pointer dmap;
1381 for(uint i=0; i<slices_s1.size(); i++) {
1382 // Compute dmap for S1 *TO PUT IN FONCTION*
1383 dmap = clitk::DistanceMap<SliceType>(slices_s1[i], BG);
1384 dmaps1.push_back(dmap);
1385 writeImage<FloatImageType>(dmap, "dmap1.mha");
1386 // Compute dmap for S2
1387 dmap = clitk::DistanceMap<SliceType>(slices_s2[i], BG);
1388 dmaps2.push_back(dmap);
1389 writeImage<FloatImageType>(dmap, "dmap2.mha");
1391 // Look in S2 for the point the closest to S1
1392 typename SliceType::PointType p = ComputeClosestPoint<SliceType>(slices_s1[i], dmaps2[i], BG);
1393 typename ImageType::PointType p3D;
1394 clitk::PointsUtils<ImageType>::Convert2DTo3D(p, S1, i, p3D);
1397 // Look in S2 for the point the closest to S1
1398 p = ComputeClosestPoint<SliceType>(slices_s2[i], dmaps1[i], BG);
1399 clitk::PointsUtils<ImageType>::Convert2DTo3D(p, S2, i, p3D);
1406 typedef itk::Image<float,3> FT;
1407 FT::Pointer f = FT::New();
1408 typename FT::Pointer d1 = clitk::JoinSlices<FT>(dmaps1, S1, 2);
1409 typename FT::Pointer d2 = clitk::JoinSlices<FT>(dmaps2, S2, 2);
1410 writeImage<FT>(d1, "d1.mha");
1411 writeImage<FT>(d2, "d2.mha");
1414 //--------------------------------------------------------------------
1417 //--------------------------------------------------------------------
1418 template<class ImageType>
1419 typename ImageType::PointType
1420 ComputeClosestPoint(const ImageType * input,
1421 const itk::Image<float, ImageType::ImageDimension> * dmap,
1422 typename ImageType::PixelType & BG)
1424 // Loop dmap + S2, if FG, get min
1425 typedef itk::Image<float,ImageType::ImageDimension> FloatImageType;
1426 typedef itk::ImageRegionConstIteratorWithIndex<ImageType> ImageIteratorType;
1427 typedef itk::ImageRegionConstIterator<FloatImageType> DMapIteratorType;
1428 ImageIteratorType iter1(input, input->GetLargestPossibleRegion());
1429 DMapIteratorType iter2(dmap, dmap->GetLargestPossibleRegion());
1433 double dmin = 100000.0;
1434 typename ImageType::IndexType indexmin;
1435 while (!iter1.IsAtEnd()) {
1436 if (iter1.Get() != BG) {
1437 double d = iter2.Get();
1439 indexmin = iter1.GetIndex();
1448 typename ImageType::PointType p;
1449 input->TransformIndexToPhysicalPoint(indexmin, p);
1452 //--------------------------------------------------------------------
1457 } // end of namespace