]> Creatis software - clitk.git/blob - itk/clitkSegmentationUtils.txx
Oups, wrong exception message
[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
38 namespace clitk {
39
40   //--------------------------------------------------------------------
41   template<class ImageType>
42   void ComputeBBFromImageRegion(const ImageType * image, 
43                                 typename ImageType::RegionType region,
44                                 typename itk::BoundingBox<unsigned long, 
45                                 ImageType::ImageDimension>::Pointer bb) {
46     typedef typename ImageType::IndexType IndexType;
47     IndexType firstIndex;
48     IndexType lastIndex;
49     for(unsigned int i=0; i<image->GetImageDimension(); i++) {
50       firstIndex[i] = region.GetIndex()[i];
51       lastIndex[i] = firstIndex[i]+region.GetSize()[i];
52     }
53
54     typedef itk::BoundingBox<unsigned long, 
55                              ImageType::ImageDimension> BBType;
56     typedef typename BBType::PointType PointType;
57     PointType lastPoint;
58     PointType firstPoint;
59     image->TransformIndexToPhysicalPoint(firstIndex, firstPoint);
60     image->TransformIndexToPhysicalPoint(lastIndex, lastPoint);
61
62     bb->SetMaximum(lastPoint);
63     bb->SetMinimum(firstPoint);
64   }
65   //--------------------------------------------------------------------
66
67
68   //--------------------------------------------------------------------
69   template<int Dimension>
70   void ComputeBBIntersection(typename itk::BoundingBox<unsigned long, Dimension>::Pointer bbo, 
71                              typename itk::BoundingBox<unsigned long, Dimension>::Pointer bbi1, 
72                              typename itk::BoundingBox<unsigned long, Dimension>::Pointer bbi2) {
73
74     typedef itk::BoundingBox<unsigned long, Dimension> BBType;
75     typedef typename BBType::PointType PointType;
76     PointType lastPoint;
77     PointType firstPoint;
78
79     for(unsigned int i=0; i<Dimension; i++) {
80       firstPoint[i] = std::max(bbi1->GetMinimum()[i], 
81                                bbi2->GetMinimum()[i]);
82       lastPoint[i] = std::min(bbi1->GetMaximum()[i], 
83                               bbi2->GetMaximum()[i]);
84     }
85
86     bbo->SetMaximum(lastPoint);
87     bbo->SetMinimum(firstPoint);
88   }
89   //--------------------------------------------------------------------
90
91
92   //--------------------------------------------------------------------
93   template<class ImageType>
94   void ComputeRegionFromBB(const ImageType * image, 
95                            const typename itk::BoundingBox<unsigned long, 
96                                                            ImageType::ImageDimension>::Pointer bb, 
97                            typename ImageType::RegionType & region) {
98     // Types
99     typedef typename ImageType::IndexType  IndexType;
100     typedef typename ImageType::PointType  PointType;
101     typedef typename ImageType::RegionType RegionType;
102     typedef typename ImageType::SizeType   SizeType;
103
104     // Region starting point
105     IndexType regionStart;
106     PointType start = bb->GetMinimum();
107     image->TransformPhysicalPointToIndex(start, regionStart);
108     
109     // Region size
110     SizeType regionSize;
111     PointType maxs = bb->GetMaximum();
112     PointType mins = bb->GetMinimum();
113     for(unsigned int i=0; i<ImageType::ImageDimension; i++) {
114       regionSize[i] = lrint((maxs[i] - mins[i])/image->GetSpacing()[i]);
115     }
116    
117     // Create region
118     region.SetIndex(regionStart);
119     region.SetSize(regionSize);
120   }
121   //--------------------------------------------------------------------
122
123   //--------------------------------------------------------------------
124   template<class ImageType, class TMaskImageType>
125   typename ImageType::Pointer
126   SetBackground(const ImageType * input, 
127                 const TMaskImageType * mask, 
128                 typename TMaskImageType::PixelType maskBG,
129                 typename ImageType::PixelType outValue, 
130                 bool inPlace) {
131     typedef SetBackgroundImageFilter<ImageType, TMaskImageType, ImageType> 
132       SetBackgroundImageFilterType;
133     typename SetBackgroundImageFilterType::Pointer setBackgroundFilter 
134       = SetBackgroundImageFilterType::New();
135     //  if (inPlace) setBackgroundFilter->ReleaseDataFlagOn(); // No seg fault
136     setBackgroundFilter->SetInPlace(inPlace); // This is important to keep memory low
137     setBackgroundFilter->SetInput(input);
138     setBackgroundFilter->SetInput2(mask);
139     setBackgroundFilter->SetMaskValue(maskBG);
140     setBackgroundFilter->SetOutsideValue(outValue);
141     setBackgroundFilter->Update();
142     return setBackgroundFilter->GetOutput();
143   }
144   //--------------------------------------------------------------------
145
146
147   //--------------------------------------------------------------------
148   template<class ImageType>
149   int GetNumberOfConnectedComponentLabels(const ImageType * input, 
150                                           typename ImageType::PixelType BG, 
151                                           bool isFullyConnected) {
152     // Connected Component label 
153     typedef itk::ConnectedComponentImageFilter<ImageType, ImageType> ConnectFilterType;
154     typename ConnectFilterType::Pointer connectFilter = ConnectFilterType::New();
155     connectFilter->SetInput(input);
156     connectFilter->SetBackgroundValue(BG);
157     connectFilter->SetFullyConnected(isFullyConnected);
158     connectFilter->Update();
159   
160     // Return result
161     return connectFilter->GetObjectCount();
162   }
163   //--------------------------------------------------------------------
164
165   //--------------------------------------------------------------------
166   /*
167     Warning : in this cas, we consider outputType like inputType, not
168     InternalImageType. Be sure it fits.
169   */
170   template<class ImageType>
171   typename ImageType::Pointer
172   Labelize(const ImageType * input, 
173            typename ImageType::PixelType BG, 
174            bool isFullyConnected, 
175            int minimalComponentSize) {
176     // InternalImageType for storing large number of component
177     typedef itk::Image<int, ImageType::ImageDimension> InternalImageType;
178   
179     // Connected Component label 
180     typedef itk::ConnectedComponentImageFilter<ImageType, InternalImageType> ConnectFilterType;
181     typename ConnectFilterType::Pointer connectFilter = ConnectFilterType::New();
182     //  connectFilter->ReleaseDataFlagOn(); 
183     connectFilter->SetInput(input);
184     connectFilter->SetBackgroundValue(BG);
185     connectFilter->SetFullyConnected(isFullyConnected);
186   
187     // Sort by size and remove too small area.
188     typedef itk::RelabelComponentImageFilter<InternalImageType, ImageType> RelabelFilterType;
189     typename RelabelFilterType::Pointer relabelFilter = RelabelFilterType::New();
190     //  relabelFilter->ReleaseDataFlagOn(); // if yes, fail when ExplosionControlledThresholdConnectedImageFilter ???
191     relabelFilter->SetInput(connectFilter->GetOutput());
192     relabelFilter->SetMinimumObjectSize(minimalComponentSize);
193     relabelFilter->Update();
194
195     // Return result
196     typename ImageType::Pointer output = relabelFilter->GetOutput();
197     return output;
198   }
199   //--------------------------------------------------------------------
200
201
202   //--------------------------------------------------------------------
203   /*
204     Warning : in this cas, we consider outputType like inputType, not
205     InternalImageType. Be sure it fits.
206   */
207   template<class ImageType>
208   typename ImageType::Pointer
209   LabelizeAndCountNumberOfObjects(const ImageType * input, 
210                                   typename ImageType::PixelType BG, 
211                                   bool isFullyConnected, 
212                                   int minimalComponentSize, 
213                                   int & nb) {
214     // InternalImageType for storing large number of component
215     typedef itk::Image<int, ImageType::ImageDimension> InternalImageType;
216   
217     // Connected Component label 
218     typedef itk::ConnectedComponentImageFilter<ImageType, InternalImageType> ConnectFilterType;
219     typename ConnectFilterType::Pointer connectFilter = ConnectFilterType::New();
220     //  connectFilter->ReleaseDataFlagOn(); 
221     connectFilter->SetInput(input);
222     connectFilter->SetBackgroundValue(BG);
223     connectFilter->SetFullyConnected(isFullyConnected);
224   
225     // Sort by size and remove too small area.
226     typedef itk::RelabelComponentImageFilter<InternalImageType, ImageType> RelabelFilterType;
227     typename RelabelFilterType::Pointer relabelFilter = RelabelFilterType::New();
228     //  relabelFilter->ReleaseDataFlagOn(); // if yes, fail when ExplosionControlledThresholdConnectedImageFilter ???
229     relabelFilter->SetInput(connectFilter->GetOutput());
230     relabelFilter->SetMinimumObjectSize(minimalComponentSize);
231     relabelFilter->Update();
232
233     nb = relabelFilter->GetNumberOfObjects();
234     // DD(relabelFilter->GetOriginalNumberOfObjects());
235     // DD(relabelFilter->GetSizeOfObjectsInPhysicalUnits()[0]);
236
237     // Return result
238     typename ImageType::Pointer output = relabelFilter->GetOutput();
239     return output;
240   }
241   //--------------------------------------------------------------------
242
243
244
245   //--------------------------------------------------------------------
246   template<class ImageType>
247   typename ImageType::Pointer
248   RemoveLabels(const ImageType * input, 
249                typename ImageType::PixelType BG,
250                std::vector<typename ImageType::PixelType> & labelsToRemove) {
251     assert(labelsToRemove.size() != 0);
252     typename ImageType::Pointer working_image;// = input;
253     for (unsigned int i=0; i <labelsToRemove.size(); i++) {
254       typedef SetBackgroundImageFilter<ImageType, ImageType> SetBackgroundImageFilterType;
255       typename SetBackgroundImageFilterType::Pointer setBackgroundFilter = SetBackgroundImageFilterType::New();
256       setBackgroundFilter->SetInput(input);
257       setBackgroundFilter->SetInput2(input);
258       setBackgroundFilter->SetMaskValue(labelsToRemove[i]);
259       setBackgroundFilter->SetOutsideValue(BG);
260       setBackgroundFilter->Update();
261       working_image = setBackgroundFilter->GetOutput();
262     }
263     return working_image;
264   }
265   //--------------------------------------------------------------------
266
267
268   //--------------------------------------------------------------------
269   template<class ImageType>
270   typename ImageType::Pointer
271   KeepLabels(const ImageType * input, 
272              typename ImageType::PixelType BG, 
273              typename ImageType::PixelType FG, 
274              typename ImageType::PixelType firstKeep, 
275              typename ImageType::PixelType lastKeep, 
276              bool useLastKeep) {
277     typedef itk::BinaryThresholdImageFilter<ImageType, ImageType> BinarizeFilterType; 
278     typename BinarizeFilterType::Pointer binarizeFilter = BinarizeFilterType::New();
279     binarizeFilter->SetInput(input);
280     binarizeFilter->SetLowerThreshold(firstKeep);
281     if (useLastKeep) binarizeFilter->SetUpperThreshold(lastKeep);
282     binarizeFilter->SetInsideValue(FG);
283     binarizeFilter->SetOutsideValue(BG);
284     binarizeFilter->Update();
285     return binarizeFilter->GetOutput();
286   }
287   //--------------------------------------------------------------------
288
289
290   //--------------------------------------------------------------------
291   template<class ImageType>
292   typename ImageType::Pointer
293   LabelizeAndSelectLabels(const ImageType * input,
294                           typename ImageType::PixelType BG, 
295                           typename ImageType::PixelType FG, 
296                           bool isFullyConnected,
297                           int minimalComponentSize,
298                           LabelizeParameters<typename ImageType::PixelType> * param)
299   {
300     typename ImageType::Pointer working_image;
301     working_image = Labelize<ImageType>(input, BG, isFullyConnected, minimalComponentSize);
302     if (param->GetLabelsToRemove().size() != 0)
303       working_image = RemoveLabels<ImageType>(working_image, BG, param->GetLabelsToRemove());
304     working_image = KeepLabels<ImageType>(working_image, 
305                                           BG, FG, 
306                                           param->GetFirstKeep(), 
307                                           param->GetLastKeep(), 
308                                           param->GetUseLastKeep());
309     return working_image;
310   }
311   //--------------------------------------------------------------------
312
313
314   //--------------------------------------------------------------------
315   template<class ImageType>
316   typename ImageType::Pointer
317   ResizeImageLike(const ImageType * input,                       
318                   const itk::ImageBase<ImageType::ImageDimension> * like, 
319                   typename ImageType::PixelType backgroundValue) 
320   {
321     typedef CropLikeImageFilter<ImageType> CropFilterType;
322     typename CropFilterType::Pointer cropFilter = CropFilterType::New();
323     cropFilter->SetInput(input);
324     cropFilter->SetCropLikeImage(like);
325     cropFilter->SetBackgroundValue(backgroundValue);
326     cropFilter->Update();
327     return cropFilter->GetOutput();  
328   }
329   //--------------------------------------------------------------------
330
331
332   //--------------------------------------------------------------------
333   template<class MaskImageType>
334   typename MaskImageType::Pointer
335   SliceBySliceRelativePosition(const MaskImageType * input,
336                                const MaskImageType * object,
337                                int direction, 
338                                double threshold, 
339                                std::string orientation, 
340                                bool uniqueConnectedComponent, 
341                                double spacing, 
342                                bool autocropFlag, 
343                                bool singleObjectCCL) 
344   {
345     typedef SliceBySliceRelativePositionFilter<MaskImageType> SliceRelPosFilterType;
346     typename SliceRelPosFilterType::Pointer sliceRelPosFilter = SliceRelPosFilterType::New();
347     sliceRelPosFilter->VerboseStepFlagOff();
348     sliceRelPosFilter->WriteStepFlagOff();
349     sliceRelPosFilter->SetInput(input);
350     sliceRelPosFilter->SetInputObject(object);
351     sliceRelPosFilter->SetDirection(direction);
352     sliceRelPosFilter->SetFuzzyThreshold(threshold);
353     sliceRelPosFilter->AddOrientationTypeString(orientation);
354     sliceRelPosFilter->SetIntermediateSpacingFlag((spacing != -1));
355     sliceRelPosFilter->SetIntermediateSpacing(spacing);
356     sliceRelPosFilter->SetUniqueConnectedComponentBySlice(uniqueConnectedComponent);
357     sliceRelPosFilter->SetUseASingleObjectConnectedComponentBySliceFlag(singleObjectCCL);
358     //    sliceRelPosFilter->SetInverseOrientationFlag(inverseflag); 
359     sliceRelPosFilter->SetAutoCropFlag(autocropFlag); 
360     sliceRelPosFilter->IgnoreEmptySliceObjectFlagOn();
361     sliceRelPosFilter->Update();
362     return sliceRelPosFilter->GetOutput();
363   }
364   //--------------------------------------------------------------------
365
366   //--------------------------------------------------------------------
367   template<class ImageType>
368   bool
369   FindExtremaPointInAGivenDirection(const ImageType * input, 
370                                     typename ImageType::PixelType bg, 
371                                     int direction, bool opposite, 
372                                     typename ImageType::PointType & point)
373   {
374     typename ImageType::PointType dummy;
375     return FindExtremaPointInAGivenDirection(input, bg, direction, 
376                                              opposite, dummy, 0, point);
377   }
378   //--------------------------------------------------------------------
379
380   //--------------------------------------------------------------------
381   template<class ImageType>
382   bool
383   FindExtremaPointInAGivenDirection(const ImageType * input, 
384                                     typename ImageType::PixelType bg, 
385                                     int direction, bool opposite, 
386                                     typename ImageType::PointType refpoint,
387                                     double distanceMax, 
388                                     typename ImageType::PointType & point)
389   {
390     /*
391       loop over input pixels, store the index in the fg that is max
392       according to the given direction. 
393     */    
394     typedef itk::ImageRegionConstIteratorWithIndex<ImageType> IteratorType;
395     IteratorType iter(input, input->GetLargestPossibleRegion());
396     iter.GoToBegin();
397     typename ImageType::IndexType max = input->GetLargestPossibleRegion().GetIndex();
398     if (opposite) max = max+input->GetLargestPossibleRegion().GetSize();
399     bool found=false;
400     while (!iter.IsAtEnd()) {
401       if (iter.Get() != bg) {
402         bool test = iter.GetIndex()[direction] >  max[direction];
403         if (opposite) test = !test;
404         if (test) {
405           typename ImageType::PointType p;
406           input->TransformIndexToPhysicalPoint(iter.GetIndex(), p);
407           if ((distanceMax==0) || (p.EuclideanDistanceTo(refpoint) < distanceMax)) {
408             max = iter.GetIndex();
409             found = true;
410           }
411         }
412       }
413       ++iter;
414     }
415     if (!found) return false;
416     input->TransformIndexToPhysicalPoint(max, point);
417     return true;
418   }
419   //--------------------------------------------------------------------
420
421
422   //--------------------------------------------------------------------
423   template<class ImageType>
424   typename ImageType::Pointer
425   CropImageRemoveGreaterThan(const ImageType * image, 
426                  int dim, double min, bool autoCrop,
427                  typename ImageType::PixelType BG) 
428   {
429     return CropImageAlongOneAxis<ImageType>(image, dim, 
430                                             image->GetOrigin()[dim], 
431                                             min,
432                                             autoCrop, BG);
433   }
434   //--------------------------------------------------------------------
435
436
437   //--------------------------------------------------------------------
438   template<class ImageType>
439   typename ImageType::Pointer
440   CropImageRemoveLowerThan(const ImageType * image, 
441                  int dim, double max, bool autoCrop,
442                  typename ImageType::PixelType BG) 
443   {
444     typename ImageType::PointType p;
445     image->TransformIndexToPhysicalPoint(image->GetLargestPossibleRegion().GetIndex()+
446                                          image->GetLargestPossibleRegion().GetSize(), p);
447     return CropImageAlongOneAxis<ImageType>(image, dim, max, p[dim], autoCrop, BG);
448   }
449   //--------------------------------------------------------------------
450
451
452   //--------------------------------------------------------------------
453   template<class ImageType>
454   typename ImageType::Pointer
455   CropImageAlongOneAxis(const ImageType * image, 
456                         int dim, double min, double max, 
457                         bool autoCrop, typename ImageType::PixelType BG) 
458   {
459     // Compute region size
460     typename ImageType::RegionType region;
461     typename ImageType::SizeType size = image->GetLargestPossibleRegion().GetSize();
462     typename ImageType::PointType p = image->GetOrigin();
463     p[dim] = min;
464     typename ImageType::IndexType start;
465     image->TransformPhysicalPointToIndex(p, start);
466     p[dim] = max;
467     typename ImageType::IndexType end;
468     image->TransformPhysicalPointToIndex(p, end);
469     size[dim] = fabs(end[dim]-start[dim]);
470     region.SetIndex(start);
471     region.SetSize(size);
472   
473     // Perform Crop
474     typedef itk::RegionOfInterestImageFilter<ImageType, ImageType> CropFilterType;
475     typename CropFilterType::Pointer cropFilter = CropFilterType::New();
476     cropFilter->SetInput(image);
477     cropFilter->SetRegionOfInterest(region);
478     cropFilter->Update();
479     typename ImageType::Pointer result = cropFilter->GetOutput();
480   
481     // Auto Crop
482     if (autoCrop) {
483       result = AutoCrop<ImageType>(result, BG);
484     }
485     return result;
486   }
487   //--------------------------------------------------------------------
488
489
490   //--------------------------------------------------------------------
491   template<class ImageType>
492   void
493   ComputeCentroids(const ImageType * image, 
494                    typename ImageType::PixelType BG, 
495                    std::vector<typename ImageType::PointType> & centroids) 
496   {
497     typedef long LabelType;
498     static const unsigned int Dim = ImageType::ImageDimension;
499     typedef itk::ShapeLabelObject< LabelType, Dim > LabelObjectType;
500     typedef itk::LabelMap< LabelObjectType > LabelMapType;
501     typedef itk::LabelImageToLabelMapFilter<ImageType, LabelMapType> ImageToMapFilterType;
502     typename ImageToMapFilterType::Pointer imageToLabelFilter = ImageToMapFilterType::New(); 
503     typedef itk::ShapeLabelMapFilter<LabelMapType, ImageType> ShapeFilterType; 
504     typename ShapeFilterType::Pointer statFilter = ShapeFilterType::New();
505     imageToLabelFilter->SetBackgroundValue(BG);
506     imageToLabelFilter->SetInput(image);
507     statFilter->SetInput(imageToLabelFilter->GetOutput());
508     statFilter->Update();
509     typename LabelMapType::Pointer labelMap = statFilter->GetOutput();
510
511     centroids.clear();
512     typename ImageType::PointType dummy;
513     centroids.push_back(dummy); // label 0 -> no centroid, use dummy point for BG 
514     //DS FIXME (not useful ! to change ..)
515     for(uint i=0; i<labelMap->GetNumberOfLabelObjects(); i++) {
516       int label = labelMap->GetLabels()[i];
517       centroids.push_back(labelMap->GetLabelObject(label)->GetCentroid());
518     } 
519   }
520   //--------------------------------------------------------------------
521
522
523   //--------------------------------------------------------------------
524   template<class ImageType, class LabelType>
525   typename itk::LabelMap< itk::ShapeLabelObject<LabelType, ImageType::ImageDimension> >::Pointer
526   ComputeLabelMap(const ImageType * image, 
527                   typename ImageType::PixelType BG, 
528                   bool computePerimeterFlag) 
529   {
530     static const unsigned int Dim = ImageType::ImageDimension;
531     typedef itk::ShapeLabelObject< LabelType, Dim > LabelObjectType;
532     typedef itk::LabelMap< LabelObjectType > LabelMapType;
533     typedef itk::LabelImageToLabelMapFilter<ImageType, LabelMapType> ImageToMapFilterType;
534     typename ImageToMapFilterType::Pointer imageToLabelFilter = ImageToMapFilterType::New(); 
535     typedef itk::ShapeLabelMapFilter<LabelMapType, ImageType> ShapeFilterType; 
536     typename ShapeFilterType::Pointer statFilter = ShapeFilterType::New();
537     imageToLabelFilter->SetBackgroundValue(BG);
538     imageToLabelFilter->SetInput(image);
539     statFilter->SetInput(imageToLabelFilter->GetOutput());
540     statFilter->SetComputePerimeter(computePerimeterFlag);
541     statFilter->Update();
542     return statFilter->GetOutput();
543   }
544   //--------------------------------------------------------------------
545
546
547   //--------------------------------------------------------------------
548   template<class ImageType>
549   void
550   ComputeCentroids2(const ImageType * image, 
551                    typename ImageType::PixelType BG, 
552                    std::vector<typename ImageType::PointType> & centroids) 
553   {
554     typedef long LabelType;
555     static const unsigned int Dim = ImageType::ImageDimension;
556     typedef itk::ShapeLabelObject< LabelType, Dim > LabelObjectType;
557     typedef itk::LabelMap< LabelObjectType > LabelMapType;
558     typedef itk::LabelImageToLabelMapFilter<ImageType, LabelMapType> ImageToMapFilterType;
559     typename ImageToMapFilterType::Pointer imageToLabelFilter = ImageToMapFilterType::New(); 
560     typedef itk::ShapeLabelMapFilter<LabelMapType, ImageType> ShapeFilterType; 
561     typename ShapeFilterType::Pointer statFilter = ShapeFilterType::New();
562     imageToLabelFilter->SetBackgroundValue(BG);
563     imageToLabelFilter->SetInput(image);
564     statFilter->SetInput(imageToLabelFilter->GetOutput());
565     statFilter->Update();
566     typename LabelMapType::Pointer labelMap = statFilter->GetOutput();
567
568     centroids.clear();
569     typename ImageType::PointType dummy;
570     centroids.push_back(dummy); // label 0 -> no centroid, use dummy point
571     for(uint i=1; i<labelMap->GetNumberOfLabelObjects()+1; i++) {
572       centroids.push_back(labelMap->GetLabelObject(i)->GetCentroid());
573     } 
574     
575     for(uint i=1; i<labelMap->GetNumberOfLabelObjects()+1; i++) {
576       DD(labelMap->GetLabelObject(i)->GetBinaryPrincipalAxes());
577       DD(labelMap->GetLabelObject(i)->GetBinaryFlatness());
578       DD(labelMap->GetLabelObject(i)->GetRoundness ());      
579
580       // search for the point on the boundary alog PA
581
582     }
583
584   }
585   //--------------------------------------------------------------------
586
587
588   //--------------------------------------------------------------------
589   template<class ImageType>
590   void
591   ExtractSlices(const ImageType * image, int direction, 
592                 std::vector<typename itk::Image<typename ImageType::PixelType, 
593                                                 ImageType::ImageDimension-1>::Pointer > & slices) 
594   {
595     typedef ExtractSliceFilter<ImageType> ExtractSliceFilterType;
596     typedef typename ExtractSliceFilterType::SliceType SliceType;
597     typename ExtractSliceFilterType::Pointer 
598       extractSliceFilter = ExtractSliceFilterType::New();
599     extractSliceFilter->SetInput(image);
600     extractSliceFilter->SetDirection(direction);
601     extractSliceFilter->Update();
602     extractSliceFilter->GetOutputSlices(slices);
603   }
604   //--------------------------------------------------------------------
605
606
607   //--------------------------------------------------------------------
608   template<class ImageType>
609   void
610   PointsUtils<ImageType>::Convert2DTo3D(const PointType2D & p2D, 
611                                         const ImageType * image, 
612                                         const int slice, 
613                                         PointType3D & p3D)  
614   {
615     IndexType3D index3D;
616     index3D[0] = index3D[1] = 0;
617     index3D[2] = image->GetLargestPossibleRegion().GetIndex()[2]+slice;
618     image->TransformIndexToPhysicalPoint(index3D, p3D);
619     p3D[0] = p2D[0]; 
620     p3D[1] = p2D[1];
621     //  p3D[2] = p[2];//(image->GetLargestPossibleRegion().GetIndex()[2]+slice)*image->GetSpacing()[2] 
622     //    + image->GetOrigin()[2];
623   }
624   //--------------------------------------------------------------------
625
626
627   //--------------------------------------------------------------------
628   template<class ImageType>
629   void 
630   PointsUtils<ImageType>::Convert2DMapTo3DList(const MapPoint2DType & map, 
631                                             const ImageType * image, 
632                                             VectorPoint3DType & list)
633   {
634     typename MapPoint2DType::const_iterator iter = map.begin();
635     while (iter != map.end()) {
636       PointType3D p;
637       Convert2DTo3D(iter->second, image, iter->first, p);
638       list.push_back(p);
639       ++iter;
640     }
641   }
642   //--------------------------------------------------------------------
643
644
645   //--------------------------------------------------------------------
646   template<class ImageType>
647   void 
648   PointsUtils<ImageType>::Convert2DListTo3DList(const VectorPoint2DType & p2D, 
649                                                 int slice,
650                                                 const ImageType * image, 
651                                                 VectorPoint3DType & list) 
652   {
653     for(uint i=0; i<p2D.size(); i++) {
654       PointType3D p;
655       Convert2DTo3D(p2D[i], image, slice, p);
656       list.push_back(p);
657     }
658   }
659   //--------------------------------------------------------------------
660
661
662   //--------------------------------------------------------------------
663   template<class ImageType>
664   void 
665   WriteListOfLandmarks(std::vector<typename ImageType::PointType> points, 
666                        std::string filename)
667   {
668     std::ofstream os; 
669     openFileForWriting(os, filename); 
670     os << "LANDMARKS1" << std::endl;  
671     for(uint i=0; i<points.size(); i++) {
672       const typename ImageType::PointType & p = points[i];
673       // Write it in the file
674       os << i << " " << p[0] << " " << p[1] << " " << p[2] << " 0 0 " << std::endl;
675     }
676     os.close();
677   }
678   //--------------------------------------------------------------------
679
680
681   //--------------------------------------------------------------------
682   template<class ImageType>
683   typename ImageType::Pointer 
684   Dilate(const ImageType * image, double radiusInMM,               
685          typename ImageType::PixelType BG,
686          typename ImageType::PixelType FG,  
687          bool extendSupport)
688   {
689     typename ImageType::SizeType r;
690     for(uint i=0; i<ImageType::ImageDimension; i++) 
691       r[i] = (uint)lrint(radiusInMM/image->GetSpacing()[i]);
692     return Dilate<ImageType>(image, r, BG, FG, extendSupport);
693   }
694   //--------------------------------------------------------------------
695
696
697   //--------------------------------------------------------------------
698   template<class ImageType>
699   typename ImageType::Pointer 
700   Dilate(const ImageType * image, typename ImageType::PointType radiusInMM, 
701          typename ImageType::PixelType BG, 
702          typename ImageType::PixelType FG, 
703          bool extendSupport)
704   {
705     typename ImageType::SizeType r;
706     for(uint i=0; i<ImageType::ImageDimension; i++) 
707       r[i] = (uint)lrint(radiusInMM[i]/image->GetSpacing()[i]);
708     return Dilate<ImageType>(image, r, BG, FG, extendSupport);
709   }
710   //--------------------------------------------------------------------
711
712
713   //--------------------------------------------------------------------
714   template<class ImageType>
715   typename ImageType::Pointer 
716   Dilate(const ImageType * image, typename ImageType::SizeType radius, 
717          typename ImageType::PixelType BG, 
718          typename ImageType::PixelType FG, 
719          bool extendSupport)
720   {
721     // Create kernel for dilatation
722     typedef itk::BinaryBallStructuringElement<typename ImageType::PixelType, 
723                                               ImageType::ImageDimension> KernelType;
724     KernelType structuringElement;
725     structuringElement.SetRadius(radius);
726     structuringElement.CreateStructuringElement();
727
728     typename ImageType::Pointer output;
729     if (extendSupport) {
730       typedef itk::ConstantPadImageFilter<ImageType, ImageType> PadFilterType;
731       typename PadFilterType::Pointer padFilter = PadFilterType::New();
732       padFilter->SetInput(image);
733       typename ImageType::SizeType lower;
734       typename ImageType::SizeType upper;
735       for(uint i=0; i<3; i++) {
736         lower[i] = upper[i] = 2*(radius[i]+1);
737       }
738       padFilter->SetPadLowerBound(lower);
739       padFilter->SetPadUpperBound(upper);
740       padFilter->Update();
741       output = padFilter->GetOutput();
742     }
743
744     // Dilate  filter
745     typedef itk::BinaryDilateImageFilter<ImageType, ImageType , KernelType> DilateFilterType;
746     typename DilateFilterType::Pointer dilateFilter = DilateFilterType::New();
747     dilateFilter->SetBackgroundValue(BG);
748     dilateFilter->SetForegroundValue(FG);
749     dilateFilter->SetBoundaryToForeground(false);
750     dilateFilter->SetKernel(structuringElement);
751     dilateFilter->SetInput(output);
752     dilateFilter->Update();
753     return dilateFilter->GetOutput();
754   }
755   //--------------------------------------------------------------------
756
757
758   //--------------------------------------------------------------------
759   template<class ImageType>
760   typename ImageType::Pointer 
761   Opening(const ImageType * image, typename ImageType::SizeType radius,
762          typename ImageType::PixelType BG,
763          typename ImageType::PixelType FG)
764   {
765     // Kernel 
766     typedef itk::BinaryBallStructuringElement<typename ImageType::PixelType, 
767                                               ImageType::ImageDimension> KernelType;    
768     KernelType structuringElement;
769     structuringElement.SetRadius(radius);
770     structuringElement.CreateStructuringElement();
771     
772     // Filter
773     typedef itk::BinaryMorphologicalOpeningImageFilter<ImageType, ImageType , KernelType> OpeningFilterType;
774     typename OpeningFilterType::Pointer open = OpeningFilterType::New();
775     open->SetInput(image);
776     open->SetBackgroundValue(BG);
777     open->SetForegroundValue(FG);
778     open->SetKernel(structuringElement);
779     open->Update();
780     return open->GetOutput();
781   }
782   //--------------------------------------------------------------------
783
784
785
786   //--------------------------------------------------------------------
787   template<class ValueType, class VectorType>
788   void ConvertOption(std::string optionName, uint given, 
789                      ValueType * values, VectorType & p, 
790                      uint dim, bool required) 
791   {
792     if (required && (given == 0)) {
793       clitkExceptionMacro("The option --" << optionName << " must be set and have 1 or " 
794                           << dim << " values.");
795     }
796     if (given == 1) {
797       for(uint i=0; i<dim; i++) p[i] = values[0];
798       return;
799     }
800     if (given == dim) {
801       for(uint i=0; i<dim; i++) p[i] = values[i];
802       return;
803     }
804     if (given == 0) return;
805     clitkExceptionMacro("The option --" << optionName << " must have 1 or " 
806                         << dim << " values.");
807   }
808   //--------------------------------------------------------------------
809
810
811   //--------------------------------------------------------------------
812   /*
813     http://www.gamedev.net/community/forums/topic.asp?topic_id=542870
814     Assuming the points are (Ax,Ay) (Bx,By) and (Cx,Cy), you need to compute:
815     (Bx - Ax) * (Cy - Ay) - (By - Ay) * (Cx - Ax)
816     This will equal zero if the point C is on the line formed by
817     points A and B, and will have a different sign depending on the
818     side. Which side this is depends on the orientation of your (x,y)
819     coordinates, but you can plug test values for A,B and C into this
820     formula to determine whether negative values are to the left or to
821     the right.
822     => to accelerate, start with formula, when change sign -> stop and fill
823
824     offsetToKeep = is used to determine which side of the line we
825     keep. The point along the mainDirection but 'offsetToKeep' mm away
826     is kept.
827   
828   */
829   template<class ImageType>
830   void 
831   SliceBySliceSetBackgroundFromLineSeparation(ImageType * input, 
832                                               std::vector<typename ImageType::PointType> & lA, 
833                                               std::vector<typename ImageType::PointType> & lB, 
834                                               typename ImageType::PixelType BG, 
835                                               int mainDirection, 
836                                               double offsetToKeep)
837   {
838     typedef itk::ImageSliceIteratorWithIndex<ImageType> SliceIteratorType;
839     SliceIteratorType siter = SliceIteratorType(input, 
840                                                 input->GetLargestPossibleRegion());
841     siter.SetFirstDirection(0);
842     siter.SetSecondDirection(1);
843     siter.GoToBegin();
844     uint i=0;
845     typename ImageType::PointType A;
846     typename ImageType::PointType B;
847     typename ImageType::PointType C;
848     assert(lA.size() == lB.size());
849     while ((i<lA.size()) && (!siter.IsAtEnd())) {
850       // Check that the current slice correspond to the current point
851       input->TransformIndexToPhysicalPoint(siter.GetIndex(), C);
852       // DD(C);
853       // DD(i);
854       // DD(lA[i]);
855       if ((fabs(C[2] - lA[i][2]))>0.01) { // is !equal with a tolerance of 0.01 mm
856       }
857       else {
858         // Define A,B,C points
859         A = lA[i];
860         B = lB[i];
861         C = A;
862         // DD(A);
863         // DD(B);
864         // DD(C);
865       
866         // Check that the line is not a point (A=B)
867         bool p = (A[0] == B[0]) && (A[1] == B[1]);
868       
869         if (!p) {
870           C[mainDirection] += offsetToKeep; // I know I must keep this point
871           double s = (B[0] - A[0]) * (C[1] - A[1]) - (B[1] - A[1]) * (C[0] - A[0]);
872           bool isPositive = s<0;
873           while (!siter.IsAtEndOfSlice()) {
874             while (!siter.IsAtEndOfLine()) {
875               // Very slow, I know ... but image should be very small
876               input->TransformIndexToPhysicalPoint(siter.GetIndex(), C);
877               double s = (B[0] - A[0]) * (C[1] - A[1]) - (B[1] - A[1]) * (C[0] - A[0]);
878               if (s == 0) siter.Set(BG); // on the line, we decide to remove
879               if (isPositive) {
880                 if (s > 0) siter.Set(BG);
881               }
882               else {
883                 if (s < 0) siter.Set(BG); 
884               }
885               ++siter;
886             }
887             siter.NextLine();
888           } // end loop slice
889         }      
890
891         ++i;
892       } // End of current slice
893       siter.NextSlice();
894     }
895   }                                                   
896   //--------------------------------------------------------------------
897
898   //--------------------------------------------------------------------
899   template<class ImageType>
900   void 
901   AndNot(ImageType * input, 
902          const ImageType * object, 
903          typename ImageType::PixelType BG)
904   {
905     typename ImageType::Pointer o;
906     bool resized=false;
907     if (!clitk::HaveSameSizeAndSpacing<ImageType, ImageType>(input, object)) {
908       o = clitk::ResizeImageLike<ImageType>(object, input, BG);
909       resized = true;
910     }
911
912     typedef clitk::BooleanOperatorLabelImageFilter<ImageType> BoolFilterType;
913     typename BoolFilterType::Pointer boolFilter = BoolFilterType::New(); 
914     boolFilter->InPlaceOn();
915     boolFilter->SetInput1(input);
916     if (resized) boolFilter->SetInput2(o);  
917     else boolFilter->SetInput2(object);
918     boolFilter->SetBackgroundValue1(BG);
919     boolFilter->SetBackgroundValue2(BG);
920     boolFilter->SetOperationType(BoolFilterType::AndNot);
921     boolFilter->Update();
922   }
923   //--------------------------------------------------------------------
924
925
926   //--------------------------------------------------------------------
927   template<class ImageType>
928   typename ImageType::Pointer
929   Binarize(const ImageType * input, 
930            typename ImageType::PixelType lower, 
931            typename ImageType::PixelType upper, 
932            typename ImageType::PixelType BG,
933            typename ImageType::PixelType FG) 
934   {
935     typedef itk::BinaryThresholdImageFilter<ImageType, ImageType> BinaryThresholdFilterType;
936     typename BinaryThresholdFilterType::Pointer binarizeFilter = BinaryThresholdFilterType::New();
937     binarizeFilter->SetInput(input);
938     binarizeFilter->InPlaceOff();
939     binarizeFilter->SetLowerThreshold(lower);
940     binarizeFilter->SetUpperThreshold(upper);
941     binarizeFilter->SetInsideValue(FG);
942     binarizeFilter->SetOutsideValue(BG);
943     binarizeFilter->Update();
944     return binarizeFilter->GetOutput();
945   }
946   //--------------------------------------------------------------------
947
948
949   //--------------------------------------------------------------------
950   template<class ImageType>
951   void
952   GetMinMaxPointPosition(const ImageType * input, 
953                          typename ImageType::PointType & min,
954                          typename ImageType::PointType & max) 
955   {
956     typename ImageType::IndexType index = input->GetLargestPossibleRegion().GetIndex();
957     input->TransformIndexToPhysicalPoint(index, min);
958     index = index+input->GetLargestPossibleRegion().GetSize();
959     input->TransformIndexToPhysicalPoint(index, max);
960   }
961   //--------------------------------------------------------------------
962
963
964   //--------------------------------------------------------------------
965   template<class ImageType>
966   typename ImageType::PointType
967   FindExtremaPointInAGivenLine(const ImageType * input, 
968                                int dimension, 
969                                bool inverse, 
970                                typename ImageType::PointType p, 
971                                typename ImageType::PixelType BG, 
972                                double distanceMax) 
973   {
974     // Which direction ?  Increasing or decreasing.
975     int d=1;
976     if (inverse) d=-1;
977   
978     // Transform to pixel index
979     typename ImageType::IndexType index;
980     input->TransformPhysicalPointToIndex(p, index);
981
982     // Loop while inside the mask;
983     while (input->GetPixel(index) != BG) {
984       index[dimension] += d;
985     }
986
987     // Transform back to Physical Units
988     typename ImageType::PointType result;
989     input->TransformIndexToPhysicalPoint(index, result);
990
991     // Check that is is not too far away
992     double distance = p.EuclideanDistanceTo(result);
993     if (distance > distanceMax) {
994       result = p; // Get back to initial value
995     }
996
997     return result;
998   }
999   //--------------------------------------------------------------------
1000
1001
1002   //--------------------------------------------------------------------
1003   template<class PointType>
1004   bool
1005   IsOnTheSameLineSide(PointType C, PointType A, PointType B, PointType like) 
1006   {
1007     // Look at the position of point 'like' according to the AB line
1008     double s = (B[0] - A[0]) * (like[1] - A[1]) - (B[1] - A[1]) * (like[0] - A[0]);
1009     bool negative = s<0;
1010   
1011     // Look the C position
1012     s = (B[0] - A[0]) * (C[1] - A[1]) - (B[1] - A[1]) * (C[0] - A[0]);
1013
1014     if (negative && (s<=0)) return true;
1015     if (!negative && (s>=0)) return true;
1016     return false;
1017   }
1018   //--------------------------------------------------------------------
1019
1020
1021   //--------------------------------------------------------------------
1022   /* Consider an input object, for each slice, find the extrema
1023      position according to a given direction and build a line segment
1024      passing throught this point in a given direction.  Output is a
1025      vector of line (from point A to B), for each slice;
1026    */
1027   template<class ImageType>
1028   void 
1029   SliceBySliceBuildLineSegmentAccordingToExtremaPosition(const ImageType * input, 
1030                                                          typename ImageType::PixelType BG, 
1031                                                          int sliceDimension, 
1032                                                          int extremaDirection, 
1033                                                          bool extremaOppositeFlag, 
1034                                                          int lineDirection,
1035                                                          double margin,
1036                                                          std::vector<typename ImageType::PointType> & A, 
1037                                                          std::vector<typename ImageType::PointType> & B)
1038   {
1039     // Type of a slice
1040     typedef typename itk::Image<typename ImageType::PixelType, ImageType::ImageDimension-1> SliceType;
1041     
1042     // Build the list of slices
1043     std::vector<typename SliceType::Pointer> slices;
1044     clitk::ExtractSlices<ImageType>(input, sliceDimension, slices);
1045
1046     // Build the list of 2D points
1047     std::map<int, typename SliceType::PointType> position2D;
1048     for(uint i=0; i<slices.size(); i++) {
1049       typename SliceType::PointType p;
1050       bool found = 
1051         clitk::FindExtremaPointInAGivenDirection<SliceType>(slices[i], BG, 
1052                                                             extremaDirection, extremaOppositeFlag, p);
1053       if (found) {
1054         position2D[i] = p;
1055       }    
1056     }
1057     
1058     // Convert 2D points in slice into 3D points
1059     clitk::PointsUtils<ImageType>::Convert2DMapTo3DList(position2D, input, A);
1060     
1061     // Create additional point just right to the previous ones, on the
1062     // given lineDirection, in order to create a horizontal/vertical line.
1063     for(uint i=0; i<A.size(); i++) {
1064       typename ImageType::PointType p = A[i];
1065       p[lineDirection] += 10;
1066       B.push_back(p);
1067       // Margins ?
1068       A[i][1] += margin;
1069       B[i][1] += margin;
1070     }
1071
1072   }
1073   //--------------------------------------------------------------------
1074
1075
1076   //--------------------------------------------------------------------
1077   template<class ImageType>
1078   typename ImageType::Pointer
1079   SliceBySliceKeepMainCCL(const ImageType * input, 
1080                           typename ImageType::PixelType BG,
1081                           typename ImageType::PixelType FG)  {
1082     
1083     // Extract slices
1084     const int d = ImageType::ImageDimension-1;
1085     typedef typename itk::Image<typename ImageType::PixelType, d> SliceType;
1086     std::vector<typename SliceType::Pointer> slices;
1087     clitk::ExtractSlices<ImageType>(input, d, slices);
1088     
1089     // Labelize and keep the main one
1090     std::vector<typename SliceType::Pointer> o;
1091     for(uint i=0; i<slices.size(); i++) {
1092       o.push_back(clitk::Labelize<SliceType>(slices[i], BG, false, 1));
1093       o[i] = clitk::KeepLabels<SliceType>(o[i], BG, FG, 1, 1, true);
1094     }
1095     
1096     // Join slices
1097     typename ImageType::Pointer output;
1098     output = clitk::JoinSlices<ImageType>(o, input, d);
1099     return output;
1100   }
1101   //--------------------------------------------------------------------
1102
1103
1104   //--------------------------------------------------------------------
1105   template<class ImageType>
1106   typename ImageType::Pointer
1107   Clone(const ImageType * input) {
1108     typedef itk::ImageDuplicator<ImageType> DuplicatorType;
1109     typename DuplicatorType::Pointer duplicator = DuplicatorType::New();
1110     duplicator->SetInputImage(input);
1111     duplicator->Update();
1112     return duplicator->GetOutput();
1113   }
1114   //--------------------------------------------------------------------
1115
1116
1117 } // end of namespace
1118