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