1 /*=========================================================================
2 Program: vv http://www.creatis.insa-lyon.fr/rio/vv
5 - University of LYON http://www.universite-lyon.fr/
6 - Léon Bérard cancer center http://www.centreleonberard.fr
7 - CREATIS CNRS laboratory http://www.creatis.insa-lyon.fr
9 This software is distributed WITHOUT ANY WARRANTY; without even
10 the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
11 PURPOSE. See the copyright notices for more information.
13 It is distributed under dual licence
15 - BSD See included LICENSE.txt file
16 - CeCILL-B http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
17 ======================================================================-====*/
20 #include "clitkSetBackgroundImageFilter.h"
21 #include "clitkSliceBySliceRelativePositionFilter.h"
22 #include "clitkCropLikeImageFilter.h"
23 #include "clitkMemoryUsage.h"
26 #include <itkConnectedComponentImageFilter.h>
27 #include <itkRelabelComponentImageFilter.h>
28 #include <itkBinaryThresholdImageFilter.h>
29 #include <itkPasteImageFilter.h>
30 #include <itkStatisticsLabelMapFilter.h>
31 #include <itkBinaryBallStructuringElement.h>
32 #include <itkBinaryDilateImageFilter.h>
33 #include <itkConstantPadImageFilter.h>
34 #include <itkImageSliceIteratorWithIndex.h>
35 #include <itkBinaryMorphologicalOpeningImageFilter.h>
36 #include <itkImageDuplicator.h>
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;
49 for(unsigned int i=0; i<image->GetImageDimension(); i++) {
50 firstIndex[i] = region.GetIndex()[i];
51 lastIndex[i] = firstIndex[i]+region.GetSize()[i];
54 typedef itk::BoundingBox<unsigned long,
55 ImageType::ImageDimension> BBType;
56 typedef typename BBType::PointType PointType;
59 image->TransformIndexToPhysicalPoint(firstIndex, firstPoint);
60 image->TransformIndexToPhysicalPoint(lastIndex, lastPoint);
62 bb->SetMaximum(lastPoint);
63 bb->SetMinimum(firstPoint);
65 //--------------------------------------------------------------------
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) {
74 typedef itk::BoundingBox<unsigned long, Dimension> BBType;
75 typedef typename BBType::PointType PointType;
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]);
86 bbo->SetMaximum(lastPoint);
87 bbo->SetMinimum(firstPoint);
89 //--------------------------------------------------------------------
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) {
99 typedef typename ImageType::IndexType IndexType;
100 typedef typename ImageType::PointType PointType;
101 typedef typename ImageType::RegionType RegionType;
102 typedef typename ImageType::SizeType SizeType;
104 // Region starting point
105 IndexType regionStart;
106 PointType start = bb->GetMinimum();
107 image->TransformPhysicalPointToIndex(start, regionStart);
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]);
118 region.SetIndex(regionStart);
119 region.SetSize(regionSize);
121 //--------------------------------------------------------------------
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,
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();
144 //--------------------------------------------------------------------
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();
161 return connectFilter->GetObjectCount();
163 //--------------------------------------------------------------------
165 //--------------------------------------------------------------------
167 Warning : in this cas, we consider outputType like inputType, not
168 InternalImageType. Be sure it fits.
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;
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);
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();
196 typename ImageType::Pointer output = relabelFilter->GetOutput();
199 //--------------------------------------------------------------------
202 //--------------------------------------------------------------------
204 Warning : in this cas, we consider outputType like inputType, not
205 InternalImageType. Be sure it fits.
207 template<class ImageType>
208 typename ImageType::Pointer
209 LabelizeAndCountNumberOfObjects(const ImageType * input,
210 typename ImageType::PixelType BG,
211 bool isFullyConnected,
212 int minimalComponentSize,
214 // InternalImageType for storing large number of component
215 typedef itk::Image<int, ImageType::ImageDimension> InternalImageType;
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);
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();
233 nb = relabelFilter->GetNumberOfObjects();
234 // DD(relabelFilter->GetOriginalNumberOfObjects());
235 // DD(relabelFilter->GetSizeOfObjectsInPhysicalUnits()[0]);
238 typename ImageType::Pointer output = relabelFilter->GetOutput();
241 //--------------------------------------------------------------------
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();
263 return working_image;
265 //--------------------------------------------------------------------
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,
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();
287 //--------------------------------------------------------------------
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)
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,
306 param->GetFirstKeep(),
307 param->GetLastKeep(),
308 param->GetUseLastKeep());
309 return working_image;
311 //--------------------------------------------------------------------
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)
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();
329 //--------------------------------------------------------------------
332 //--------------------------------------------------------------------
333 template<class MaskImageType>
334 typename MaskImageType::Pointer
335 SliceBySliceRelativePosition(const MaskImageType * input,
336 const MaskImageType * object,
339 std::string orientation,
340 bool uniqueConnectedComponent,
343 bool singleObjectCCL)
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->SetUniqueConnectedComponentBySliceFlag(uniqueConnectedComponent);
357 sliceRelPosFilter->ObjectCCLSelectionFlagOff();
358 sliceRelPosFilter->SetUseTheLargestObjectCCLFlag(singleObjectCCL);
359 // sliceRelPosFilter->SetInverseOrientationFlag(inverseflag);
360 sliceRelPosFilter->SetAutoCropFlag(autocropFlag);
361 sliceRelPosFilter->IgnoreEmptySliceObjectFlagOn();
362 sliceRelPosFilter->Update();
363 return sliceRelPosFilter->GetOutput();
365 //--------------------------------------------------------------------
368 //--------------------------------------------------------------------
369 template<class ImageType>
371 FindExtremaPointInAGivenDirection(const ImageType * input,
372 typename ImageType::PixelType bg,
373 int direction, bool opposite,
374 typename ImageType::PointType & point)
376 typename ImageType::PointType dummy;
377 return FindExtremaPointInAGivenDirection(input, bg, direction,
378 opposite, dummy, 0, point);
380 //--------------------------------------------------------------------
383 //--------------------------------------------------------------------
384 template<class ImageType>
386 FindExtremaPointInAGivenDirection(const ImageType * input,
387 typename ImageType::PixelType bg,
388 int direction, bool opposite,
389 typename ImageType::PointType refpoint,
391 typename ImageType::PointType & point)
394 loop over input pixels, store the index in the fg that is max
395 according to the given direction.
397 typedef itk::ImageRegionConstIteratorWithIndex<ImageType> IteratorType;
398 IteratorType iter(input, input->GetLargestPossibleRegion());
400 typename ImageType::IndexType max = input->GetLargestPossibleRegion().GetIndex();
401 if (opposite) max = max+input->GetLargestPossibleRegion().GetSize();
403 while (!iter.IsAtEnd()) {
404 if (iter.Get() != bg) {
405 bool test = iter.GetIndex()[direction] > max[direction];
406 if (opposite) test = !test;
408 typename ImageType::PointType p;
409 input->TransformIndexToPhysicalPoint(iter.GetIndex(), p);
410 if ((distanceMax==0) || (p.EuclideanDistanceTo(refpoint) < distanceMax)) {
411 max = iter.GetIndex();
418 if (!found) return false;
419 input->TransformIndexToPhysicalPoint(max, point);
422 //--------------------------------------------------------------------
425 //--------------------------------------------------------------------
426 template<class ImageType>
427 typename ImageType::Pointer
428 CropImageRemoveGreaterThan(const ImageType * image,
429 int dim, double min, bool autoCrop,
430 typename ImageType::PixelType BG)
432 return CropImageAlongOneAxis<ImageType>(image, dim,
433 image->GetOrigin()[dim],
437 //--------------------------------------------------------------------
440 //--------------------------------------------------------------------
441 template<class ImageType>
442 typename ImageType::Pointer
443 CropImageRemoveLowerThan(const ImageType * image,
444 int dim, double max, bool autoCrop,
445 typename ImageType::PixelType BG)
447 typename ImageType::PointType p;
448 image->TransformIndexToPhysicalPoint(image->GetLargestPossibleRegion().GetIndex()+
449 image->GetLargestPossibleRegion().GetSize(), p);
450 return CropImageAlongOneAxis<ImageType>(image, dim, max, p[dim], autoCrop, BG);
452 //--------------------------------------------------------------------
455 //--------------------------------------------------------------------
456 template<class ImageType>
457 typename ImageType::Pointer
458 CropImageAlongOneAxis(const ImageType * image,
459 int dim, double min, double max,
460 bool autoCrop, typename ImageType::PixelType BG)
462 // Compute region size
463 typename ImageType::RegionType region;
464 typename ImageType::SizeType size = image->GetLargestPossibleRegion().GetSize();
465 typename ImageType::PointType p = image->GetOrigin();
467 typename ImageType::IndexType start;
468 image->TransformPhysicalPointToIndex(p, start);
470 typename ImageType::IndexType end;
471 image->TransformPhysicalPointToIndex(p, end);
472 size[dim] = abs(end[dim]-start[dim]);
473 region.SetIndex(start);
474 region.SetSize(size);
477 typedef itk::RegionOfInterestImageFilter<ImageType, ImageType> CropFilterType;
478 typename CropFilterType::Pointer cropFilter = CropFilterType::New();
479 cropFilter->SetInput(image);
480 cropFilter->SetRegionOfInterest(region);
481 cropFilter->Update();
482 typename ImageType::Pointer result = cropFilter->GetOutput();
486 result = AutoCrop<ImageType>(result, BG);
490 //--------------------------------------------------------------------
493 //--------------------------------------------------------------------
494 template<class ImageType>
496 ComputeCentroids(const ImageType * image,
497 typename ImageType::PixelType BG,
498 std::vector<typename ImageType::PointType> & centroids)
500 typedef long LabelType;
501 static const unsigned int Dim = ImageType::ImageDimension;
502 typedef itk::ShapeLabelObject< LabelType, Dim > LabelObjectType;
503 typedef itk::LabelMap< LabelObjectType > LabelMapType;
504 typedef itk::LabelImageToLabelMapFilter<ImageType, LabelMapType> ImageToMapFilterType;
505 typename ImageToMapFilterType::Pointer imageToLabelFilter = ImageToMapFilterType::New();
506 typedef itk::ShapeLabelMapFilter<LabelMapType, ImageType> ShapeFilterType;
507 typename ShapeFilterType::Pointer statFilter = ShapeFilterType::New();
508 imageToLabelFilter->SetBackgroundValue(BG);
509 imageToLabelFilter->SetInput(image);
510 statFilter->SetInput(imageToLabelFilter->GetOutput());
511 statFilter->Update();
512 typename LabelMapType::Pointer labelMap = statFilter->GetOutput();
515 typename ImageType::PointType dummy;
516 centroids.push_back(dummy); // label 0 -> no centroid, use dummy point for BG
517 //DS FIXME (not useful ! to change ..)
518 for(uint i=0; i<labelMap->GetNumberOfLabelObjects(); i++) {
519 int label = labelMap->GetLabels()[i];
520 centroids.push_back(labelMap->GetLabelObject(label)->GetCentroid());
523 //--------------------------------------------------------------------
526 //--------------------------------------------------------------------
527 template<class ImageType, class LabelType>
528 typename itk::LabelMap< itk::ShapeLabelObject<LabelType, ImageType::ImageDimension> >::Pointer
529 ComputeLabelMap(const ImageType * image,
530 typename ImageType::PixelType BG,
531 bool computePerimeterFlag)
533 static const unsigned int Dim = ImageType::ImageDimension;
534 typedef itk::ShapeLabelObject< LabelType, Dim > LabelObjectType;
535 typedef itk::LabelMap< LabelObjectType > LabelMapType;
536 typedef itk::LabelImageToLabelMapFilter<ImageType, LabelMapType> ImageToMapFilterType;
537 typename ImageToMapFilterType::Pointer imageToLabelFilter = ImageToMapFilterType::New();
538 typedef itk::ShapeLabelMapFilter<LabelMapType, ImageType> ShapeFilterType;
539 typename ShapeFilterType::Pointer statFilter = ShapeFilterType::New();
540 imageToLabelFilter->SetBackgroundValue(BG);
541 imageToLabelFilter->SetInput(image);
542 statFilter->SetInput(imageToLabelFilter->GetOutput());
543 statFilter->SetComputePerimeter(computePerimeterFlag);
544 statFilter->Update();
545 return statFilter->GetOutput();
547 //--------------------------------------------------------------------
550 //--------------------------------------------------------------------
551 template<class ImageType>
553 ComputeCentroids2(const ImageType * image,
554 typename ImageType::PixelType BG,
555 std::vector<typename ImageType::PointType> & centroids)
557 typedef long LabelType;
558 static const unsigned int Dim = ImageType::ImageDimension;
559 typedef itk::ShapeLabelObject< LabelType, Dim > LabelObjectType;
560 typedef itk::LabelMap< LabelObjectType > LabelMapType;
561 typedef itk::LabelImageToLabelMapFilter<ImageType, LabelMapType> ImageToMapFilterType;
562 typename ImageToMapFilterType::Pointer imageToLabelFilter = ImageToMapFilterType::New();
563 typedef itk::ShapeLabelMapFilter<LabelMapType, ImageType> ShapeFilterType;
564 typename ShapeFilterType::Pointer statFilter = ShapeFilterType::New();
565 imageToLabelFilter->SetBackgroundValue(BG);
566 imageToLabelFilter->SetInput(image);
567 statFilter->SetInput(imageToLabelFilter->GetOutput());
568 statFilter->Update();
569 typename LabelMapType::Pointer labelMap = statFilter->GetOutput();
572 typename ImageType::PointType dummy;
573 centroids.push_back(dummy); // label 0 -> no centroid, use dummy point
574 for(uint i=1; i<labelMap->GetNumberOfLabelObjects()+1; i++) {
575 centroids.push_back(labelMap->GetLabelObject(i)->GetCentroid());
578 for(uint i=1; i<labelMap->GetNumberOfLabelObjects()+1; i++) {
579 DD(labelMap->GetLabelObject(i)->GetBinaryPrincipalAxes());
580 DD(labelMap->GetLabelObject(i)->GetBinaryFlatness());
581 DD(labelMap->GetLabelObject(i)->GetRoundness ());
583 // search for the point on the boundary alog PA
588 //--------------------------------------------------------------------
591 //--------------------------------------------------------------------
592 template<class ImageType>
594 PointsUtils<ImageType>::Convert2DTo3D(const PointType2D & p2D,
595 const ImageType * image,
600 index3D[0] = index3D[1] = 0;
601 index3D[2] = image->GetLargestPossibleRegion().GetIndex()[2]+slice;
602 image->TransformIndexToPhysicalPoint(index3D, p3D);
605 // p3D[2] = p[2];//(image->GetLargestPossibleRegion().GetIndex()[2]+slice)*image->GetSpacing()[2]
606 // + image->GetOrigin()[2];
608 //--------------------------------------------------------------------
611 //--------------------------------------------------------------------
612 template<class ImageType>
614 PointsUtils<ImageType>::Convert2DMapTo3DList(const MapPoint2DType & map,
615 const ImageType * image,
616 VectorPoint3DType & list)
618 typename MapPoint2DType::const_iterator iter = map.begin();
619 while (iter != map.end()) {
621 Convert2DTo3D(iter->second, image, iter->first, p);
626 //--------------------------------------------------------------------
629 //--------------------------------------------------------------------
630 template<class ImageType>
632 PointsUtils<ImageType>::Convert2DListTo3DList(const VectorPoint2DType & p2D,
634 const ImageType * image,
635 VectorPoint3DType & list)
637 for(uint i=0; i<p2D.size(); i++) {
639 Convert2DTo3D(p2D[i], image, slice, p);
643 //--------------------------------------------------------------------
646 //--------------------------------------------------------------------
647 template<class ImageType>
649 WriteListOfLandmarks(std::vector<typename ImageType::PointType> points,
650 std::string filename)
653 openFileForWriting(os, filename);
654 os << "LANDMARKS1" << std::endl;
655 for(uint i=0; i<points.size(); i++) {
656 const typename ImageType::PointType & p = points[i];
657 // Write it in the file
658 os << i << " " << p[0] << " " << p[1] << " " << p[2] << " 0 0 " << std::endl;
662 //--------------------------------------------------------------------
665 //--------------------------------------------------------------------
666 template<class ImageType>
667 typename ImageType::Pointer
668 Dilate(const ImageType * image, double radiusInMM,
669 typename ImageType::PixelType BG,
670 typename ImageType::PixelType FG,
673 typename ImageType::SizeType r;
674 for(uint i=0; i<ImageType::ImageDimension; i++)
675 r[i] = (uint)lrint(radiusInMM/image->GetSpacing()[i]);
676 return Dilate<ImageType>(image, r, BG, FG, extendSupport);
678 //--------------------------------------------------------------------
681 //--------------------------------------------------------------------
682 template<class ImageType>
683 typename ImageType::Pointer
684 Dilate(const ImageType * image, typename ImageType::PointType radiusInMM,
685 typename ImageType::PixelType BG,
686 typename ImageType::PixelType FG,
689 typename ImageType::SizeType r;
690 for(uint i=0; i<ImageType::ImageDimension; i++)
691 r[i] = (uint)lrint(radiusInMM[i]/image->GetSpacing()[i]);
692 return Dilate<ImageType>(image, r, BG, FG, extendSupport);
694 //--------------------------------------------------------------------
697 //--------------------------------------------------------------------
698 template<class ImageType>
699 typename ImageType::Pointer
700 Dilate(const ImageType * image, typename ImageType::SizeType radius,
701 typename ImageType::PixelType BG,
702 typename ImageType::PixelType FG,
705 // Create kernel for dilatation
706 typedef itk::BinaryBallStructuringElement<typename ImageType::PixelType,
707 ImageType::ImageDimension> KernelType;
708 KernelType structuringElement;
709 structuringElement.SetRadius(radius);
710 structuringElement.CreateStructuringElement();
712 typename ImageType::Pointer output;
714 typedef itk::ConstantPadImageFilter<ImageType, ImageType> PadFilterType;
715 typename PadFilterType::Pointer padFilter = PadFilterType::New();
716 padFilter->SetInput(image);
717 typename ImageType::SizeType lower;
718 typename ImageType::SizeType upper;
719 for(uint i=0; i<3; i++) {
720 lower[i] = upper[i] = 2*(radius[i]+1);
722 padFilter->SetPadLowerBound(lower);
723 padFilter->SetPadUpperBound(upper);
725 output = padFilter->GetOutput();
729 typedef itk::BinaryDilateImageFilter<ImageType, ImageType , KernelType> DilateFilterType;
730 typename DilateFilterType::Pointer dilateFilter = DilateFilterType::New();
731 dilateFilter->SetBackgroundValue(BG);
732 dilateFilter->SetForegroundValue(FG);
733 dilateFilter->SetBoundaryToForeground(false);
734 dilateFilter->SetKernel(structuringElement);
735 if (extendSupport) dilateFilter->SetInput(output);
736 else dilateFilter->SetInput(image);
737 dilateFilter->Update();
738 return dilateFilter->GetOutput();
740 //--------------------------------------------------------------------
743 //--------------------------------------------------------------------
744 template<class ImageType>
745 typename ImageType::Pointer
746 Opening(const ImageType * image, typename ImageType::SizeType radius,
747 typename ImageType::PixelType BG,
748 typename ImageType::PixelType FG)
751 typedef itk::BinaryBallStructuringElement<typename ImageType::PixelType,
752 ImageType::ImageDimension> KernelType;
753 KernelType structuringElement;
754 structuringElement.SetRadius(radius);
755 structuringElement.CreateStructuringElement();
758 typedef itk::BinaryMorphologicalOpeningImageFilter<ImageType, ImageType , KernelType> OpeningFilterType;
759 typename OpeningFilterType::Pointer open = OpeningFilterType::New();
760 open->SetInput(image);
761 open->SetBackgroundValue(BG);
762 open->SetForegroundValue(FG);
763 open->SetKernel(structuringElement);
765 return open->GetOutput();
767 //--------------------------------------------------------------------
771 //--------------------------------------------------------------------
772 template<class ValueType, class VectorType>
773 void ConvertOption(std::string optionName, uint given,
774 ValueType * values, VectorType & p,
775 uint dim, bool required)
777 if (required && (given == 0)) {
778 clitkExceptionMacro("The option --" << optionName << " must be set and have 1 or "
779 << dim << " values.");
782 for(uint i=0; i<dim; i++) p[i] = values[0];
786 for(uint i=0; i<dim; i++) p[i] = values[i];
789 if (given == 0) return;
790 clitkExceptionMacro("The option --" << optionName << " must have 1 or "
791 << dim << " values.");
793 //--------------------------------------------------------------------
796 //--------------------------------------------------------------------
798 http://www.gamedev.net/community/forums/topic.asp?topic_id=542870
799 Assuming the points are (Ax,Ay) (Bx,By) and (Cx,Cy), you need to compute:
800 (Bx - Ax) * (Cy - Ay) - (By - Ay) * (Cx - Ax)
801 This will equal zero if the point C is on the line formed by
802 points A and B, and will have a different sign depending on the
803 side. Which side this is depends on the orientation of your (x,y)
804 coordinates, but you can plug test values for A,B and C into this
805 formula to determine whether negative values are to the left or to
807 => to accelerate, start with formula, when change sign -> stop and fill
809 offsetToKeep = is used to determine which side of the line we
810 keep. The point along the mainDirection but 'offsetToKeep' mm away
814 template<class ImageType>
816 SliceBySliceSetBackgroundFromLineSeparation(ImageType * input,
817 std::vector<typename ImageType::PointType> & lA,
818 std::vector<typename ImageType::PointType> & lB,
819 typename ImageType::PixelType BG,
823 assert((mainDirection==0) || (mainDirection==1));
824 typedef itk::ImageSliceIteratorWithIndex<ImageType> SliceIteratorType;
825 SliceIteratorType siter = SliceIteratorType(input,
826 input->GetLargestPossibleRegion());
827 siter.SetFirstDirection(0);
828 siter.SetSecondDirection(1);
831 typename ImageType::PointType A;
832 typename ImageType::PointType B;
833 typename ImageType::PointType C;
834 assert(lA.size() == lB.size());
835 while ((i<lA.size()) && (!siter.IsAtEnd())) {
836 // Check that the current slice correspond to the current point
837 input->TransformIndexToPhysicalPoint(siter.GetIndex(), C);
838 if ((fabs(C[2] - lA[i][2]))>0.01) { // is !equal with a tolerance of 0.01 mm
841 // Define A,B,C points
846 // Check that the line is not a point (A=B)
847 bool p = (A[0] == B[0]) && (A[1] == B[1]);
850 C[mainDirection] += offsetToKeep; // I know I must keep this point
851 double s = (B[0] - A[0]) * (C[1] - A[1]) - (B[1] - A[1]) * (C[0] - A[0]);
852 bool isPositive = s<0;
853 while (!siter.IsAtEndOfSlice()) {
854 while (!siter.IsAtEndOfLine()) {
855 // Very slow, I know ... but image should be very small
856 input->TransformIndexToPhysicalPoint(siter.GetIndex(), C);
857 double s = (B[0] - A[0]) * (C[1] - A[1]) - (B[1] - A[1]) * (C[0] - A[0]);
858 if (s == 0) siter.Set(BG); // on the line, we decide to remove
860 if (s > 0) siter.Set(BG);
863 if (s < 0) siter.Set(BG);
872 } // End of current slice
876 //--------------------------------------------------------------------
879 //--------------------------------------------------------------------
880 template<class ImageType>
882 AndNot(ImageType * input,
883 const ImageType * object,
884 typename ImageType::PixelType BG)
886 typename ImageType::Pointer o;
888 if (!clitk::HaveSameSizeAndSpacing<ImageType, ImageType>(input, object)) {
889 o = clitk::ResizeImageLike<ImageType>(object, input, BG);
893 typedef clitk::BooleanOperatorLabelImageFilter<ImageType> BoolFilterType;
894 typename BoolFilterType::Pointer boolFilter = BoolFilterType::New();
895 boolFilter->InPlaceOn();
896 boolFilter->SetInput1(input);
897 if (resized) boolFilter->SetInput2(o);
898 else boolFilter->SetInput2(object);
899 boolFilter->SetBackgroundValue1(BG);
900 boolFilter->SetBackgroundValue2(BG);
901 boolFilter->SetOperationType(BoolFilterType::AndNot);
902 boolFilter->Update();
904 //--------------------------------------------------------------------
907 //--------------------------------------------------------------------
908 template<class ImageType>
910 And(ImageType * input,
911 const ImageType * object,
912 typename ImageType::PixelType BG)
914 typename ImageType::Pointer o;
916 if (!clitk::HaveSameSizeAndSpacing<ImageType, ImageType>(input, object)) {
917 o = clitk::ResizeImageLike<ImageType>(object, input, BG);
921 typedef clitk::BooleanOperatorLabelImageFilter<ImageType> BoolFilterType;
922 typename BoolFilterType::Pointer boolFilter = BoolFilterType::New();
923 boolFilter->InPlaceOn();
924 boolFilter->SetInput1(input);
925 if (resized) boolFilter->SetInput2(o);
926 else boolFilter->SetInput2(object);
927 boolFilter->SetBackgroundValue1(BG);
928 boolFilter->SetBackgroundValue2(BG);
929 boolFilter->SetOperationType(BoolFilterType::And);
930 boolFilter->Update();
932 //--------------------------------------------------------------------
935 //--------------------------------------------------------------------
936 template<class ImageType>
938 Or(ImageType * input,
939 const ImageType * object,
940 typename ImageType::PixelType BG)
942 typename ImageType::Pointer o;
944 if (!clitk::HaveSameSizeAndSpacing<ImageType, ImageType>(input, object)) {
945 o = clitk::ResizeImageLike<ImageType>(object, input, BG);
949 typedef clitk::BooleanOperatorLabelImageFilter<ImageType> BoolFilterType;
950 typename BoolFilterType::Pointer boolFilter = BoolFilterType::New();
951 boolFilter->InPlaceOn();
952 boolFilter->SetInput1(input);
953 if (resized) boolFilter->SetInput2(o);
954 else boolFilter->SetInput2(object);
955 boolFilter->SetBackgroundValue1(BG);
956 boolFilter->SetBackgroundValue2(BG);
957 boolFilter->SetOperationType(BoolFilterType::Or);
958 boolFilter->Update();
960 //--------------------------------------------------------------------
963 //--------------------------------------------------------------------
964 template<class ImageType>
965 typename ImageType::Pointer
966 Binarize(const ImageType * input,
967 typename ImageType::PixelType lower,
968 typename ImageType::PixelType upper,
969 typename ImageType::PixelType BG,
970 typename ImageType::PixelType FG)
972 typedef itk::BinaryThresholdImageFilter<ImageType, ImageType> BinaryThresholdFilterType;
973 typename BinaryThresholdFilterType::Pointer binarizeFilter = BinaryThresholdFilterType::New();
974 binarizeFilter->SetInput(input);
975 binarizeFilter->InPlaceOff();
976 binarizeFilter->SetLowerThreshold(lower);
977 binarizeFilter->SetUpperThreshold(upper);
978 binarizeFilter->SetInsideValue(FG);
979 binarizeFilter->SetOutsideValue(BG);
980 binarizeFilter->Update();
981 return binarizeFilter->GetOutput();
983 //--------------------------------------------------------------------
986 //--------------------------------------------------------------------
987 template<class ImageType>
989 GetMinMaxPointPosition(const ImageType * input,
990 typename ImageType::PointType & min,
991 typename ImageType::PointType & max)
993 typename ImageType::IndexType index = input->GetLargestPossibleRegion().GetIndex();
994 input->TransformIndexToPhysicalPoint(index, min);
995 index = index+input->GetLargestPossibleRegion().GetSize();
996 input->TransformIndexToPhysicalPoint(index, max);
998 //--------------------------------------------------------------------
1001 //--------------------------------------------------------------------
1002 template<class ImageType>
1003 typename ImageType::PointType
1004 FindExtremaPointInAGivenLine(const ImageType * input,
1007 typename ImageType::PointType p,
1008 typename ImageType::PixelType BG,
1011 // Which direction ? Increasing or decreasing.
1015 // Transform to pixel index
1016 typename ImageType::IndexType index;
1017 input->TransformPhysicalPointToIndex(p, index);
1019 // Loop while inside the mask;
1020 while (input->GetPixel(index) != BG) {
1021 index[dimension] += d;
1024 // Transform back to Physical Units
1025 typename ImageType::PointType result;
1026 input->TransformIndexToPhysicalPoint(index, result);
1028 // Check that is is not too far away
1029 double distance = p.EuclideanDistanceTo(result);
1030 if (distance > distanceMax) {
1031 result = p; // Get back to initial value
1036 //--------------------------------------------------------------------
1039 //--------------------------------------------------------------------
1040 template<class PointType>
1042 IsOnTheSameLineSide(PointType C, PointType A, PointType B, PointType like)
1044 // Look at the position of point 'like' according to the AB line
1045 double s = (B[0] - A[0]) * (like[1] - A[1]) - (B[1] - A[1]) * (like[0] - A[0]);
1046 bool negative = s<0;
1048 // Look the C position
1049 s = (B[0] - A[0]) * (C[1] - A[1]) - (B[1] - A[1]) * (C[0] - A[0]);
1051 if (negative && (s<=0)) return true;
1052 if (!negative && (s>=0)) return true;
1055 //--------------------------------------------------------------------
1058 //--------------------------------------------------------------------
1059 /* Consider an input object, for each slice, find the extrema
1060 position according to a given direction and build a line segment
1061 passing throught this point in a given direction. Output is a
1062 vector of line (from point A to B), for each slice;
1064 template<class ImageType>
1066 SliceBySliceBuildLineSegmentAccordingToExtremaPosition(const ImageType * input,
1067 typename ImageType::PixelType BG,
1069 int extremaDirection,
1070 bool extremaOppositeFlag,
1073 std::vector<typename ImageType::PointType> & A,
1074 std::vector<typename ImageType::PointType> & B)
1077 typedef typename itk::Image<typename ImageType::PixelType, ImageType::ImageDimension-1> SliceType;
1079 // Build the list of slices
1080 std::vector<typename SliceType::Pointer> slices;
1081 clitk::ExtractSlices<ImageType>(input, sliceDimension, slices);
1083 // Build the list of 2D points
1084 std::map<int, typename SliceType::PointType> position2D;
1085 for(uint i=0; i<slices.size(); i++) {
1086 typename SliceType::PointType p;
1088 clitk::FindExtremaPointInAGivenDirection<SliceType>(slices[i], BG,
1089 extremaDirection, extremaOppositeFlag, p);
1095 // Convert 2D points in slice into 3D points
1096 clitk::PointsUtils<ImageType>::Convert2DMapTo3DList(position2D, input, A);
1098 // Create additional point just right to the previous ones, on the
1099 // given lineDirection, in order to create a horizontal/vertical line.
1100 for(uint i=0; i<A.size(); i++) {
1101 typename ImageType::PointType p = A[i];
1102 p[lineDirection] += 10;
1105 A[i][extremaDirection] += margin;
1106 B[i][extremaDirection] += margin;
1110 //--------------------------------------------------------------------
1113 //--------------------------------------------------------------------
1114 template<class ImageType>
1115 typename ImageType::Pointer
1116 SliceBySliceKeepMainCCL(const ImageType * input,
1117 typename ImageType::PixelType BG,
1118 typename ImageType::PixelType FG) {
1121 const int d = ImageType::ImageDimension-1;
1122 typedef typename itk::Image<typename ImageType::PixelType, d> SliceType;
1123 std::vector<typename SliceType::Pointer> slices;
1124 clitk::ExtractSlices<ImageType>(input, d, slices);
1126 // Labelize and keep the main one
1127 std::vector<typename SliceType::Pointer> o;
1128 for(uint i=0; i<slices.size(); i++) {
1129 o.push_back(clitk::Labelize<SliceType>(slices[i], BG, false, 1));
1130 o[i] = clitk::KeepLabels<SliceType>(o[i], BG, FG, 1, 1, true);
1134 typename ImageType::Pointer output;
1135 output = clitk::JoinSlices<ImageType>(o, input, d);
1138 //--------------------------------------------------------------------
1141 //--------------------------------------------------------------------
1142 template<class ImageType>
1143 typename ImageType::Pointer
1144 Clone(const ImageType * input) {
1145 typedef itk::ImageDuplicator<ImageType> DuplicatorType;
1146 typename DuplicatorType::Pointer duplicator = DuplicatorType::New();
1147 duplicator->SetInputImage(input);
1148 duplicator->Update();
1149 return duplicator->GetOutput();
1151 //--------------------------------------------------------------------
1154 //--------------------------------------------------------------------
1155 /* Consider an input object, start at A, for each slice (dim1):
1156 - compute the intersection between the AB line and the current slice
1157 - remove what is at lower or greater according to dim2 of this point
1160 template<class ImageType>
1161 typename ImageType::Pointer
1162 SliceBySliceSetBackgroundFromSingleLine(const ImageType * input,
1163 typename ImageType::PixelType BG,
1164 typename ImageType::PointType & A,
1165 typename ImageType::PointType & B,
1166 int dim1, int dim2, bool removeLowerPartFlag)
1170 typedef typename itk::Image<typename ImageType::PixelType, ImageType::ImageDimension-1> SliceType;
1171 typedef typename SliceType::Pointer SlicePointer;
1172 std::vector<SlicePointer> slices;
1173 clitk::ExtractSlices<ImageType>(input, dim1, slices);
1175 // Start at slice that contains A, and stop at B
1176 typename ImageType::IndexType Ap;
1177 typename ImageType::IndexType Bp;
1178 input->TransformPhysicalPointToIndex(A, Ap);
1179 input->TransformPhysicalPointToIndex(B, Bp);
1181 // Determine slice largest region
1182 typename SliceType::RegionType region = slices[0]->GetLargestPossibleRegion();
1183 typename SliceType::SizeType size = region.GetSize();
1184 typename SliceType::IndexType index = region.GetIndex();
1187 double a = (Bp[dim2]-Ap[dim2])/(Bp[dim1]-Ap[dim1]);
1188 double b = Ap[dim2];
1190 // Loop from slice A to slice B
1191 for(uint i=0; i<(Bp[dim1]-Ap[dim1]); i++) {
1192 // Compute intersection between line AB and current slice for the dim2
1194 // Change region (lower than dim2)
1195 if (removeLowerPartFlag) {
1196 size[dim2] = p-Ap[dim2];
1199 size[dim2] = slices[0]->GetLargestPossibleRegion().GetSize()[dim2]-p;
1202 region.SetSize(size);
1203 region.SetIndex(index);
1204 // Fill region with BG (simple region iterator)
1205 FillRegionWithValue<SliceType>(slices[i+Ap[dim1]], BG, region);
1207 typedef itk::ImageRegionIterator<SliceType> IteratorType;
1208 IteratorType iter(slices[i+Ap[dim1]], region);
1210 while (!iter.IsAtEnd()) {
1219 typename ImageType::Pointer output;
1220 output = clitk::JoinSlices<ImageType>(slices, input, dim1);
1223 //--------------------------------------------------------------------
1225 //--------------------------------------------------------------------
1226 /* Consider an input object, slice by slice, use the point A and set
1227 pixel to BG according to their position relatively to A
1229 template<class ImageType>
1230 typename ImageType::Pointer
1231 SliceBySliceSetBackgroundFromPoints(const ImageType * input,
1232 typename ImageType::PixelType BG,
1234 std::vector<typename ImageType::PointType> & A,
1235 bool removeGreaterThanXFlag,
1236 bool removeGreaterThanYFlag)
1240 typedef typename itk::Image<typename ImageType::PixelType, ImageType::ImageDimension-1> SliceType;
1241 typedef typename SliceType::Pointer SlicePointer;
1242 std::vector<SlicePointer> slices;
1243 clitk::ExtractSlices<ImageType>(input, sliceDim, slices);
1245 // Start at slice that contains A
1246 typename ImageType::IndexType Ap;
1248 // Determine slice largest region
1249 typename SliceType::RegionType region = slices[0]->GetLargestPossibleRegion();
1250 typename SliceType::SizeType size = region.GetSize();
1251 typename SliceType::IndexType index = region.GetIndex();
1253 // Loop from slice A to slice B
1254 for(uint i=0; i<A.size(); i++) {
1255 input->TransformPhysicalPointToIndex(A[i], Ap);
1256 uint sliceIndex = Ap[2] - input->GetLargestPossibleRegion().GetIndex()[2];
1257 if ((sliceIndex < 0) || (sliceIndex >= slices.size())) {
1258 continue; // do not consider this slice
1261 // Compute region for BG
1262 if (removeGreaterThanXFlag) {
1264 size[0] = region.GetSize()[0]-(index[0]-region.GetIndex()[0]);
1267 index[0] = region.GetIndex()[0];
1268 size[0] = Ap[0] - index[0];
1271 if (removeGreaterThanYFlag) {
1273 size[1] = region.GetSize()[1]-(index[1]-region.GetIndex()[1]);
1276 index[1] = region.GetIndex()[1];
1277 size[1] = Ap[1] - index[1];
1281 region.SetSize(size);
1282 region.SetIndex(index);
1284 // Fill region with BG (simple region iterator)
1285 FillRegionWithValue<SliceType>(slices[sliceIndex], BG, region);
1290 typename ImageType::Pointer output;
1291 output = clitk::JoinSlices<ImageType>(slices, input, sliceDim);
1294 //--------------------------------------------------------------------
1297 //--------------------------------------------------------------------
1298 template<class ImageType>
1300 FillRegionWithValue(ImageType * input, typename ImageType::PixelType value, typename ImageType::RegionType & region)
1302 typedef itk::ImageRegionIterator<ImageType> IteratorType;
1303 IteratorType iter(input, region);
1305 while (!iter.IsAtEnd()) {
1310 //--------------------------------------------------------------------
1313 //--------------------------------------------------------------------
1314 template<class ImageType>
1316 GetMinMaxBoundary(ImageType * input, typename ImageType::PointType & min,
1317 typename ImageType::PointType & max)
1319 typedef typename ImageType::PointType PointType;
1320 typedef typename ImageType::IndexType IndexType;
1321 IndexType min_i, max_i;
1322 min_i = input->GetLargestPossibleRegion().GetIndex();
1323 for(uint i=0; i<ImageType::ImageDimension; i++)
1324 max_i[i] = input->GetLargestPossibleRegion().GetSize()[i] + min_i[i];
1325 input->TransformIndexToPhysicalPoint(min_i, min);
1326 input->TransformIndexToPhysicalPoint(max_i, max);
1328 //--------------------------------------------------------------------
1330 } // end of namespace