]> Creatis software - clitk.git/blob - itk/clitkSegmentationUtils.txx
Merge branch 'master' of git.creatis.insa-lyon.fr:clitk
[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);
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     image->TransformIndexToPhysicalPoint(image->GetLargestPossibleRegion().GetIndex()+
387                                          image->GetLargestPossibleRegion().GetSize(), p);
388     // Add GetSpacing because remove Lower or equal than
389     // DD(max);
390     // DD(p);
391     // DD(max+image->GetSpacing()[dim]);
392     return CropImageAlongOneAxis<ImageType>(image, dim, max+image->GetSpacing()[dim], p[dim], autoCrop, BG);
393   }
394   //--------------------------------------------------------------------
395
396
397   //--------------------------------------------------------------------
398   template<class ImageType>
399   typename ImageType::Pointer
400   CropImageAlongOneAxis(const ImageType * image, 
401                         int dim, double min, double max, 
402                         bool autoCrop, typename ImageType::PixelType BG) 
403   {
404     // Compute region size
405     typename ImageType::RegionType region;
406     typename ImageType::SizeType size = image->GetLargestPossibleRegion().GetSize();
407     typename ImageType::PointType p = image->GetOrigin();
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     double m = image->GetOrigin()[dim] + size[dim]*image->GetSpacing()[dim];
412     if (max > m) p[dim] = m; // Check if not outside the image
413     else p[dim] = max;
414     typename ImageType::IndexType end;
415     image->TransformPhysicalPointToIndex(p, end);
416     size[dim] = abs(end[dim]-start[dim]);
417     region.SetIndex(start);
418     region.SetSize(size);
419   
420     // Perform Crop
421     typedef itk::RegionOfInterestImageFilter<ImageType, ImageType> CropFilterType;
422     typename CropFilterType::Pointer cropFilter = CropFilterType::New();
423     cropFilter->SetInput(image);
424     cropFilter->SetRegionOfInterest(region);
425     cropFilter->Update();
426     typename ImageType::Pointer result = cropFilter->GetOutput();
427   
428     // Auto Crop
429     if (autoCrop) {
430       result = AutoCrop<ImageType>(result, BG);
431     }
432     return result;
433   }
434   //--------------------------------------------------------------------
435
436
437   //--------------------------------------------------------------------
438   template<class ImageType>
439   void
440   ComputeCentroids(const ImageType * image, 
441                    typename ImageType::PixelType BG, 
442                    std::vector<typename ImageType::PointType> & centroids) 
443   {
444     typedef long LabelType;
445     static const unsigned int Dim = ImageType::ImageDimension;
446     typedef itk::ShapeLabelObject< LabelType, Dim > LabelObjectType;
447     typedef itk::LabelMap< LabelObjectType > LabelMapType;
448     typedef itk::LabelImageToLabelMapFilter<ImageType, LabelMapType> ImageToMapFilterType;
449     typename ImageToMapFilterType::Pointer imageToLabelFilter = ImageToMapFilterType::New(); 
450     typedef itk::ShapeLabelMapFilter<LabelMapType, ImageType> ShapeFilterType; 
451     typename ShapeFilterType::Pointer statFilter = ShapeFilterType::New();
452     imageToLabelFilter->SetBackgroundValue(BG);
453     imageToLabelFilter->SetInput(image);
454     statFilter->SetInput(imageToLabelFilter->GetOutput());
455     statFilter->Update();
456     typename LabelMapType::Pointer labelMap = statFilter->GetOutput();
457
458     centroids.clear();
459     typename ImageType::PointType dummy;
460     centroids.push_back(dummy); // label 0 -> no centroid, use dummy point for BG 
461     //DS FIXME (not useful ! to change ..)
462     for(uint i=0; i<labelMap->GetNumberOfLabelObjects(); i++) {
463       int label = labelMap->GetLabels()[i];
464       centroids.push_back(labelMap->GetLabelObject(label)->GetCentroid());
465     } 
466   }
467   //--------------------------------------------------------------------
468
469
470   //--------------------------------------------------------------------
471   template<class ImageType, class LabelType>
472   typename itk::LabelMap< itk::ShapeLabelObject<LabelType, ImageType::ImageDimension> >::Pointer
473   ComputeLabelMap(const ImageType * image, 
474                   typename ImageType::PixelType BG, 
475                   bool computePerimeterFlag) 
476   {
477     static const unsigned int Dim = ImageType::ImageDimension;
478     typedef itk::ShapeLabelObject< LabelType, Dim > LabelObjectType;
479     typedef itk::LabelMap< LabelObjectType > LabelMapType;
480     typedef itk::LabelImageToLabelMapFilter<ImageType, LabelMapType> ImageToMapFilterType;
481     typename ImageToMapFilterType::Pointer imageToLabelFilter = ImageToMapFilterType::New(); 
482     typedef itk::ShapeLabelMapFilter<LabelMapType, ImageType> ShapeFilterType; 
483     typename ShapeFilterType::Pointer statFilter = ShapeFilterType::New();
484     imageToLabelFilter->SetBackgroundValue(BG);
485     imageToLabelFilter->SetInput(image);
486     statFilter->SetInput(imageToLabelFilter->GetOutput());
487     statFilter->SetComputePerimeter(computePerimeterFlag);
488     statFilter->Update();
489     return statFilter->GetOutput();
490   }
491   //--------------------------------------------------------------------
492
493
494   //--------------------------------------------------------------------
495   template<class ImageType>
496   void
497   ComputeCentroids2(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
518     for(uint i=1; i<labelMap->GetNumberOfLabelObjects()+1; i++) {
519       centroids.push_back(labelMap->GetLabelObject(i)->GetCentroid());
520     } 
521     
522     for(uint i=1; i<labelMap->GetNumberOfLabelObjects()+1; i++) {
523       DD(labelMap->GetLabelObject(i)->GetBinaryPrincipalAxes());
524       DD(labelMap->GetLabelObject(i)->GetBinaryFlatness());
525       DD(labelMap->GetLabelObject(i)->GetRoundness ());      
526
527       // search for the point on the boundary alog PA
528
529     }
530
531   }
532   //--------------------------------------------------------------------
533
534
535   //--------------------------------------------------------------------
536   template<class ImageType>
537   void
538   PointsUtils<ImageType>::Convert2DTo3D(const PointType2D & p2D, 
539                                         const ImageType * image, 
540                                         const int slice, 
541                                         PointType3D & p3D)  
542   {
543     IndexType3D index3D;
544     index3D[0] = index3D[1] = 0;
545     index3D[2] = image->GetLargestPossibleRegion().GetIndex()[2]+slice;
546     image->TransformIndexToPhysicalPoint(index3D, p3D);
547     p3D[0] = p2D[0]; 
548     p3D[1] = p2D[1];
549     //  p3D[2] = p[2];//(image->GetLargestPossibleRegion().GetIndex()[2]+slice)*image->GetSpacing()[2] 
550     //    + image->GetOrigin()[2];
551   }
552   //--------------------------------------------------------------------
553
554
555   //--------------------------------------------------------------------
556   template<class ImageType>
557   void 
558   PointsUtils<ImageType>::Convert2DMapTo3DList(const MapPoint2DType & map, 
559                                             const ImageType * image, 
560                                             VectorPoint3DType & list)
561   {
562     typename MapPoint2DType::const_iterator iter = map.begin();
563     while (iter != map.end()) {
564       PointType3D p;
565       Convert2DTo3D(iter->second, image, iter->first, p);
566       list.push_back(p);
567       ++iter;
568     }
569   }
570   //--------------------------------------------------------------------
571
572
573   //--------------------------------------------------------------------
574   template<class ImageType>
575   void 
576   PointsUtils<ImageType>::Convert2DListTo3DList(const VectorPoint2DType & p2D, 
577                                                 int slice,
578                                                 const ImageType * image, 
579                                                 VectorPoint3DType & list) 
580   {
581     for(uint i=0; i<p2D.size(); i++) {
582       PointType3D p;
583       Convert2DTo3D(p2D[i], image, slice, p);
584       list.push_back(p);
585     }
586   }
587   //--------------------------------------------------------------------
588
589
590   //--------------------------------------------------------------------
591   template<class ImageType>
592   void 
593   WriteListOfLandmarks(std::vector<typename ImageType::PointType> points, 
594                        std::string filename)
595   {
596     std::ofstream os; 
597     openFileForWriting(os, filename); 
598     os << "LANDMARKS1" << std::endl;  
599     for(uint i=0; i<points.size(); i++) {
600       const typename ImageType::PointType & p = points[i];
601       // Write it in the file
602       os << i << " " << p[0] << " " << p[1] << " " << p[2] << " 0 0 " << std::endl;
603     }
604     os.close();
605   }
606   //--------------------------------------------------------------------
607
608
609   //--------------------------------------------------------------------
610   template<class ImageType>
611   typename ImageType::Pointer 
612   Dilate(const ImageType * image, double radiusInMM,               
613          typename ImageType::PixelType BG,
614          typename ImageType::PixelType FG,  
615          bool extendSupport)
616   {
617     typename ImageType::SizeType r;
618     for(uint i=0; i<ImageType::ImageDimension; i++) 
619       r[i] = (uint)lrint(radiusInMM/image->GetSpacing()[i]);
620     return Dilate<ImageType>(image, r, BG, FG, extendSupport);
621   }
622   //--------------------------------------------------------------------
623
624
625   //--------------------------------------------------------------------
626   template<class ImageType>
627   typename ImageType::Pointer 
628   Dilate(const ImageType * image, typename ImageType::PointType radiusInMM, 
629          typename ImageType::PixelType BG, 
630          typename ImageType::PixelType FG, 
631          bool extendSupport)
632   {
633     typename ImageType::SizeType r;
634     for(uint i=0; i<ImageType::ImageDimension; i++) 
635       r[i] = (uint)lrint(radiusInMM[i]/image->GetSpacing()[i]);
636     return Dilate<ImageType>(image, r, BG, FG, extendSupport);
637   }
638   //--------------------------------------------------------------------
639
640
641   //--------------------------------------------------------------------
642   template<class ImageType>
643   typename ImageType::Pointer 
644   Dilate(const ImageType * image, typename ImageType::SizeType radius, 
645          typename ImageType::PixelType BG, 
646          typename ImageType::PixelType FG, 
647          bool extendSupport)
648   {
649     // Create kernel for dilatation
650     typedef itk::BinaryBallStructuringElement<typename ImageType::PixelType, 
651                                               ImageType::ImageDimension> KernelType;
652     KernelType structuringElement;
653     structuringElement.SetRadius(radius);
654     structuringElement.CreateStructuringElement();
655
656     typename ImageType::Pointer output;
657     if (extendSupport) {
658       typedef itk::ConstantPadImageFilter<ImageType, ImageType> PadFilterType;
659       typename PadFilterType::Pointer padFilter = PadFilterType::New();
660       padFilter->SetInput(image);
661       typename ImageType::SizeType lower;
662       typename ImageType::SizeType upper;
663       for(uint i=0; i<3; i++) {
664         lower[i] = upper[i] = 2*(radius[i]+1);
665       }
666       padFilter->SetPadLowerBound(lower);
667       padFilter->SetPadUpperBound(upper);
668       padFilter->Update();
669       output = padFilter->GetOutput();
670     }
671
672     // Dilate  filter
673     typedef itk::BinaryDilateImageFilter<ImageType, ImageType , KernelType> DilateFilterType;
674     typename DilateFilterType::Pointer dilateFilter = DilateFilterType::New();
675     dilateFilter->SetBackgroundValue(BG);
676     dilateFilter->SetForegroundValue(FG);
677     dilateFilter->SetBoundaryToForeground(false);
678     dilateFilter->SetKernel(structuringElement);
679     if (extendSupport) dilateFilter->SetInput(output);
680     else dilateFilter->SetInput(image);
681     dilateFilter->Update();
682     return dilateFilter->GetOutput();
683   }
684   //--------------------------------------------------------------------
685
686
687   //--------------------------------------------------------------------
688   template<class ImageType>
689   typename ImageType::Pointer 
690   Opening(const ImageType * image, typename ImageType::SizeType radius,
691          typename ImageType::PixelType BG,
692          typename ImageType::PixelType FG)
693   {
694     // Kernel 
695     typedef itk::BinaryBallStructuringElement<typename ImageType::PixelType, 
696                                               ImageType::ImageDimension> KernelType;    
697     KernelType structuringElement;
698     structuringElement.SetRadius(radius);
699     structuringElement.CreateStructuringElement();
700     
701     // Filter
702     typedef itk::BinaryMorphologicalOpeningImageFilter<ImageType, ImageType , KernelType> OpeningFilterType;
703     typename OpeningFilterType::Pointer open = OpeningFilterType::New();
704     open->SetInput(image);
705     open->SetBackgroundValue(BG);
706     open->SetForegroundValue(FG);
707     open->SetKernel(structuringElement);
708     open->Update();
709     return open->GetOutput();
710   }
711   //--------------------------------------------------------------------
712
713
714
715   //--------------------------------------------------------------------
716   template<class ValueType, class VectorType>
717   void ConvertOption(std::string optionName, uint given, 
718                      ValueType * values, VectorType & p, 
719                      uint dim, bool required) 
720   {
721     if (required && (given == 0)) {
722       clitkExceptionMacro("The option --" << optionName << " must be set and have 1 or " 
723                           << dim << " values.");
724     }
725     if (given == 1) {
726       for(uint i=0; i<dim; i++) p[i] = values[0];
727       return;
728     }
729     if (given == dim) {
730       for(uint i=0; i<dim; i++) p[i] = values[i];
731       return;
732     }
733     if (given == 0) return;
734     clitkExceptionMacro("The option --" << optionName << " must have 1 or " 
735                         << dim << " values.");
736   }
737   //--------------------------------------------------------------------
738
739
740   //--------------------------------------------------------------------
741   /*
742     http://www.gamedev.net/community/forums/topic.asp?topic_id=542870
743     Assuming the points are (Ax,Ay) (Bx,By) and (Cx,Cy), you need to compute:
744     (Bx - Ax) * (Cy - Ay) - (By - Ay) * (Cx - Ax)
745     This will equal zero if the point C is on the line formed by
746     points A and B, and will have a different sign depending on the
747     side. Which side this is depends on the orientation of your (x,y)
748     coordinates, but you can plug test values for A,B and C into this
749     formula to determine whether negative values are to the left or to
750     the right.
751     => to accelerate, start with formula, when change sign -> stop and fill
752
753     offsetToKeep = is used to determine which side of the line we
754     keep. The point along the mainDirection but 'offsetToKeep' mm away
755     is kept.
756   
757   */
758   template<class ImageType>
759   void 
760   SliceBySliceSetBackgroundFromLineSeparation(ImageType * input, 
761                                               std::vector<typename ImageType::PointType> & lA, 
762                                               std::vector<typename ImageType::PointType> & lB, 
763                                               typename ImageType::PixelType BG, 
764                                               int mainDirection, 
765                                               double offsetToKeep)
766   {
767     assert((mainDirection==0) || (mainDirection==1));
768     typedef itk::ImageSliceIteratorWithIndex<ImageType> SliceIteratorType;
769     SliceIteratorType siter = SliceIteratorType(input, input->GetLargestPossibleRegion());
770     siter.SetFirstDirection(0);
771     siter.SetSecondDirection(1);
772     siter.GoToBegin();
773     uint i=0;
774     typename ImageType::PointType A;
775     typename ImageType::PointType B;
776     typename ImageType::PointType C;
777     assert(lA.size() == lB.size());
778     while ((i<lA.size()) && (!siter.IsAtEnd())) {
779       // Check that the current slice correspond to the current point
780       input->TransformIndexToPhysicalPoint(siter.GetIndex(), C);
781       if ((fabs(C[2] - lA[i][2]))>0.01) { // is !equal with a tolerance of 0.01 mm
782       }
783       else {
784         // Define A,B,C points
785         A = lA[i];
786         B = lB[i];
787         C = A;
788         // Check that the line is not a point (A=B)
789         bool p = (A[0] == B[0]) && (A[1] == B[1]);
790       
791         if (!p) {
792           C[mainDirection] += offsetToKeep; // I know I must keep this point
793           double s = (B[0] - A[0]) * (C[1] - A[1]) - (B[1] - A[1]) * (C[0] - A[0]);
794           bool isPositive = s<0;
795           while (!siter.IsAtEndOfSlice()) {
796             while (!siter.IsAtEndOfLine()) {
797               // Very slow, I know ... but image should be very small
798               input->TransformIndexToPhysicalPoint(siter.GetIndex(), C);
799               double s = (B[0] - A[0]) * (C[1] - A[1]) - (B[1] - A[1]) * (C[0] - A[0]);
800               if (s == 0) siter.Set(BG); // on the line, we decide to remove
801               if (isPositive) {
802                 if (s > 0) siter.Set(BG);
803               }
804               else {
805                 if (s < 0) siter.Set(BG); 
806               }
807               ++siter;
808             }
809             siter.NextLine();
810           } // end loop slice
811         }      
812
813         ++i;
814       } // End of current slice
815       siter.NextSlice();
816     }
817   }                                                   
818   //--------------------------------------------------------------------
819
820
821   //--------------------------------------------------------------------
822   template<class ImageType>
823   void 
824   AndNot(ImageType * input, 
825          const ImageType * object, 
826          typename ImageType::PixelType BG)
827   {
828     typename ImageType::Pointer o;
829     bool resized=false;
830     if (!clitk::HaveSameSizeAndSpacing<ImageType, ImageType>(input, object)) {
831       o = clitk::ResizeImageLike<ImageType>(object, input, BG);
832       resized = true;
833     }
834
835     typedef clitk::BooleanOperatorLabelImageFilter<ImageType> BoolFilterType;
836     typename BoolFilterType::Pointer boolFilter = BoolFilterType::New(); 
837     boolFilter->InPlaceOn();
838     boolFilter->SetInput1(input);
839     if (resized) boolFilter->SetInput2(o);  
840     else boolFilter->SetInput2(object);
841     boolFilter->SetBackgroundValue1(BG);
842     boolFilter->SetBackgroundValue2(BG);
843     boolFilter->SetOperationType(BoolFilterType::AndNot);
844     boolFilter->Update();
845   }
846   //--------------------------------------------------------------------
847
848
849   //--------------------------------------------------------------------
850   template<class ImageType>
851   void 
852   And(ImageType * input, 
853       const ImageType * object, 
854       typename ImageType::PixelType BG)
855   {
856     typename ImageType::Pointer o;
857     bool resized=false;
858     if (!clitk::HaveSameSizeAndSpacing<ImageType, ImageType>(input, object)) {
859       o = clitk::ResizeImageLike<ImageType>(object, input, BG);
860       resized = true;
861     }
862
863     typedef clitk::BooleanOperatorLabelImageFilter<ImageType> BoolFilterType;
864     typename BoolFilterType::Pointer boolFilter = BoolFilterType::New(); 
865     boolFilter->InPlaceOn();
866     boolFilter->SetInput1(input);
867     if (resized) boolFilter->SetInput2(o);  
868     else boolFilter->SetInput2(object);
869     boolFilter->SetBackgroundValue1(BG);
870     boolFilter->SetBackgroundValue2(BG);
871     boolFilter->SetOperationType(BoolFilterType::And);
872     boolFilter->Update();
873   }
874   //--------------------------------------------------------------------
875
876
877   //--------------------------------------------------------------------
878   template<class ImageType>
879   void 
880   Or(ImageType * input, 
881      const ImageType * object, 
882      typename ImageType::PixelType BG)
883   {
884     typename ImageType::Pointer o;
885     bool resized=false;
886     if (!clitk::HaveSameSizeAndSpacing<ImageType, ImageType>(input, object)) {
887       o = clitk::ResizeImageLike<ImageType>(object, input, BG);
888       resized = true;
889     }
890
891     typedef clitk::BooleanOperatorLabelImageFilter<ImageType> BoolFilterType;
892     typename BoolFilterType::Pointer boolFilter = BoolFilterType::New(); 
893     boolFilter->InPlaceOn();
894     boolFilter->SetInput1(input);
895     if (resized) boolFilter->SetInput2(o);  
896     else boolFilter->SetInput2(object);
897     boolFilter->SetBackgroundValue1(BG);
898     boolFilter->SetBackgroundValue2(BG);
899     boolFilter->SetOperationType(BoolFilterType::Or);
900     boolFilter->Update();
901   }
902   //--------------------------------------------------------------------
903
904
905   //--------------------------------------------------------------------
906   template<class ImageType>
907   typename ImageType::Pointer
908   Binarize(const ImageType * input, 
909            typename ImageType::PixelType lower, 
910            typename ImageType::PixelType upper, 
911            typename ImageType::PixelType BG,
912            typename ImageType::PixelType FG) 
913   {
914     typedef itk::BinaryThresholdImageFilter<ImageType, ImageType> BinaryThresholdFilterType;
915     typename BinaryThresholdFilterType::Pointer binarizeFilter = BinaryThresholdFilterType::New();
916     binarizeFilter->SetInput(input);
917     binarizeFilter->InPlaceOff();
918     binarizeFilter->SetLowerThreshold(lower);
919     binarizeFilter->SetUpperThreshold(upper);
920     binarizeFilter->SetInsideValue(FG);
921     binarizeFilter->SetOutsideValue(BG);
922     binarizeFilter->Update();
923     return binarizeFilter->GetOutput();
924   }
925   //--------------------------------------------------------------------
926
927
928   //--------------------------------------------------------------------
929   template<class ImageType>
930   void
931   GetMinMaxPointPosition(const ImageType * input, 
932                          typename ImageType::PointType & min,
933                          typename ImageType::PointType & max) 
934   {
935     typename ImageType::IndexType index = input->GetLargestPossibleRegion().GetIndex();
936     input->TransformIndexToPhysicalPoint(index, min);
937     index = index+input->GetLargestPossibleRegion().GetSize();
938     input->TransformIndexToPhysicalPoint(index, max);
939   }
940   //--------------------------------------------------------------------
941
942
943   //--------------------------------------------------------------------
944   template<class ImageType>
945   typename ImageType::PointType
946   FindExtremaPointInAGivenLine(const ImageType * input, 
947                                int dimension, 
948                                bool inverse, 
949                                typename ImageType::PointType p, 
950                                typename ImageType::PixelType BG, 
951                                double distanceMax) 
952   {
953     // Which direction ?  Increasing or decreasing.
954     int d=1;
955     if (inverse) d=-1;
956   
957     // Transform to pixel index
958     typename ImageType::IndexType index;
959     input->TransformPhysicalPointToIndex(p, index);
960
961     // Loop while inside the mask;
962     while (input->GetPixel(index) != BG) {
963       index[dimension] += d;
964     }
965
966     // Transform back to Physical Units
967     typename ImageType::PointType result;
968     input->TransformIndexToPhysicalPoint(index, result);
969
970     // Check that is is not too far away
971     double distance = p.EuclideanDistanceTo(result);
972     if (distance > distanceMax) {
973       result = p; // Get back to initial value
974     }
975
976     return result;
977   }
978   //--------------------------------------------------------------------
979
980
981   //--------------------------------------------------------------------
982   template<class PointType>
983   bool
984   IsOnTheSameLineSide(PointType C, PointType A, PointType B, PointType like) 
985   {
986     // Look at the position of point 'like' according to the AB line
987     double s = (B[0] - A[0]) * (like[1] - A[1]) - (B[1] - A[1]) * (like[0] - A[0]);
988     bool negative = s<0;
989   
990     // Look the C position
991     s = (B[0] - A[0]) * (C[1] - A[1]) - (B[1] - A[1]) * (C[0] - A[0]);
992
993     if (negative && (s<=0)) return true;
994     if (!negative && (s>=0)) return true;
995     return false;
996   }
997   //--------------------------------------------------------------------
998
999
1000   //--------------------------------------------------------------------
1001   /* Consider an input object, for each slice, find the extrema
1002      position according to a given direction and build a line segment
1003      passing throught this point in a given direction.  Output is a
1004      vector of line (from point A to B), for each slice;
1005    */
1006   template<class ImageType>
1007   void 
1008   SliceBySliceBuildLineSegmentAccordingToExtremaPosition(const ImageType * input, 
1009                                                          typename ImageType::PixelType BG, 
1010                                                          int sliceDimension, 
1011                                                          int extremaDirection, 
1012                                                          bool extremaOppositeFlag, 
1013                                                          int lineDirection,
1014                                                          double margin,
1015                                                          std::vector<typename ImageType::PointType> & A, 
1016                                                          std::vector<typename ImageType::PointType> & B)
1017   {
1018     // Type of a slice
1019     typedef typename itk::Image<typename ImageType::PixelType, ImageType::ImageDimension-1> SliceType;
1020     
1021     // Build the list of slices
1022     std::vector<typename SliceType::Pointer> slices;
1023     clitk::ExtractSlices<ImageType>(input, sliceDimension, slices);
1024
1025     // Build the list of 2D points
1026     std::map<int, typename SliceType::PointType> position2D;
1027     for(uint i=0; i<slices.size(); i++) {
1028       typename SliceType::PointType p;
1029       bool found = 
1030         clitk::FindExtremaPointInAGivenDirection<SliceType>(slices[i], BG, 
1031                                                             extremaDirection, extremaOppositeFlag, p);
1032       if (found) {
1033         position2D[i] = p;
1034       }
1035     }
1036     
1037     // Convert 2D points in slice into 3D points
1038     clitk::PointsUtils<ImageType>::Convert2DMapTo3DList(position2D, input, A);
1039     
1040     // Create additional point just right to the previous ones, on the
1041     // given lineDirection, in order to create a horizontal/vertical line.
1042     for(uint i=0; i<A.size(); i++) {
1043       typename ImageType::PointType p = A[i];
1044       p[lineDirection] += 10;
1045       B.push_back(p);
1046       // Margins ?
1047       A[i][extremaDirection] += margin;
1048       B[i][extremaDirection] += margin;
1049     }
1050
1051   }
1052   //--------------------------------------------------------------------
1053
1054
1055   //--------------------------------------------------------------------
1056   template<class ImageType>
1057   typename ImageType::Pointer
1058   SliceBySliceKeepMainCCL(const ImageType * input, 
1059                           typename ImageType::PixelType BG,
1060                           typename ImageType::PixelType FG)  {
1061     
1062     // Extract slices
1063     const int d = ImageType::ImageDimension-1;
1064     typedef typename itk::Image<typename ImageType::PixelType, d> SliceType;
1065     std::vector<typename SliceType::Pointer> slices;
1066     clitk::ExtractSlices<ImageType>(input, d, slices);
1067     
1068     // Labelize and keep the main one
1069     std::vector<typename SliceType::Pointer> o;
1070     for(uint i=0; i<slices.size(); i++) {
1071       o.push_back(clitk::Labelize<SliceType>(slices[i], BG, false, 1));
1072       o[i] = clitk::KeepLabels<SliceType>(o[i], BG, FG, 1, 1, true);
1073     }
1074     
1075     // Join slices
1076     typename ImageType::Pointer output;
1077     output = clitk::JoinSlices<ImageType>(o, input, d);
1078     return output;
1079   }
1080   //--------------------------------------------------------------------
1081
1082
1083   //--------------------------------------------------------------------
1084   template<class ImageType>
1085   typename ImageType::Pointer
1086   Clone(const ImageType * input) {
1087     typedef itk::ImageDuplicator<ImageType> DuplicatorType;
1088     typename DuplicatorType::Pointer duplicator = DuplicatorType::New();
1089     duplicator->SetInputImage(input);
1090     duplicator->Update();
1091     return duplicator->GetOutput();
1092   }
1093   //--------------------------------------------------------------------
1094
1095
1096   //--------------------------------------------------------------------
1097   /* Consider an input object, start at A, for each slice (dim1): 
1098      - compute the intersection between the AB line and the current slice
1099      - remove what is at lower or greater according to dim2 of this point
1100      - stop at B
1101   */
1102   template<class ImageType>
1103   typename ImageType::Pointer
1104   SliceBySliceSetBackgroundFromSingleLine(const ImageType * input, 
1105                                           typename ImageType::PixelType BG, 
1106                                           typename ImageType::PointType & A, 
1107                                           typename ImageType::PointType & B, 
1108                                           int dim1, int dim2, bool removeLowerPartFlag)
1109     
1110   {
1111     // Extract slices
1112     typedef typename itk::Image<typename ImageType::PixelType, ImageType::ImageDimension-1> SliceType;
1113     typedef typename SliceType::Pointer SlicePointer;
1114     std::vector<SlicePointer> slices;
1115     clitk::ExtractSlices<ImageType>(input, dim1, slices);
1116
1117     // Start at slice that contains A, and stop at B
1118     typename ImageType::IndexType Ap;
1119     typename ImageType::IndexType Bp;
1120     input->TransformPhysicalPointToIndex(A, Ap);
1121     input->TransformPhysicalPointToIndex(B, Bp);
1122     
1123     // Determine slice largest region
1124     typename SliceType::RegionType region = slices[0]->GetLargestPossibleRegion();
1125     typename SliceType::SizeType size = region.GetSize();
1126     typename SliceType::IndexType index = region.GetIndex();
1127
1128     // Line slope
1129     double a = (Bp[dim2]-Ap[dim2])/(Bp[dim1]-Ap[dim1]);
1130     double b = Ap[dim2];
1131
1132     // Loop from slice A to slice B
1133     for(uint i=0; i<(Bp[dim1]-Ap[dim1]); i++) {
1134       // Compute intersection between line AB and current slice for the dim2
1135       double p = a*i+b;
1136       // Change region (lower than dim2)
1137       if (removeLowerPartFlag) {
1138         size[dim2] = p-Ap[dim2];
1139       }
1140       else {
1141         size[dim2] = slices[0]->GetLargestPossibleRegion().GetSize()[dim2]-p;
1142         index[dim2] = p;
1143       }
1144       region.SetSize(size);
1145       region.SetIndex(index);
1146       // Fill region with BG (simple region iterator)
1147       FillRegionWithValue<SliceType>(slices[i+Ap[dim1]], BG, region);
1148       /*
1149       typedef itk::ImageRegionIterator<SliceType> IteratorType;
1150       IteratorType iter(slices[i+Ap[dim1]], region);
1151       iter.GoToBegin();
1152       while (!iter.IsAtEnd()) {
1153         iter.Set(BG);
1154         ++iter;
1155       }
1156       */
1157       // Loop
1158     }
1159     
1160     // Merge slices
1161     typename ImageType::Pointer output;
1162     output = clitk::JoinSlices<ImageType>(slices, input, dim1);
1163     return output;
1164   }
1165   //--------------------------------------------------------------------
1166
1167   //--------------------------------------------------------------------
1168   /* Consider an input object, slice by slice, use the point A and set
1169      pixel to BG according to their position relatively to A
1170   */
1171   template<class ImageType>
1172   typename ImageType::Pointer
1173   SliceBySliceSetBackgroundFromPoints(const ImageType * input, 
1174                                       typename ImageType::PixelType BG, 
1175                                       int sliceDim,
1176                                       std::vector<typename ImageType::PointType> & A, 
1177                                       bool removeGreaterThanXFlag,
1178                                       bool removeGreaterThanYFlag)
1179     
1180   {
1181     // Extract slices
1182     typedef typename itk::Image<typename ImageType::PixelType, ImageType::ImageDimension-1> SliceType;
1183     typedef typename SliceType::Pointer SlicePointer;
1184     std::vector<SlicePointer> slices;
1185     clitk::ExtractSlices<ImageType>(input, sliceDim, slices);
1186
1187     // Start at slice that contains A
1188     typename ImageType::IndexType Ap;
1189     
1190     // Determine slice largest region
1191     typename SliceType::RegionType region = slices[0]->GetLargestPossibleRegion();
1192     typename SliceType::SizeType size = region.GetSize();
1193     typename SliceType::IndexType index = region.GetIndex();
1194
1195     // Loop from slice A to slice B
1196     for(uint i=0; i<A.size(); i++) {
1197       input->TransformPhysicalPointToIndex(A[i], Ap);
1198       uint sliceIndex = Ap[2] - input->GetLargestPossibleRegion().GetIndex()[2];
1199       if ((sliceIndex < 0) || (sliceIndex >= slices.size())) {
1200         continue; // do not consider this slice
1201       }
1202       
1203       // Compute region for BG
1204       if (removeGreaterThanXFlag) {
1205         index[0] = Ap[0];
1206         size[0] = region.GetSize()[0]-(index[0]-region.GetIndex()[0]);
1207       }
1208       else {
1209         index[0] = region.GetIndex()[0];
1210         size[0] = Ap[0] - index[0];
1211       }
1212
1213       if (removeGreaterThanYFlag) {
1214         index[1] = Ap[1];
1215         size[1] = region.GetSize()[1]-(index[1]-region.GetIndex()[1]);
1216       }
1217       else {
1218         index[1] = region.GetIndex()[1];
1219         size[1] = Ap[1] - index[1];
1220       }
1221
1222       // Set region
1223       region.SetSize(size);
1224       region.SetIndex(index);
1225
1226       // Fill region with BG (simple region iterator)
1227       FillRegionWithValue<SliceType>(slices[sliceIndex], BG, region);
1228       // Loop
1229     }
1230     
1231     // Merge slices
1232     typename ImageType::Pointer output;
1233     output = clitk::JoinSlices<ImageType>(slices, input, sliceDim);
1234     return output;
1235   }
1236   //--------------------------------------------------------------------
1237
1238
1239   //--------------------------------------------------------------------
1240   template<class ImageType>
1241   void
1242   FillRegionWithValue(ImageType * input, typename ImageType::PixelType value, typename ImageType::RegionType & region)
1243   {
1244     typedef itk::ImageRegionIterator<ImageType> IteratorType;
1245     IteratorType iter(input, region);
1246     iter.GoToBegin();
1247     while (!iter.IsAtEnd()) {
1248       iter.Set(value);
1249       ++iter;
1250     }    
1251   }
1252   //--------------------------------------------------------------------
1253
1254
1255   //--------------------------------------------------------------------
1256   template<class ImageType>
1257   void
1258   GetMinMaxBoundary(ImageType * input, typename ImageType::PointType & min, 
1259                     typename ImageType::PointType & max)
1260   {
1261     typedef typename ImageType::PointType PointType;
1262     typedef typename ImageType::IndexType IndexType;
1263     IndexType min_i, max_i;
1264     min_i = input->GetLargestPossibleRegion().GetIndex();
1265     for(uint i=0; i<ImageType::ImageDimension; i++)
1266       max_i[i] = input->GetLargestPossibleRegion().GetSize()[i] + min_i[i];
1267     input->TransformIndexToPhysicalPoint(min_i, min);
1268     input->TransformIndexToPhysicalPoint(max_i, max);  
1269   }
1270   //--------------------------------------------------------------------
1271
1272
1273   //--------------------------------------------------------------------
1274   template<class ImageType>
1275   typename itk::Image<float, ImageType::ImageDimension>::Pointer
1276   DistanceMap(const ImageType * input, typename ImageType::PixelType BG)//, 
1277   //              typename itk::Image<float, ImageType::ImageDimension>::Pointer dmap) 
1278   {
1279     typedef itk::Image<float,ImageType::ImageDimension> FloatImageType;
1280     typedef itk::SignedMaurerDistanceMapImageFilter<ImageType, FloatImageType> DistanceMapFilterType;
1281     typename DistanceMapFilterType::Pointer filter = DistanceMapFilterType::New();
1282     filter->SetInput(input);
1283     filter->SetUseImageSpacing(true);
1284     filter->SquaredDistanceOff();
1285     filter->SetBackgroundValue(BG);
1286     filter->Update();
1287     return filter->GetOutput();
1288   }
1289   //--------------------------------------------------------------------
1290
1291
1292   //--------------------------------------------------------------------
1293   template<class ImageType>
1294   void 
1295   SliceBySliceBuildLineSegmentAccordingToMinimalDistanceBetweenStructures(const ImageType * S1, 
1296                                                                           const ImageType * S2, 
1297                                                                           typename ImageType::PixelType BG, 
1298                                                                           int sliceDimension, 
1299                                                                           std::vector<typename ImageType::PointType> & A, 
1300                                                                           std::vector<typename ImageType::PointType> & B)
1301   {
1302     // Extract slices
1303     typedef typename itk::Image<typename ImageType::PixelType, 2> SliceType;
1304     typedef typename SliceType::Pointer SlicePointer;
1305     std::vector<SlicePointer> slices_s1;
1306     std::vector<SlicePointer> slices_s2;
1307     clitk::ExtractSlices<ImageType>(S1, sliceDimension, slices_s1);
1308     clitk::ExtractSlices<ImageType>(S2, sliceDimension, slices_s2);
1309
1310     assert(slices_s1.size() == slices_s2.size());
1311
1312     // Prepare dmap
1313     typedef itk::Image<float,2> FloatImageType;
1314     typedef itk::SignedMaurerDistanceMapImageFilter<SliceType, FloatImageType> DistanceMapFilterType;
1315     std::vector<typename FloatImageType::Pointer> dmaps1;
1316     std::vector<typename FloatImageType::Pointer> dmaps2;
1317     typename FloatImageType::Pointer dmap;
1318
1319     // loop on slices
1320     for(uint i=0; i<slices_s1.size(); i++) {
1321       // Compute dmap for S1 *TO PUT IN FONCTION*
1322       dmap = clitk::DistanceMap<SliceType>(slices_s1[i], BG);
1323       dmaps1.push_back(dmap);
1324       //writeImage<FloatImageType>(dmap, "dmap1.mha");
1325       // Compute dmap for S2
1326       dmap = clitk::DistanceMap<SliceType>(slices_s2[i], BG);
1327       dmaps2.push_back(dmap);
1328       //writeImage<FloatImageType>(dmap, "dmap2.mha");
1329       
1330       // Look in S2 for the point the closest to S1
1331       typename SliceType::PointType p = ComputeClosestPoint<SliceType>(slices_s1[i], dmaps2[i], BG);
1332       typename ImageType::PointType p3D;
1333       clitk::PointsUtils<ImageType>::Convert2DTo3D(p, S1, i, p3D);
1334       A.push_back(p3D);
1335
1336       // Look in S2 for the point the closest to S1
1337       p = ComputeClosestPoint<SliceType>(slices_s2[i], dmaps1[i], BG);
1338       clitk::PointsUtils<ImageType>::Convert2DTo3D(p, S2, i, p3D);
1339       B.push_back(p3D);
1340
1341     }
1342
1343     // Debug dmap
1344     /*
1345       typedef itk::Image<float,3> FT;
1346       FT::Pointer f = FT::New();
1347       typename FT::Pointer d1 = clitk::JoinSlices<FT>(dmaps1, S1, 2);
1348       typename FT::Pointer d2 = clitk::JoinSlices<FT>(dmaps2, S2, 2);
1349       writeImage<FT>(d1, "d1.mha");
1350       writeImage<FT>(d2, "d2.mha");
1351     */
1352   }
1353   //--------------------------------------------------------------------
1354
1355
1356   //--------------------------------------------------------------------
1357   template<class ImageType>
1358   typename ImageType::PointType
1359   ComputeClosestPoint(const ImageType * input, 
1360                       const itk::Image<float, ImageType::ImageDimension> * dmap, 
1361                       typename ImageType::PixelType & BG) 
1362   {
1363     // Loop dmap + S2, if FG, get min
1364     typedef itk::Image<float,ImageType::ImageDimension> FloatImageType;
1365     typedef itk::ImageRegionConstIteratorWithIndex<ImageType> ImageIteratorType;
1366     typedef itk::ImageRegionConstIterator<FloatImageType> DMapIteratorType;
1367     ImageIteratorType iter1(input, input->GetLargestPossibleRegion());
1368     DMapIteratorType iter2(dmap, dmap->GetLargestPossibleRegion());
1369     
1370     iter1.GoToBegin();
1371     iter2.GoToBegin();
1372     double dmin = 100000.0;
1373     typename ImageType::IndexType indexmin;
1374     indexmin.Fill(0);
1375     while (!iter1.IsAtEnd()) {
1376       if (iter1.Get() != BG) {
1377         double d = iter2.Get();
1378         if (d<dmin) {
1379           indexmin = iter1.GetIndex();
1380           dmin = d;
1381         }
1382       }
1383       ++iter1;
1384       ++iter2;
1385     }
1386     
1387     // Convert in Point
1388     typename ImageType::PointType p;
1389     input->TransformIndexToPhysicalPoint(indexmin, p);
1390     return p;
1391   }
1392   //--------------------------------------------------------------------
1393      
1394
1395   //--------------------------------------------------------------------
1396   template<class ImageType>
1397   typename ImageType::Pointer
1398   RemoveNegativeIndexFromRegion(ImageType * input) {
1399     typedef itk::ChangeInformationImageFilter< ImageType > InfoFilterType; 
1400     typename InfoFilterType::Pointer indexChangeFilter = InfoFilterType::New(); 
1401     indexChangeFilter->ChangeRegionOn(); 
1402     // The next line is commented because not exist in itk 3
1403     // typename InfoFilterType::OutputImageOffsetValueType indexShift[3];
1404     long indexShift[3];
1405     typename ImageType::IndexType index = input->GetLargestPossibleRegion().GetIndex();
1406     for(uint i=0;i<ImageType::ImageDimension; i++)
1407       indexShift[i] = (index[i]<0 ? -index[i]:0);
1408     typename ImageType::PointType origin;
1409     for(uint i=0;i<ImageType::ImageDimension; i++)
1410     origin[i] = input->GetOrigin()[i] - indexShift[i]*input->GetSpacing()[i];
1411     indexChangeFilter->SetOutputOffset( indexShift ); 
1412     indexChangeFilter->SetInput(input); 
1413     indexChangeFilter->SetOutputOrigin(origin);
1414     indexChangeFilter->ChangeOriginOn();
1415     indexChangeFilter->Update();
1416     return indexChangeFilter->GetOutput();
1417   }
1418   //--------------------------------------------------------------------
1419
1420
1421 } // end of namespace
1422