]> Creatis software - clitk.git/blob - itk/clitkSegmentationUtils.txx
adca988098ad8c78fee5cfb132b890fb0bef0c09
[clitk.git] / itk / clitkSegmentationUtils.txx
1 /*=========================================================================
2   Program:   vv                     http://www.creatis.insa-lyon.fr/rio/vv
3
4   Authors belong to: 
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
8
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.
12
13   It is distributed under dual licence
14
15   - BSD        See included LICENSE.txt file
16   - CeCILL-B   http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
17   ======================================================================-====*/
18
19 // clitk
20 #include "clitkSetBackgroundImageFilter.h"
21 #include "clitkSliceBySliceRelativePositionFilter.h"
22 #include "clitkCropLikeImageFilter.h"
23 #include "clitkMemoryUsage.h"
24
25 // itk
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>
38
39 namespace clitk {
40
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;
48     IndexType firstIndex;
49     IndexType lastIndex;
50     for(unsigned int i=0; i<image->GetImageDimension(); i++) {
51       firstIndex[i] = region.GetIndex()[i];
52       lastIndex[i] = firstIndex[i]+region.GetSize()[i];
53     }
54
55     typedef itk::BoundingBox<unsigned long, 
56                              ImageType::ImageDimension> BBType;
57     typedef typename BBType::PointType PointType;
58     PointType lastPoint;
59     PointType firstPoint;
60     image->TransformIndexToPhysicalPoint(firstIndex, firstPoint);
61     image->TransformIndexToPhysicalPoint(lastIndex, lastPoint);
62
63     bb->SetMaximum(lastPoint);
64     bb->SetMinimum(firstPoint);
65   }
66   //--------------------------------------------------------------------
67
68
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) {
74
75     typedef itk::BoundingBox<unsigned long, Dimension> BBType;
76     typedef typename BBType::PointType PointType;
77     PointType lastPoint;
78     PointType firstPoint;
79
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]);
85     }
86
87     bbo->SetMaximum(lastPoint);
88     bbo->SetMinimum(firstPoint);
89   }
90   //--------------------------------------------------------------------
91
92
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) {
99     // Types
100     typedef typename ImageType::IndexType  IndexType;
101     typedef typename ImageType::PointType  PointType;
102     typedef typename ImageType::RegionType RegionType;
103     typedef typename ImageType::SizeType   SizeType;
104
105     // Region starting point
106     IndexType regionStart;
107     PointType start = bb->GetMinimum();
108     image->TransformPhysicalPointToIndex(start, regionStart);
109     
110     // Region size
111     SizeType regionSize;
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]);
116     }
117    
118     // Create region
119     region.SetIndex(regionStart);
120     region.SetSize(regionSize);
121   }
122   //--------------------------------------------------------------------
123
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, 
131                 bool inPlace) {
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();
144   }
145   //--------------------------------------------------------------------
146
147
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();
160   
161     // Return result
162     return connectFilter->GetObjectCount();
163   }
164   //--------------------------------------------------------------------
165
166   //--------------------------------------------------------------------
167   /*
168     Warning : in this cas, we consider outputType like inputType, not
169     InternalImageType. Be sure it fits.
170   */
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;
179   
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);
187   
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();
195
196     // Return result
197     typename ImageType::Pointer output = relabelFilter->GetOutput();
198     return output;
199   }
200   //--------------------------------------------------------------------
201
202
203   //--------------------------------------------------------------------
204   /*
205     Warning : in this cas, we consider outputType like inputType, not
206     InternalImageType. Be sure it fits.
207   */
208   template<class ImageType>
209   typename ImageType::Pointer
210   LabelizeAndCountNumberOfObjects(const ImageType * input, 
211                                   typename ImageType::PixelType BG, 
212                                   bool isFullyConnected, 
213                                   int minimalComponentSize, 
214                                   int & nb) {
215     // InternalImageType for storing large number of component
216     typedef itk::Image<int, ImageType::ImageDimension> InternalImageType;
217   
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);
225   
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();
233
234     nb = relabelFilter->GetNumberOfObjects();
235     // DD(relabelFilter->GetOriginalNumberOfObjects());
236     // DD(relabelFilter->GetSizeOfObjectsInPhysicalUnits()[0]);
237
238     // Return result
239     typename ImageType::Pointer output = relabelFilter->GetOutput();
240     return output;
241   }
242   //--------------------------------------------------------------------
243
244
245
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();
263     }
264     return working_image;
265   }
266   //--------------------------------------------------------------------
267
268
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, 
277              bool useLastKeep) {
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();
287   }
288   //--------------------------------------------------------------------
289
290
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)
300   {
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, 
306                                           BG, FG, 
307                                           param->GetFirstKeep(), 
308                                           param->GetLastKeep(), 
309                                           param->GetUseLastKeep());
310     return working_image;
311   }
312   //--------------------------------------------------------------------
313
314
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) 
321   {
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();  
329   }
330   //--------------------------------------------------------------------
331
332
333   //--------------------------------------------------------------------
334   template<class MaskImageType>
335   typename MaskImageType::Pointer
336   SliceBySliceRelativePosition(const MaskImageType * input,
337                                const MaskImageType * object,
338                                int direction, 
339                                double threshold, 
340                                std::string orientation, 
341                                bool uniqueConnectedComponent, 
342                                double spacing, 
343                                bool autocropFlag, 
344                                bool singleObjectCCL) 
345   {
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();
365   }
366   //--------------------------------------------------------------------
367
368
369   //--------------------------------------------------------------------
370   template<class ImageType>
371   bool
372   FindExtremaPointInAGivenDirection(const ImageType * input, 
373                                     typename ImageType::PixelType bg, 
374                                     int direction, bool opposite, 
375                                     typename ImageType::PointType & point)
376   {
377     typename ImageType::PointType dummy;
378     return FindExtremaPointInAGivenDirection(input, bg, direction, 
379                                              opposite, dummy, 0, point);
380   }
381   //--------------------------------------------------------------------
382
383
384   //--------------------------------------------------------------------
385   template<class ImageType>
386   bool
387   FindExtremaPointInAGivenDirection(const ImageType * input, 
388                                     typename ImageType::PixelType bg, 
389                                     int direction, bool opposite, 
390                                     typename ImageType::PointType refpoint,
391                                     double distanceMax, 
392                                     typename ImageType::PointType & point)
393   {
394     /*
395       loop over input pixels, store the index in the fg that is max
396       according to the given direction. 
397     */    
398     typedef itk::ImageRegionConstIteratorWithIndex<ImageType> IteratorType;
399     IteratorType iter(input, input->GetLargestPossibleRegion());
400     iter.GoToBegin();
401     typename ImageType::IndexType max = input->GetLargestPossibleRegion().GetIndex();
402     if (opposite) max = max+input->GetLargestPossibleRegion().GetSize();
403     bool found=false;
404     while (!iter.IsAtEnd()) {
405       if (iter.Get() != bg) {
406         bool test = iter.GetIndex()[direction] >  max[direction];
407         if (opposite) test = !test;
408         if (test) {
409           typename ImageType::PointType p;
410           input->TransformIndexToPhysicalPoint(iter.GetIndex(), p);
411           if ((distanceMax==0) || (p.EuclideanDistanceTo(refpoint) < distanceMax)) {
412             max = iter.GetIndex();
413             found = true;
414           }
415         }
416       }
417       ++iter;
418     }
419     if (!found) return false;
420     input->TransformIndexToPhysicalPoint(max, point);
421     return true;
422   }
423   //--------------------------------------------------------------------
424
425
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) 
432   {
433     return CropImageAlongOneAxis<ImageType>(image, dim, 
434                                             image->GetOrigin()[dim], 
435                                             min,
436                                             autoCrop, BG);
437   }
438   //--------------------------------------------------------------------
439
440
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) 
447   {
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);
452   }
453   //--------------------------------------------------------------------
454
455
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) 
462   {
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     p[dim] = min;
468     typename ImageType::IndexType start;
469     image->TransformPhysicalPointToIndex(p, start);
470     p[dim] = max;
471     typename ImageType::IndexType end;
472     image->TransformPhysicalPointToIndex(p, end);
473     size[dim] = abs(end[dim]-start[dim]);
474     region.SetIndex(start);
475     region.SetSize(size);
476   
477     // Perform Crop
478     typedef itk::RegionOfInterestImageFilter<ImageType, ImageType> CropFilterType;
479     typename CropFilterType::Pointer cropFilter = CropFilterType::New();
480     cropFilter->SetInput(image);
481     cropFilter->SetRegionOfInterest(region);
482     cropFilter->Update();
483     typename ImageType::Pointer result = cropFilter->GetOutput();
484   
485     // Auto Crop
486     if (autoCrop) {
487       result = AutoCrop<ImageType>(result, BG);
488     }
489     return result;
490   }
491   //--------------------------------------------------------------------
492
493
494   //--------------------------------------------------------------------
495   template<class ImageType>
496   void
497   ComputeCentroids(const ImageType * image, 
498                    typename ImageType::PixelType BG, 
499                    std::vector<typename ImageType::PointType> & centroids) 
500   {
501     typedef long LabelType;
502     static const unsigned int Dim = ImageType::ImageDimension;
503     typedef itk::ShapeLabelObject< LabelType, Dim > LabelObjectType;
504     typedef itk::LabelMap< LabelObjectType > LabelMapType;
505     typedef itk::LabelImageToLabelMapFilter<ImageType, LabelMapType> ImageToMapFilterType;
506     typename ImageToMapFilterType::Pointer imageToLabelFilter = ImageToMapFilterType::New(); 
507     typedef itk::ShapeLabelMapFilter<LabelMapType, ImageType> ShapeFilterType; 
508     typename ShapeFilterType::Pointer statFilter = ShapeFilterType::New();
509     imageToLabelFilter->SetBackgroundValue(BG);
510     imageToLabelFilter->SetInput(image);
511     statFilter->SetInput(imageToLabelFilter->GetOutput());
512     statFilter->Update();
513     typename LabelMapType::Pointer labelMap = statFilter->GetOutput();
514
515     centroids.clear();
516     typename ImageType::PointType dummy;
517     centroids.push_back(dummy); // label 0 -> no centroid, use dummy point for BG 
518     //DS FIXME (not useful ! to change ..)
519     for(uint i=0; i<labelMap->GetNumberOfLabelObjects(); i++) {
520       int label = labelMap->GetLabels()[i];
521       centroids.push_back(labelMap->GetLabelObject(label)->GetCentroid());
522     } 
523   }
524   //--------------------------------------------------------------------
525
526
527   //--------------------------------------------------------------------
528   template<class ImageType, class LabelType>
529   typename itk::LabelMap< itk::ShapeLabelObject<LabelType, ImageType::ImageDimension> >::Pointer
530   ComputeLabelMap(const ImageType * image, 
531                   typename ImageType::PixelType BG, 
532                   bool computePerimeterFlag) 
533   {
534     static const unsigned int Dim = ImageType::ImageDimension;
535     typedef itk::ShapeLabelObject< LabelType, Dim > LabelObjectType;
536     typedef itk::LabelMap< LabelObjectType > LabelMapType;
537     typedef itk::LabelImageToLabelMapFilter<ImageType, LabelMapType> ImageToMapFilterType;
538     typename ImageToMapFilterType::Pointer imageToLabelFilter = ImageToMapFilterType::New(); 
539     typedef itk::ShapeLabelMapFilter<LabelMapType, ImageType> ShapeFilterType; 
540     typename ShapeFilterType::Pointer statFilter = ShapeFilterType::New();
541     imageToLabelFilter->SetBackgroundValue(BG);
542     imageToLabelFilter->SetInput(image);
543     statFilter->SetInput(imageToLabelFilter->GetOutput());
544     statFilter->SetComputePerimeter(computePerimeterFlag);
545     statFilter->Update();
546     return statFilter->GetOutput();
547   }
548   //--------------------------------------------------------------------
549
550
551   //--------------------------------------------------------------------
552   template<class ImageType>
553   void
554   ComputeCentroids2(const ImageType * image, 
555                    typename ImageType::PixelType BG, 
556                    std::vector<typename ImageType::PointType> & centroids) 
557   {
558     typedef long LabelType;
559     static const unsigned int Dim = ImageType::ImageDimension;
560     typedef itk::ShapeLabelObject< LabelType, Dim > LabelObjectType;
561     typedef itk::LabelMap< LabelObjectType > LabelMapType;
562     typedef itk::LabelImageToLabelMapFilter<ImageType, LabelMapType> ImageToMapFilterType;
563     typename ImageToMapFilterType::Pointer imageToLabelFilter = ImageToMapFilterType::New(); 
564     typedef itk::ShapeLabelMapFilter<LabelMapType, ImageType> ShapeFilterType; 
565     typename ShapeFilterType::Pointer statFilter = ShapeFilterType::New();
566     imageToLabelFilter->SetBackgroundValue(BG);
567     imageToLabelFilter->SetInput(image);
568     statFilter->SetInput(imageToLabelFilter->GetOutput());
569     statFilter->Update();
570     typename LabelMapType::Pointer labelMap = statFilter->GetOutput();
571
572     centroids.clear();
573     typename ImageType::PointType dummy;
574     centroids.push_back(dummy); // label 0 -> no centroid, use dummy point
575     for(uint i=1; i<labelMap->GetNumberOfLabelObjects()+1; i++) {
576       centroids.push_back(labelMap->GetLabelObject(i)->GetCentroid());
577     } 
578     
579     for(uint i=1; i<labelMap->GetNumberOfLabelObjects()+1; i++) {
580       DD(labelMap->GetLabelObject(i)->GetBinaryPrincipalAxes());
581       DD(labelMap->GetLabelObject(i)->GetBinaryFlatness());
582       DD(labelMap->GetLabelObject(i)->GetRoundness ());      
583
584       // search for the point on the boundary alog PA
585
586     }
587
588   }
589   //--------------------------------------------------------------------
590
591
592   //--------------------------------------------------------------------
593   template<class ImageType>
594   void
595   PointsUtils<ImageType>::Convert2DTo3D(const PointType2D & p2D, 
596                                         const ImageType * image, 
597                                         const int slice, 
598                                         PointType3D & p3D)  
599   {
600     IndexType3D index3D;
601     index3D[0] = index3D[1] = 0;
602     index3D[2] = image->GetLargestPossibleRegion().GetIndex()[2]+slice;
603     image->TransformIndexToPhysicalPoint(index3D, p3D);
604     p3D[0] = p2D[0]; 
605     p3D[1] = p2D[1];
606     //  p3D[2] = p[2];//(image->GetLargestPossibleRegion().GetIndex()[2]+slice)*image->GetSpacing()[2] 
607     //    + image->GetOrigin()[2];
608   }
609   //--------------------------------------------------------------------
610
611
612   //--------------------------------------------------------------------
613   template<class ImageType>
614   void 
615   PointsUtils<ImageType>::Convert2DMapTo3DList(const MapPoint2DType & map, 
616                                             const ImageType * image, 
617                                             VectorPoint3DType & list)
618   {
619     typename MapPoint2DType::const_iterator iter = map.begin();
620     while (iter != map.end()) {
621       PointType3D p;
622       Convert2DTo3D(iter->second, image, iter->first, p);
623       list.push_back(p);
624       ++iter;
625     }
626   }
627   //--------------------------------------------------------------------
628
629
630   //--------------------------------------------------------------------
631   template<class ImageType>
632   void 
633   PointsUtils<ImageType>::Convert2DListTo3DList(const VectorPoint2DType & p2D, 
634                                                 int slice,
635                                                 const ImageType * image, 
636                                                 VectorPoint3DType & list) 
637   {
638     for(uint i=0; i<p2D.size(); i++) {
639       PointType3D p;
640       Convert2DTo3D(p2D[i], image, slice, p);
641       list.push_back(p);
642     }
643   }
644   //--------------------------------------------------------------------
645
646
647   //--------------------------------------------------------------------
648   template<class ImageType>
649   void 
650   WriteListOfLandmarks(std::vector<typename ImageType::PointType> points, 
651                        std::string filename)
652   {
653     std::ofstream os; 
654     openFileForWriting(os, filename); 
655     os << "LANDMARKS1" << std::endl;  
656     for(uint i=0; i<points.size(); i++) {
657       const typename ImageType::PointType & p = points[i];
658       // Write it in the file
659       os << i << " " << p[0] << " " << p[1] << " " << p[2] << " 0 0 " << std::endl;
660     }
661     os.close();
662   }
663   //--------------------------------------------------------------------
664
665
666   //--------------------------------------------------------------------
667   template<class ImageType>
668   typename ImageType::Pointer 
669   Dilate(const ImageType * image, double radiusInMM,               
670          typename ImageType::PixelType BG,
671          typename ImageType::PixelType FG,  
672          bool extendSupport)
673   {
674     typename ImageType::SizeType r;
675     for(uint i=0; i<ImageType::ImageDimension; i++) 
676       r[i] = (uint)lrint(radiusInMM/image->GetSpacing()[i]);
677     return Dilate<ImageType>(image, r, BG, FG, extendSupport);
678   }
679   //--------------------------------------------------------------------
680
681
682   //--------------------------------------------------------------------
683   template<class ImageType>
684   typename ImageType::Pointer 
685   Dilate(const ImageType * image, typename ImageType::PointType radiusInMM, 
686          typename ImageType::PixelType BG, 
687          typename ImageType::PixelType FG, 
688          bool extendSupport)
689   {
690     typename ImageType::SizeType r;
691     for(uint i=0; i<ImageType::ImageDimension; i++) 
692       r[i] = (uint)lrint(radiusInMM[i]/image->GetSpacing()[i]);
693     return Dilate<ImageType>(image, r, BG, FG, extendSupport);
694   }
695   //--------------------------------------------------------------------
696
697
698   //--------------------------------------------------------------------
699   template<class ImageType>
700   typename ImageType::Pointer 
701   Dilate(const ImageType * image, typename ImageType::SizeType radius, 
702          typename ImageType::PixelType BG, 
703          typename ImageType::PixelType FG, 
704          bool extendSupport)
705   {
706     // Create kernel for dilatation
707     typedef itk::BinaryBallStructuringElement<typename ImageType::PixelType, 
708                                               ImageType::ImageDimension> KernelType;
709     KernelType structuringElement;
710     structuringElement.SetRadius(radius);
711     structuringElement.CreateStructuringElement();
712
713     typename ImageType::Pointer output;
714     if (extendSupport) {
715       typedef itk::ConstantPadImageFilter<ImageType, ImageType> PadFilterType;
716       typename PadFilterType::Pointer padFilter = PadFilterType::New();
717       padFilter->SetInput(image);
718       typename ImageType::SizeType lower;
719       typename ImageType::SizeType upper;
720       for(uint i=0; i<3; i++) {
721         lower[i] = upper[i] = 2*(radius[i]+1);
722       }
723       padFilter->SetPadLowerBound(lower);
724       padFilter->SetPadUpperBound(upper);
725       padFilter->Update();
726       output = padFilter->GetOutput();
727     }
728
729     // Dilate  filter
730     typedef itk::BinaryDilateImageFilter<ImageType, ImageType , KernelType> DilateFilterType;
731     typename DilateFilterType::Pointer dilateFilter = DilateFilterType::New();
732     dilateFilter->SetBackgroundValue(BG);
733     dilateFilter->SetForegroundValue(FG);
734     dilateFilter->SetBoundaryToForeground(false);
735     dilateFilter->SetKernel(structuringElement);
736     if (extendSupport) dilateFilter->SetInput(output);
737     else dilateFilter->SetInput(image);
738     dilateFilter->Update();
739     return dilateFilter->GetOutput();
740   }
741   //--------------------------------------------------------------------
742
743
744   //--------------------------------------------------------------------
745   template<class ImageType>
746   typename ImageType::Pointer 
747   Opening(const ImageType * image, typename ImageType::SizeType radius,
748          typename ImageType::PixelType BG,
749          typename ImageType::PixelType FG)
750   {
751     // Kernel 
752     typedef itk::BinaryBallStructuringElement<typename ImageType::PixelType, 
753                                               ImageType::ImageDimension> KernelType;    
754     KernelType structuringElement;
755     structuringElement.SetRadius(radius);
756     structuringElement.CreateStructuringElement();
757     
758     // Filter
759     typedef itk::BinaryMorphologicalOpeningImageFilter<ImageType, ImageType , KernelType> OpeningFilterType;
760     typename OpeningFilterType::Pointer open = OpeningFilterType::New();
761     open->SetInput(image);
762     open->SetBackgroundValue(BG);
763     open->SetForegroundValue(FG);
764     open->SetKernel(structuringElement);
765     open->Update();
766     return open->GetOutput();
767   }
768   //--------------------------------------------------------------------
769
770
771
772   //--------------------------------------------------------------------
773   template<class ValueType, class VectorType>
774   void ConvertOption(std::string optionName, uint given, 
775                      ValueType * values, VectorType & p, 
776                      uint dim, bool required) 
777   {
778     if (required && (given == 0)) {
779       clitkExceptionMacro("The option --" << optionName << " must be set and have 1 or " 
780                           << dim << " values.");
781     }
782     if (given == 1) {
783       for(uint i=0; i<dim; i++) p[i] = values[0];
784       return;
785     }
786     if (given == dim) {
787       for(uint i=0; i<dim; i++) p[i] = values[i];
788       return;
789     }
790     if (given == 0) return;
791     clitkExceptionMacro("The option --" << optionName << " must have 1 or " 
792                         << dim << " values.");
793   }
794   //--------------------------------------------------------------------
795
796
797   //--------------------------------------------------------------------
798   /*
799     http://www.gamedev.net/community/forums/topic.asp?topic_id=542870
800     Assuming the points are (Ax,Ay) (Bx,By) and (Cx,Cy), you need to compute:
801     (Bx - Ax) * (Cy - Ay) - (By - Ay) * (Cx - Ax)
802     This will equal zero if the point C is on the line formed by
803     points A and B, and will have a different sign depending on the
804     side. Which side this is depends on the orientation of your (x,y)
805     coordinates, but you can plug test values for A,B and C into this
806     formula to determine whether negative values are to the left or to
807     the right.
808     => to accelerate, start with formula, when change sign -> stop and fill
809
810     offsetToKeep = is used to determine which side of the line we
811     keep. The point along the mainDirection but 'offsetToKeep' mm away
812     is kept.
813   
814   */
815   template<class ImageType>
816   void 
817   SliceBySliceSetBackgroundFromLineSeparation(ImageType * input, 
818                                               std::vector<typename ImageType::PointType> & lA, 
819                                               std::vector<typename ImageType::PointType> & lB, 
820                                               typename ImageType::PixelType BG, 
821                                               int mainDirection, 
822                                               double offsetToKeep)
823   {
824     assert((mainDirection==0) || (mainDirection==1));
825     typedef itk::ImageSliceIteratorWithIndex<ImageType> SliceIteratorType;
826     SliceIteratorType siter = SliceIteratorType(input, 
827                                                 input->GetLargestPossibleRegion());
828     siter.SetFirstDirection(0);
829     siter.SetSecondDirection(1);
830     siter.GoToBegin();
831     uint i=0;
832     typename ImageType::PointType A;
833     typename ImageType::PointType B;
834     typename ImageType::PointType C;
835     assert(lA.size() == lB.size());
836     while ((i<lA.size()) && (!siter.IsAtEnd())) {
837       // Check that the current slice correspond to the current point
838       input->TransformIndexToPhysicalPoint(siter.GetIndex(), C);
839       if ((fabs(C[2] - lA[i][2]))>0.01) { // is !equal with a tolerance of 0.01 mm
840       }
841       else {
842         // Define A,B,C points
843         A = lA[i];
844         B = lB[i];
845         C = A;
846       
847         // Check that the line is not a point (A=B)
848         bool p = (A[0] == B[0]) && (A[1] == B[1]);
849       
850         if (!p) {
851           C[mainDirection] += offsetToKeep; // I know I must keep this point
852           double s = (B[0] - A[0]) * (C[1] - A[1]) - (B[1] - A[1]) * (C[0] - A[0]);
853           bool isPositive = s<0;
854           while (!siter.IsAtEndOfSlice()) {
855             while (!siter.IsAtEndOfLine()) {
856               // Very slow, I know ... but image should be very small
857               input->TransformIndexToPhysicalPoint(siter.GetIndex(), C);
858               double s = (B[0] - A[0]) * (C[1] - A[1]) - (B[1] - A[1]) * (C[0] - A[0]);
859               if (s == 0) siter.Set(BG); // on the line, we decide to remove
860               if (isPositive) {
861                 if (s > 0) siter.Set(BG);
862               }
863               else {
864                 if (s < 0) siter.Set(BG); 
865               }
866               ++siter;
867             }
868             siter.NextLine();
869           } // end loop slice
870         }      
871
872         ++i;
873       } // End of current slice
874       siter.NextSlice();
875     }
876   }                                                   
877   //--------------------------------------------------------------------
878
879
880   //--------------------------------------------------------------------
881   template<class ImageType>
882   void 
883   AndNot(ImageType * input, 
884          const ImageType * object, 
885          typename ImageType::PixelType BG)
886   {
887     typename ImageType::Pointer o;
888     bool resized=false;
889     if (!clitk::HaveSameSizeAndSpacing<ImageType, ImageType>(input, object)) {
890       o = clitk::ResizeImageLike<ImageType>(object, input, BG);
891       resized = true;
892     }
893
894     typedef clitk::BooleanOperatorLabelImageFilter<ImageType> BoolFilterType;
895     typename BoolFilterType::Pointer boolFilter = BoolFilterType::New(); 
896     boolFilter->InPlaceOn();
897     boolFilter->SetInput1(input);
898     if (resized) boolFilter->SetInput2(o);  
899     else boolFilter->SetInput2(object);
900     boolFilter->SetBackgroundValue1(BG);
901     boolFilter->SetBackgroundValue2(BG);
902     boolFilter->SetOperationType(BoolFilterType::AndNot);
903     boolFilter->Update();
904   }
905   //--------------------------------------------------------------------
906
907
908   //--------------------------------------------------------------------
909   template<class ImageType>
910   void 
911   And(ImageType * input, 
912       const ImageType * object, 
913       typename ImageType::PixelType BG)
914   {
915     typename ImageType::Pointer o;
916     bool resized=false;
917     if (!clitk::HaveSameSizeAndSpacing<ImageType, ImageType>(input, object)) {
918       o = clitk::ResizeImageLike<ImageType>(object, input, BG);
919       resized = true;
920     }
921
922     typedef clitk::BooleanOperatorLabelImageFilter<ImageType> BoolFilterType;
923     typename BoolFilterType::Pointer boolFilter = BoolFilterType::New(); 
924     boolFilter->InPlaceOn();
925     boolFilter->SetInput1(input);
926     if (resized) boolFilter->SetInput2(o);  
927     else boolFilter->SetInput2(object);
928     boolFilter->SetBackgroundValue1(BG);
929     boolFilter->SetBackgroundValue2(BG);
930     boolFilter->SetOperationType(BoolFilterType::And);
931     boolFilter->Update();
932   }
933   //--------------------------------------------------------------------
934
935
936   //--------------------------------------------------------------------
937   template<class ImageType>
938   void 
939   Or(ImageType * input, 
940      const ImageType * object, 
941      typename ImageType::PixelType BG)
942   {
943     typename ImageType::Pointer o;
944     bool resized=false;
945     if (!clitk::HaveSameSizeAndSpacing<ImageType, ImageType>(input, object)) {
946       o = clitk::ResizeImageLike<ImageType>(object, input, BG);
947       resized = true;
948     }
949
950     typedef clitk::BooleanOperatorLabelImageFilter<ImageType> BoolFilterType;
951     typename BoolFilterType::Pointer boolFilter = BoolFilterType::New(); 
952     boolFilter->InPlaceOn();
953     boolFilter->SetInput1(input);
954     if (resized) boolFilter->SetInput2(o);  
955     else boolFilter->SetInput2(object);
956     boolFilter->SetBackgroundValue1(BG);
957     boolFilter->SetBackgroundValue2(BG);
958     boolFilter->SetOperationType(BoolFilterType::Or);
959     boolFilter->Update();
960   }
961   //--------------------------------------------------------------------
962
963
964   //--------------------------------------------------------------------
965   template<class ImageType>
966   typename ImageType::Pointer
967   Binarize(const ImageType * input, 
968            typename ImageType::PixelType lower, 
969            typename ImageType::PixelType upper, 
970            typename ImageType::PixelType BG,
971            typename ImageType::PixelType FG) 
972   {
973     typedef itk::BinaryThresholdImageFilter<ImageType, ImageType> BinaryThresholdFilterType;
974     typename BinaryThresholdFilterType::Pointer binarizeFilter = BinaryThresholdFilterType::New();
975     binarizeFilter->SetInput(input);
976     binarizeFilter->InPlaceOff();
977     binarizeFilter->SetLowerThreshold(lower);
978     binarizeFilter->SetUpperThreshold(upper);
979     binarizeFilter->SetInsideValue(FG);
980     binarizeFilter->SetOutsideValue(BG);
981     binarizeFilter->Update();
982     return binarizeFilter->GetOutput();
983   }
984   //--------------------------------------------------------------------
985
986
987   //--------------------------------------------------------------------
988   template<class ImageType>
989   void
990   GetMinMaxPointPosition(const ImageType * input, 
991                          typename ImageType::PointType & min,
992                          typename ImageType::PointType & max) 
993   {
994     typename ImageType::IndexType index = input->GetLargestPossibleRegion().GetIndex();
995     input->TransformIndexToPhysicalPoint(index, min);
996     index = index+input->GetLargestPossibleRegion().GetSize();
997     input->TransformIndexToPhysicalPoint(index, max);
998   }
999   //--------------------------------------------------------------------
1000
1001
1002   //--------------------------------------------------------------------
1003   template<class ImageType>
1004   typename ImageType::PointType
1005   FindExtremaPointInAGivenLine(const ImageType * input, 
1006                                int dimension, 
1007                                bool inverse, 
1008                                typename ImageType::PointType p, 
1009                                typename ImageType::PixelType BG, 
1010                                double distanceMax) 
1011   {
1012     // Which direction ?  Increasing or decreasing.
1013     int d=1;
1014     if (inverse) d=-1;
1015   
1016     // Transform to pixel index
1017     typename ImageType::IndexType index;
1018     input->TransformPhysicalPointToIndex(p, index);
1019
1020     // Loop while inside the mask;
1021     while (input->GetPixel(index) != BG) {
1022       index[dimension] += d;
1023     }
1024
1025     // Transform back to Physical Units
1026     typename ImageType::PointType result;
1027     input->TransformIndexToPhysicalPoint(index, result);
1028
1029     // Check that is is not too far away
1030     double distance = p.EuclideanDistanceTo(result);
1031     if (distance > distanceMax) {
1032       result = p; // Get back to initial value
1033     }
1034
1035     return result;
1036   }
1037   //--------------------------------------------------------------------
1038
1039
1040   //--------------------------------------------------------------------
1041   template<class PointType>
1042   bool
1043   IsOnTheSameLineSide(PointType C, PointType A, PointType B, PointType like) 
1044   {
1045     // Look at the position of point 'like' according to the AB line
1046     double s = (B[0] - A[0]) * (like[1] - A[1]) - (B[1] - A[1]) * (like[0] - A[0]);
1047     bool negative = s<0;
1048   
1049     // Look the C position
1050     s = (B[0] - A[0]) * (C[1] - A[1]) - (B[1] - A[1]) * (C[0] - A[0]);
1051
1052     if (negative && (s<=0)) return true;
1053     if (!negative && (s>=0)) return true;
1054     return false;
1055   }
1056   //--------------------------------------------------------------------
1057
1058
1059   //--------------------------------------------------------------------
1060   /* Consider an input object, for each slice, find the extrema
1061      position according to a given direction and build a line segment
1062      passing throught this point in a given direction.  Output is a
1063      vector of line (from point A to B), for each slice;
1064    */
1065   template<class ImageType>
1066   void 
1067   SliceBySliceBuildLineSegmentAccordingToExtremaPosition(const ImageType * input, 
1068                                                          typename ImageType::PixelType BG, 
1069                                                          int sliceDimension, 
1070                                                          int extremaDirection, 
1071                                                          bool extremaOppositeFlag, 
1072                                                          int lineDirection,
1073                                                          double margin,
1074                                                          std::vector<typename ImageType::PointType> & A, 
1075                                                          std::vector<typename ImageType::PointType> & B)
1076   {
1077     // Type of a slice
1078     typedef typename itk::Image<typename ImageType::PixelType, ImageType::ImageDimension-1> SliceType;
1079     
1080     // Build the list of slices
1081     std::vector<typename SliceType::Pointer> slices;
1082     clitk::ExtractSlices<ImageType>(input, sliceDimension, slices);
1083
1084     // Build the list of 2D points
1085     std::map<int, typename SliceType::PointType> position2D;
1086     for(uint i=0; i<slices.size(); i++) {
1087       typename SliceType::PointType p;
1088       bool found = 
1089         clitk::FindExtremaPointInAGivenDirection<SliceType>(slices[i], BG, 
1090                                                             extremaDirection, extremaOppositeFlag, p);
1091       if (found) {
1092         position2D[i] = p;
1093       }
1094     }
1095     
1096     // Convert 2D points in slice into 3D points
1097     clitk::PointsUtils<ImageType>::Convert2DMapTo3DList(position2D, input, A);
1098     
1099     // Create additional point just right to the previous ones, on the
1100     // given lineDirection, in order to create a horizontal/vertical line.
1101     for(uint i=0; i<A.size(); i++) {
1102       typename ImageType::PointType p = A[i];
1103       p[lineDirection] += 10;
1104       B.push_back(p);
1105       // Margins ?
1106       A[i][extremaDirection] += margin;
1107       B[i][extremaDirection] += margin;
1108     }
1109
1110   }
1111   //--------------------------------------------------------------------
1112
1113
1114   //--------------------------------------------------------------------
1115   template<class ImageType>
1116   typename ImageType::Pointer
1117   SliceBySliceKeepMainCCL(const ImageType * input, 
1118                           typename ImageType::PixelType BG,
1119                           typename ImageType::PixelType FG)  {
1120     
1121     // Extract slices
1122     const int d = ImageType::ImageDimension-1;
1123     typedef typename itk::Image<typename ImageType::PixelType, d> SliceType;
1124     std::vector<typename SliceType::Pointer> slices;
1125     clitk::ExtractSlices<ImageType>(input, d, slices);
1126     
1127     // Labelize and keep the main one
1128     std::vector<typename SliceType::Pointer> o;
1129     for(uint i=0; i<slices.size(); i++) {
1130       o.push_back(clitk::Labelize<SliceType>(slices[i], BG, false, 1));
1131       o[i] = clitk::KeepLabels<SliceType>(o[i], BG, FG, 1, 1, true);
1132     }
1133     
1134     // Join slices
1135     typename ImageType::Pointer output;
1136     output = clitk::JoinSlices<ImageType>(o, input, d);
1137     return output;
1138   }
1139   //--------------------------------------------------------------------
1140
1141
1142   //--------------------------------------------------------------------
1143   template<class ImageType>
1144   typename ImageType::Pointer
1145   Clone(const ImageType * input) {
1146     typedef itk::ImageDuplicator<ImageType> DuplicatorType;
1147     typename DuplicatorType::Pointer duplicator = DuplicatorType::New();
1148     duplicator->SetInputImage(input);
1149     duplicator->Update();
1150     return duplicator->GetOutput();
1151   }
1152   //--------------------------------------------------------------------
1153
1154
1155   //--------------------------------------------------------------------
1156   /* Consider an input object, start at A, for each slice (dim1): 
1157      - compute the intersection between the AB line and the current slice
1158      - remove what is at lower or greater according to dim2 of this point
1159      - stop at B
1160   */
1161   template<class ImageType>
1162   typename ImageType::Pointer
1163   SliceBySliceSetBackgroundFromSingleLine(const ImageType * input, 
1164                                           typename ImageType::PixelType BG, 
1165                                           typename ImageType::PointType & A, 
1166                                           typename ImageType::PointType & B, 
1167                                           int dim1, int dim2, bool removeLowerPartFlag)
1168     
1169   {
1170     // Extract slices
1171     typedef typename itk::Image<typename ImageType::PixelType, ImageType::ImageDimension-1> SliceType;
1172     typedef typename SliceType::Pointer SlicePointer;
1173     std::vector<SlicePointer> slices;
1174     clitk::ExtractSlices<ImageType>(input, dim1, slices);
1175
1176     // Start at slice that contains A, and stop at B
1177     typename ImageType::IndexType Ap;
1178     typename ImageType::IndexType Bp;
1179     input->TransformPhysicalPointToIndex(A, Ap);
1180     input->TransformPhysicalPointToIndex(B, Bp);
1181     
1182     // Determine slice largest region
1183     typename SliceType::RegionType region = slices[0]->GetLargestPossibleRegion();
1184     typename SliceType::SizeType size = region.GetSize();
1185     typename SliceType::IndexType index = region.GetIndex();
1186
1187     // Line slope
1188     double a = (Bp[dim2]-Ap[dim2])/(Bp[dim1]-Ap[dim1]);
1189     double b = Ap[dim2];
1190
1191     // Loop from slice A to slice B
1192     for(uint i=0; i<(Bp[dim1]-Ap[dim1]); i++) {
1193       // Compute intersection between line AB and current slice for the dim2
1194       double p = a*i+b;
1195       // Change region (lower than dim2)
1196       if (removeLowerPartFlag) {
1197         size[dim2] = p-Ap[dim2];
1198       }
1199       else {
1200         size[dim2] = slices[0]->GetLargestPossibleRegion().GetSize()[dim2]-p;
1201         index[dim2] = p;
1202       }
1203       region.SetSize(size);
1204       region.SetIndex(index);
1205       // Fill region with BG (simple region iterator)
1206       FillRegionWithValue<SliceType>(slices[i+Ap[dim1]], BG, region);
1207       /*
1208       typedef itk::ImageRegionIterator<SliceType> IteratorType;
1209       IteratorType iter(slices[i+Ap[dim1]], region);
1210       iter.GoToBegin();
1211       while (!iter.IsAtEnd()) {
1212         iter.Set(BG);
1213         ++iter;
1214       }
1215       */
1216       // Loop
1217     }
1218     
1219     // Merge slices
1220     typename ImageType::Pointer output;
1221     output = clitk::JoinSlices<ImageType>(slices, input, dim1);
1222     return output;
1223   }
1224   //--------------------------------------------------------------------
1225
1226   //--------------------------------------------------------------------
1227   /* Consider an input object, slice by slice, use the point A and set
1228      pixel to BG according to their position relatively to A
1229   */
1230   template<class ImageType>
1231   typename ImageType::Pointer
1232   SliceBySliceSetBackgroundFromPoints(const ImageType * input, 
1233                                       typename ImageType::PixelType BG, 
1234                                       int sliceDim,
1235                                       std::vector<typename ImageType::PointType> & A, 
1236                                       bool removeGreaterThanXFlag,
1237                                       bool removeGreaterThanYFlag)
1238     
1239   {
1240     // Extract slices
1241     typedef typename itk::Image<typename ImageType::PixelType, ImageType::ImageDimension-1> SliceType;
1242     typedef typename SliceType::Pointer SlicePointer;
1243     std::vector<SlicePointer> slices;
1244     clitk::ExtractSlices<ImageType>(input, sliceDim, slices);
1245
1246     // Start at slice that contains A
1247     typename ImageType::IndexType Ap;
1248     
1249     // Determine slice largest region
1250     typename SliceType::RegionType region = slices[0]->GetLargestPossibleRegion();
1251     typename SliceType::SizeType size = region.GetSize();
1252     typename SliceType::IndexType index = region.GetIndex();
1253
1254     // Loop from slice A to slice B
1255     for(uint i=0; i<A.size(); i++) {
1256       input->TransformPhysicalPointToIndex(A[i], Ap);
1257       uint sliceIndex = Ap[2] - input->GetLargestPossibleRegion().GetIndex()[2];
1258       if ((sliceIndex < 0) || (sliceIndex >= slices.size())) {
1259         continue; // do not consider this slice
1260       }
1261       
1262       // Compute region for BG
1263       if (removeGreaterThanXFlag) {
1264         index[0] = Ap[0];
1265         size[0] = region.GetSize()[0]-(index[0]-region.GetIndex()[0]);
1266       }
1267       else {
1268         index[0] = region.GetIndex()[0];
1269         size[0] = Ap[0] - index[0];
1270       }
1271
1272       if (removeGreaterThanYFlag) {
1273         index[1] = Ap[1];
1274         size[1] = region.GetSize()[1]-(index[1]-region.GetIndex()[1]);
1275       }
1276       else {
1277         index[1] = region.GetIndex()[1];
1278         size[1] = Ap[1] - index[1];
1279       }
1280
1281       // Set region
1282       region.SetSize(size);
1283       region.SetIndex(index);
1284
1285       // Fill region with BG (simple region iterator)
1286       FillRegionWithValue<SliceType>(slices[sliceIndex], BG, region);
1287       // Loop
1288     }
1289     
1290     // Merge slices
1291     typename ImageType::Pointer output;
1292     output = clitk::JoinSlices<ImageType>(slices, input, sliceDim);
1293     return output;
1294   }
1295   //--------------------------------------------------------------------
1296
1297
1298   //--------------------------------------------------------------------
1299   template<class ImageType>
1300   void
1301   FillRegionWithValue(ImageType * input, typename ImageType::PixelType value, typename ImageType::RegionType & region)
1302   {
1303     typedef itk::ImageRegionIterator<ImageType> IteratorType;
1304     IteratorType iter(input, region);
1305     iter.GoToBegin();
1306     while (!iter.IsAtEnd()) {
1307       iter.Set(value);
1308       ++iter;
1309     }    
1310   }
1311   //--------------------------------------------------------------------
1312
1313
1314   //--------------------------------------------------------------------
1315   template<class ImageType>
1316   void
1317   GetMinMaxBoundary(ImageType * input, typename ImageType::PointType & min, 
1318                     typename ImageType::PointType & max)
1319   {
1320     typedef typename ImageType::PointType PointType;
1321     typedef typename ImageType::IndexType IndexType;
1322     IndexType min_i, max_i;
1323     min_i = input->GetLargestPossibleRegion().GetIndex();
1324     for(uint i=0; i<ImageType::ImageDimension; i++)
1325       max_i[i] = input->GetLargestPossibleRegion().GetSize()[i] + min_i[i];
1326     input->TransformIndexToPhysicalPoint(min_i, min);
1327     input->TransformIndexToPhysicalPoint(max_i, max);  
1328   }
1329   //--------------------------------------------------------------------
1330
1331
1332   //--------------------------------------------------------------------
1333   template<class ImageType>
1334   typename itk::Image<float, ImageType::ImageDimension>::Pointer
1335   DistanceMap(const ImageType * input, typename ImageType::PixelType BG)//, 
1336   //              typename itk::Image<float, ImageType::ImageDimension>::Pointer dmap) 
1337   {
1338     typedef itk::Image<float,ImageType::ImageDimension> FloatImageType;
1339     typedef itk::SignedMaurerDistanceMapImageFilter<ImageType, FloatImageType> DistanceMapFilterType;
1340     typename DistanceMapFilterType::Pointer filter = DistanceMapFilterType::New();
1341     filter->SetInput(input);
1342     filter->SetUseImageSpacing(true);
1343     filter->SquaredDistanceOff();
1344     filter->SetBackgroundValue(BG);
1345     filter->Update();
1346     return filter->GetOutput();
1347   }
1348   //--------------------------------------------------------------------
1349
1350
1351   //--------------------------------------------------------------------
1352   template<class ImageType>
1353   void 
1354   SliceBySliceBuildLineSegmentAccordingToMinimalDistanceBetweenStructures(const ImageType * S1, 
1355                                                                           const ImageType * S2, 
1356                                                                           typename ImageType::PixelType BG, 
1357                                                                           int sliceDimension, 
1358                                                                           std::vector<typename ImageType::PointType> & A, 
1359                                                                           std::vector<typename ImageType::PointType> & B)
1360   {
1361     // Extract slices
1362     typedef typename itk::Image<typename ImageType::PixelType, 2> SliceType;
1363     typedef typename SliceType::Pointer SlicePointer;
1364     std::vector<SlicePointer> slices_s1;
1365     std::vector<SlicePointer> slices_s2;
1366     clitk::ExtractSlices<ImageType>(S1, sliceDimension, slices_s1);
1367     clitk::ExtractSlices<ImageType>(S2, sliceDimension, slices_s2);
1368
1369     assert(slices_s1.size() == slices_s2.size());
1370
1371     // Prepare dmap
1372     typedef itk::Image<float,2> FloatImageType;
1373     typedef itk::SignedMaurerDistanceMapImageFilter<SliceType, FloatImageType> DistanceMapFilterType;
1374     std::vector<typename FloatImageType::Pointer> dmaps1;
1375     std::vector<typename FloatImageType::Pointer> dmaps2;
1376     typename FloatImageType::Pointer dmap;
1377
1378     // loop on slices
1379     for(uint i=0; i<slices_s1.size(); i++) {
1380       // Compute dmap for S1 *TO PUT IN FONCTION*
1381       dmap = clitk::DistanceMap<SliceType>(slices_s1[i], BG);
1382       dmaps1.push_back(dmap);
1383       writeImage<FloatImageType>(dmap, "dmap1.mha");
1384       // Compute dmap for S2
1385       dmap = clitk::DistanceMap<SliceType>(slices_s2[i], BG);
1386       dmaps2.push_back(dmap);
1387       writeImage<FloatImageType>(dmap, "dmap2.mha");
1388       
1389       // Look in S2 for the point the closest to S1
1390       typename SliceType::PointType p = ComputeClosestPoint<SliceType>(slices_s1[i], dmaps2[i], BG);
1391       typename ImageType::PointType p3D;
1392       clitk::PointsUtils<ImageType>::Convert2DTo3D(p, S1, i, p3D);
1393       A.push_back(p3D);
1394
1395       // Look in S2 for the point the closest to S1
1396       p = ComputeClosestPoint<SliceType>(slices_s2[i], dmaps1[i], BG);
1397       clitk::PointsUtils<ImageType>::Convert2DTo3D(p, S2, i, p3D);
1398       B.push_back(p3D);
1399
1400     }
1401
1402     // Debug dmap
1403     /*
1404       typedef itk::Image<float,3> FT;
1405       FT::Pointer f = FT::New();
1406       typename FT::Pointer d1 = clitk::JoinSlices<FT>(dmaps1, S1, 2);
1407       typename FT::Pointer d2 = clitk::JoinSlices<FT>(dmaps2, S2, 2);
1408       writeImage<FT>(d1, "d1.mha");
1409       writeImage<FT>(d2, "d2.mha");
1410     */
1411   }
1412   //--------------------------------------------------------------------
1413
1414
1415   //--------------------------------------------------------------------
1416   template<class ImageType>
1417   typename ImageType::PointType
1418   ComputeClosestPoint(const ImageType * input, 
1419                       const itk::Image<float, ImageType::ImageDimension> * dmap, 
1420                       typename ImageType::PixelType & BG) 
1421   {
1422     // Loop dmap + S2, if FG, get min
1423     typedef itk::Image<float,ImageType::ImageDimension> FloatImageType;
1424     typedef itk::ImageRegionConstIteratorWithIndex<ImageType> ImageIteratorType;
1425     typedef itk::ImageRegionConstIterator<FloatImageType> DMapIteratorType;
1426     ImageIteratorType iter1(input, input->GetLargestPossibleRegion());
1427     DMapIteratorType iter2(dmap, dmap->GetLargestPossibleRegion());
1428     
1429     iter1.GoToBegin();
1430     iter2.GoToBegin();
1431     double dmin = 100000.0;
1432     typename ImageType::IndexType indexmin;
1433     while (!iter1.IsAtEnd()) {
1434       if (iter1.Get() != BG) {
1435         double d = iter2.Get();
1436         if (d<dmin) {
1437           indexmin = iter1.GetIndex();
1438           dmin = d;
1439         }
1440       }
1441       ++iter1;
1442       ++iter2;
1443     }
1444     
1445     // Convert in Point
1446     typename ImageType::PointType p;
1447     input->TransformIndexToPhysicalPoint(indexmin, p);
1448     return p;
1449   }
1450   //--------------------------------------------------------------------
1451      
1452
1453
1454
1455 } // end of namespace
1456