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>
39 //--------------------------------------------------------------------
40 template<class ImageType>
41 void ComputeBBFromImageRegion(const ImageType * image,
42 typename ImageType::RegionType region,
43 typename itk::BoundingBox<unsigned long,
44 ImageType::ImageDimension>::Pointer bb) {
45 typedef typename ImageType::IndexType IndexType;
48 for(unsigned int i=0; i<image->GetImageDimension(); i++) {
49 firstIndex[i] = region.GetIndex()[i];
50 lastIndex[i] = firstIndex[i]+region.GetSize()[i];
53 typedef itk::BoundingBox<unsigned long,
54 ImageType::ImageDimension> BBType;
55 typedef typename BBType::PointType PointType;
58 image->TransformIndexToPhysicalPoint(firstIndex, firstPoint);
59 image->TransformIndexToPhysicalPoint(lastIndex, lastPoint);
61 bb->SetMaximum(lastPoint);
62 bb->SetMinimum(firstPoint);
64 //--------------------------------------------------------------------
67 //--------------------------------------------------------------------
68 template<int Dimension>
69 void ComputeBBIntersection(typename itk::BoundingBox<unsigned long, Dimension>::Pointer bbo,
70 typename itk::BoundingBox<unsigned long, Dimension>::Pointer bbi1,
71 typename itk::BoundingBox<unsigned long, Dimension>::Pointer bbi2) {
73 typedef itk::BoundingBox<unsigned long, Dimension> BBType;
74 typedef typename BBType::PointType PointType;
78 for(unsigned int i=0; i<Dimension; i++) {
79 firstPoint[i] = std::max(bbi1->GetMinimum()[i],
80 bbi2->GetMinimum()[i]);
81 lastPoint[i] = std::min(bbi1->GetMaximum()[i],
82 bbi2->GetMaximum()[i]);
85 bbo->SetMaximum(lastPoint);
86 bbo->SetMinimum(firstPoint);
88 //--------------------------------------------------------------------
91 //--------------------------------------------------------------------
92 template<class ImageType>
93 void ComputeRegionFromBB(const ImageType * image,
94 const typename itk::BoundingBox<unsigned long,
95 ImageType::ImageDimension>::Pointer bb,
96 typename ImageType::RegionType & region) {
98 typedef typename ImageType::IndexType IndexType;
99 typedef typename ImageType::PointType PointType;
100 typedef typename ImageType::RegionType RegionType;
101 typedef typename ImageType::SizeType SizeType;
103 // Region starting point
104 IndexType regionStart;
105 PointType start = bb->GetMinimum();
106 image->TransformPhysicalPointToIndex(start, regionStart);
110 PointType maxs = bb->GetMaximum();
111 PointType mins = bb->GetMinimum();
112 for(unsigned int i=0; i<ImageType::ImageDimension; i++) {
113 regionSize[i] = lrint((maxs[i] - mins[i])/image->GetSpacing()[i]);
117 region.SetIndex(regionStart);
118 region.SetSize(regionSize);
120 //--------------------------------------------------------------------
122 //--------------------------------------------------------------------
123 template<class ImageType, class TMaskImageType>
124 typename ImageType::Pointer
125 SetBackground(const ImageType * input,
126 const TMaskImageType * mask,
127 typename TMaskImageType::PixelType maskBG,
128 typename ImageType::PixelType outValue,
130 typedef SetBackgroundImageFilter<ImageType, TMaskImageType, ImageType>
131 SetBackgroundImageFilterType;
132 typename SetBackgroundImageFilterType::Pointer setBackgroundFilter
133 = SetBackgroundImageFilterType::New();
134 // if (inPlace) setBackgroundFilter->ReleaseDataFlagOn(); // No seg fault
135 setBackgroundFilter->SetInPlace(inPlace); // This is important to keep memory low
136 setBackgroundFilter->SetInput(input);
137 setBackgroundFilter->SetInput2(mask);
138 setBackgroundFilter->SetMaskValue(maskBG);
139 setBackgroundFilter->SetOutsideValue(outValue);
140 setBackgroundFilter->Update();
141 return setBackgroundFilter->GetOutput();
143 //--------------------------------------------------------------------
146 //--------------------------------------------------------------------
147 template<class ImageType>
148 int GetNumberOfConnectedComponentLabels(const ImageType * input,
149 typename ImageType::PixelType BG,
150 bool isFullyConnected) {
151 // Connected Component label
152 typedef itk::ConnectedComponentImageFilter<ImageType, ImageType> ConnectFilterType;
153 typename ConnectFilterType::Pointer connectFilter = ConnectFilterType::New();
154 connectFilter->SetInput(input);
155 connectFilter->SetBackgroundValue(BG);
156 connectFilter->SetFullyConnected(isFullyConnected);
157 connectFilter->Update();
160 return connectFilter->GetObjectCount();
162 //--------------------------------------------------------------------
164 //--------------------------------------------------------------------
166 Warning : in this cas, we consider outputType like inputType, not
167 InternalImageType. Be sure it fits.
169 template<class ImageType>
170 typename ImageType::Pointer
171 Labelize(const ImageType * input,
172 typename ImageType::PixelType BG,
173 bool isFullyConnected,
174 int minimalComponentSize) {
175 // InternalImageType for storing large number of component
176 typedef itk::Image<int, ImageType::ImageDimension> InternalImageType;
178 // Connected Component label
179 typedef itk::ConnectedComponentImageFilter<ImageType, InternalImageType> ConnectFilterType;
180 typename ConnectFilterType::Pointer connectFilter = ConnectFilterType::New();
181 // connectFilter->ReleaseDataFlagOn();
182 connectFilter->SetInput(input);
183 connectFilter->SetBackgroundValue(BG);
184 connectFilter->SetFullyConnected(isFullyConnected);
186 // Sort by size and remove too small area.
187 typedef itk::RelabelComponentImageFilter<InternalImageType, ImageType> RelabelFilterType;
188 typename RelabelFilterType::Pointer relabelFilter = RelabelFilterType::New();
189 // relabelFilter->ReleaseDataFlagOn(); // if yes, fail when ExplosionControlledThresholdConnectedImageFilter ???
190 relabelFilter->SetInput(connectFilter->GetOutput());
191 relabelFilter->SetMinimumObjectSize(minimalComponentSize);
192 relabelFilter->Update();
195 typename ImageType::Pointer output = relabelFilter->GetOutput();
198 //--------------------------------------------------------------------
201 //--------------------------------------------------------------------
203 Warning : in this cas, we consider outputType like inputType, not
204 InternalImageType. Be sure it fits.
206 template<class ImageType>
207 typename ImageType::Pointer
208 LabelizeAndCountNumberOfObjects(const ImageType * input,
209 typename ImageType::PixelType BG,
210 bool isFullyConnected,
211 int minimalComponentSize,
213 // InternalImageType for storing large number of component
214 typedef itk::Image<int, ImageType::ImageDimension> InternalImageType;
216 // Connected Component label
217 typedef itk::ConnectedComponentImageFilter<ImageType, InternalImageType> ConnectFilterType;
218 typename ConnectFilterType::Pointer connectFilter = ConnectFilterType::New();
219 // connectFilter->ReleaseDataFlagOn();
220 connectFilter->SetInput(input);
221 connectFilter->SetBackgroundValue(BG);
222 connectFilter->SetFullyConnected(isFullyConnected);
224 // Sort by size and remove too small area.
225 typedef itk::RelabelComponentImageFilter<InternalImageType, ImageType> RelabelFilterType;
226 typename RelabelFilterType::Pointer relabelFilter = RelabelFilterType::New();
227 // relabelFilter->ReleaseDataFlagOn(); // if yes, fail when ExplosionControlledThresholdConnectedImageFilter ???
228 relabelFilter->SetInput(connectFilter->GetOutput());
229 relabelFilter->SetMinimumObjectSize(minimalComponentSize);
230 relabelFilter->Update();
232 nb = relabelFilter->GetNumberOfObjects();
233 // DD(relabelFilter->GetOriginalNumberOfObjects());
234 // DD(relabelFilter->GetSizeOfObjectsInPhysicalUnits()[0]);
237 typename ImageType::Pointer output = relabelFilter->GetOutput();
240 //--------------------------------------------------------------------
244 //--------------------------------------------------------------------
245 template<class ImageType>
246 typename ImageType::Pointer
247 RemoveLabels(const ImageType * input,
248 typename ImageType::PixelType BG,
249 std::vector<typename ImageType::PixelType> & labelsToRemove) {
250 assert(labelsToRemove.size() != 0);
251 typename ImageType::Pointer working_image;// = input;
252 for (unsigned int i=0; i <labelsToRemove.size(); i++) {
253 typedef SetBackgroundImageFilter<ImageType, ImageType> SetBackgroundImageFilterType;
254 typename SetBackgroundImageFilterType::Pointer setBackgroundFilter = SetBackgroundImageFilterType::New();
255 setBackgroundFilter->SetInput(input);
256 setBackgroundFilter->SetInput2(input);
257 setBackgroundFilter->SetMaskValue(labelsToRemove[i]);
258 setBackgroundFilter->SetOutsideValue(BG);
259 setBackgroundFilter->Update();
260 working_image = setBackgroundFilter->GetOutput();
262 return working_image;
264 //--------------------------------------------------------------------
267 //--------------------------------------------------------------------
268 template<class ImageType>
269 typename ImageType::Pointer
270 KeepLabels(const ImageType * input,
271 typename ImageType::PixelType BG,
272 typename ImageType::PixelType FG,
273 typename ImageType::PixelType firstKeep,
274 typename ImageType::PixelType lastKeep,
276 typedef itk::BinaryThresholdImageFilter<ImageType, ImageType> BinarizeFilterType;
277 typename BinarizeFilterType::Pointer binarizeFilter = BinarizeFilterType::New();
278 binarizeFilter->SetInput(input);
279 binarizeFilter->SetLowerThreshold(firstKeep);
280 if (useLastKeep) binarizeFilter->SetUpperThreshold(lastKeep);
281 binarizeFilter->SetInsideValue(FG);
282 binarizeFilter->SetOutsideValue(BG);
283 binarizeFilter->Update();
284 return binarizeFilter->GetOutput();
286 //--------------------------------------------------------------------
289 //--------------------------------------------------------------------
290 template<class ImageType>
291 typename ImageType::Pointer
292 LabelizeAndSelectLabels(const ImageType * input,
293 typename ImageType::PixelType BG,
294 typename ImageType::PixelType FG,
295 bool isFullyConnected,
296 int minimalComponentSize,
297 LabelizeParameters<typename ImageType::PixelType> * param)
299 typename ImageType::Pointer working_image;
300 working_image = Labelize<ImageType>(input, BG, isFullyConnected, minimalComponentSize);
301 if (param->GetLabelsToRemove().size() != 0)
302 working_image = RemoveLabels<ImageType>(working_image, BG, param->GetLabelsToRemove());
303 working_image = KeepLabels<ImageType>(working_image,
305 param->GetFirstKeep(),
306 param->GetLastKeep(),
307 param->GetUseLastKeep());
308 return working_image;
310 //--------------------------------------------------------------------
313 //--------------------------------------------------------------------
314 template<class ImageType>
315 typename ImageType::Pointer
316 ResizeImageLike(const ImageType * input,
317 const itk::ImageBase<ImageType::ImageDimension> * like,
318 typename ImageType::PixelType backgroundValue)
320 typedef CropLikeImageFilter<ImageType> CropFilterType;
321 typename CropFilterType::Pointer cropFilter = CropFilterType::New();
322 cropFilter->SetInput(input);
323 cropFilter->SetCropLikeImage(like);
324 cropFilter->SetBackgroundValue(backgroundValue);
325 cropFilter->Update();
326 return cropFilter->GetOutput();
328 //--------------------------------------------------------------------
331 //--------------------------------------------------------------------
332 template<class MaskImageType>
333 typename MaskImageType::Pointer
334 SliceBySliceRelativePosition(const MaskImageType * input,
335 const MaskImageType * object,
338 std::string orientation,
339 bool uniqueConnectedComponent,
342 bool singleObjectCCL)
344 typedef SliceBySliceRelativePositionFilter<MaskImageType> SliceRelPosFilterType;
345 typename SliceRelPosFilterType::Pointer sliceRelPosFilter = SliceRelPosFilterType::New();
346 sliceRelPosFilter->VerboseStepFlagOff();
347 sliceRelPosFilter->WriteStepFlagOff();
348 sliceRelPosFilter->SetInput(input);
349 sliceRelPosFilter->SetInputObject(object);
350 sliceRelPosFilter->SetDirection(direction);
351 sliceRelPosFilter->SetFuzzyThreshold(threshold);
352 sliceRelPosFilter->AddOrientationTypeString(orientation);
353 sliceRelPosFilter->SetIntermediateSpacingFlag((spacing != -1));
354 sliceRelPosFilter->SetIntermediateSpacing(spacing);
355 sliceRelPosFilter->SetUniqueConnectedComponentBySlice(uniqueConnectedComponent);
356 sliceRelPosFilter->SetUseASingleObjectConnectedComponentBySliceFlag(singleObjectCCL);
357 // sliceRelPosFilter->SetInverseOrientationFlag(inverseflag);
358 sliceRelPosFilter->SetAutoCropFlag(autocropFlag);
359 sliceRelPosFilter->IgnoreEmptySliceObjectFlagOn();
360 sliceRelPosFilter->Update();
361 return sliceRelPosFilter->GetOutput();
363 //--------------------------------------------------------------------
365 //--------------------------------------------------------------------
366 template<class ImageType>
368 FindExtremaPointInAGivenDirection(const ImageType * input,
369 typename ImageType::PixelType bg,
370 int direction, bool opposite,
371 typename ImageType::PointType & point)
373 typename ImageType::PointType dummy;
374 return FindExtremaPointInAGivenDirection(input, bg, direction,
375 opposite, dummy, 0, point);
377 //--------------------------------------------------------------------
379 //--------------------------------------------------------------------
380 template<class ImageType>
382 FindExtremaPointInAGivenDirection(const ImageType * input,
383 typename ImageType::PixelType bg,
384 int direction, bool opposite,
385 typename ImageType::PointType refpoint,
387 typename ImageType::PointType & point)
390 loop over input pixels, store the index in the fg that is max
391 according to the given direction.
393 typedef itk::ImageRegionConstIteratorWithIndex<ImageType> IteratorType;
394 IteratorType iter(input, input->GetLargestPossibleRegion());
396 typename ImageType::IndexType max = input->GetLargestPossibleRegion().GetIndex();
397 if (opposite) max = max+input->GetLargestPossibleRegion().GetSize();
399 while (!iter.IsAtEnd()) {
400 if (iter.Get() != bg) {
401 bool test = iter.GetIndex()[direction] > max[direction];
402 if (opposite) test = !test;
404 typename ImageType::PointType p;
405 input->TransformIndexToPhysicalPoint(iter.GetIndex(), p);
406 if ((distanceMax==0) || (p.EuclideanDistanceTo(refpoint) < distanceMax)) {
407 max = iter.GetIndex();
414 if (!found) return false;
415 input->TransformIndexToPhysicalPoint(max, point);
418 //--------------------------------------------------------------------
421 //--------------------------------------------------------------------
422 template<class ImageType>
423 typename ImageType::Pointer
424 CropImageRemoveGreaterThan(const ImageType * image,
425 int dim, double min, bool autoCrop,
426 typename ImageType::PixelType BG)
428 return CropImageAlongOneAxis<ImageType>(image, dim,
429 image->GetOrigin()[dim],
433 //--------------------------------------------------------------------
436 //--------------------------------------------------------------------
437 template<class ImageType>
438 typename ImageType::Pointer
439 CropImageRemoveLowerThan(const ImageType * image,
440 int dim, double max, bool autoCrop,
441 typename ImageType::PixelType BG)
443 typename ImageType::PointType p;
444 image->TransformIndexToPhysicalPoint(image->GetLargestPossibleRegion().GetIndex()+
445 image->GetLargestPossibleRegion().GetSize(), p);
446 return CropImageAlongOneAxis<ImageType>(image, dim, max, p[dim], autoCrop, BG);
448 //--------------------------------------------------------------------
451 //--------------------------------------------------------------------
452 template<class ImageType>
453 typename ImageType::Pointer
454 CropImageAlongOneAxis(const ImageType * image,
455 int dim, double min, double max,
456 bool autoCrop, typename ImageType::PixelType BG)
458 // Compute region size
459 typename ImageType::RegionType region;
460 typename ImageType::SizeType size = image->GetLargestPossibleRegion().GetSize();
461 typename ImageType::PointType p = image->GetOrigin();
463 typename ImageType::IndexType start;
464 image->TransformPhysicalPointToIndex(p, start);
466 typename ImageType::IndexType end;
467 image->TransformPhysicalPointToIndex(p, end);
468 size[dim] = fabs(end[dim]-start[dim]);
469 region.SetIndex(start);
470 region.SetSize(size);
473 typedef itk::RegionOfInterestImageFilter<ImageType, ImageType> CropFilterType;
474 typename CropFilterType::Pointer cropFilter = CropFilterType::New();
475 cropFilter->SetInput(image);
476 cropFilter->SetRegionOfInterest(region);
477 cropFilter->Update();
478 typename ImageType::Pointer result = cropFilter->GetOutput();
482 result = AutoCrop<ImageType>(result, BG);
486 //--------------------------------------------------------------------
489 //--------------------------------------------------------------------
490 template<class ImageType>
492 ComputeCentroids(const ImageType * image,
493 typename ImageType::PixelType BG,
494 std::vector<typename ImageType::PointType> & centroids)
496 typedef long LabelType;
497 static const unsigned int Dim = ImageType::ImageDimension;
498 typedef itk::ShapeLabelObject< LabelType, Dim > LabelObjectType;
499 typedef itk::LabelMap< LabelObjectType > LabelMapType;
500 typedef itk::LabelImageToLabelMapFilter<ImageType, LabelMapType> ImageToMapFilterType;
501 typename ImageToMapFilterType::Pointer imageToLabelFilter = ImageToMapFilterType::New();
502 typedef itk::ShapeLabelMapFilter<LabelMapType, ImageType> ShapeFilterType;
503 typename ShapeFilterType::Pointer statFilter = ShapeFilterType::New();
504 imageToLabelFilter->SetBackgroundValue(BG);
505 imageToLabelFilter->SetInput(image);
506 statFilter->SetInput(imageToLabelFilter->GetOutput());
507 statFilter->Update();
508 typename LabelMapType::Pointer labelMap = statFilter->GetOutput();
511 typename ImageType::PointType dummy;
512 centroids.push_back(dummy); // label 0 -> no centroid, use dummy point for BG
513 //DS FIXME (not useful ! to change ..)
514 DD(labelMap->GetNumberOfLabelObjects());
515 for(uint i=0; i<labelMap->GetNumberOfLabelObjects(); i++) {
516 int label = labelMap->GetLabels()[i];
517 centroids.push_back(labelMap->GetLabelObject(label)->GetCentroid());
520 //--------------------------------------------------------------------
523 //--------------------------------------------------------------------
524 template<class ImageType>
526 ComputeCentroids2(const ImageType * image,
527 typename ImageType::PixelType BG,
528 std::vector<typename ImageType::PointType> & centroids)
530 typedef long LabelType;
531 static const unsigned int Dim = ImageType::ImageDimension;
532 typedef itk::ShapeLabelObject< LabelType, Dim > LabelObjectType;
533 typedef itk::LabelMap< LabelObjectType > LabelMapType;
534 typedef itk::LabelImageToLabelMapFilter<ImageType, LabelMapType> ImageToMapFilterType;
535 typename ImageToMapFilterType::Pointer imageToLabelFilter = ImageToMapFilterType::New();
536 typedef itk::ShapeLabelMapFilter<LabelMapType, ImageType> ShapeFilterType;
537 typename ShapeFilterType::Pointer statFilter = ShapeFilterType::New();
538 imageToLabelFilter->SetBackgroundValue(BG);
539 imageToLabelFilter->SetInput(image);
540 statFilter->SetInput(imageToLabelFilter->GetOutput());
541 statFilter->Update();
542 typename LabelMapType::Pointer labelMap = statFilter->GetOutput();
545 typename ImageType::PointType dummy;
546 centroids.push_back(dummy); // label 0 -> no centroid, use dummy point
547 for(uint i=1; i<labelMap->GetNumberOfLabelObjects()+1; i++) {
548 centroids.push_back(labelMap->GetLabelObject(i)->GetCentroid());
551 for(uint i=1; i<labelMap->GetNumberOfLabelObjects()+1; i++) {
552 DD(labelMap->GetLabelObject(i)->GetBinaryPrincipalAxes());
553 DD(labelMap->GetLabelObject(i)->GetBinaryFlatness());
554 DD(labelMap->GetLabelObject(i)->GetRoundness ());
556 // search for the point on the boundary alog PA
561 //--------------------------------------------------------------------
564 //--------------------------------------------------------------------
565 template<class ImageType>
567 ExtractSlices(const ImageType * image, int direction,
568 std::vector<typename itk::Image<typename ImageType::PixelType,
569 ImageType::ImageDimension-1>::Pointer > & slices)
571 typedef ExtractSliceFilter<ImageType> ExtractSliceFilterType;
572 typedef typename ExtractSliceFilterType::SliceType SliceType;
573 typename ExtractSliceFilterType::Pointer
574 extractSliceFilter = ExtractSliceFilterType::New();
575 extractSliceFilter->SetInput(image);
576 extractSliceFilter->SetDirection(direction);
577 extractSliceFilter->Update();
578 extractSliceFilter->GetOutputSlices(slices);
580 //--------------------------------------------------------------------
583 //--------------------------------------------------------------------
584 template<class ImageType>
586 PointsUtils<ImageType>::Convert2DTo3D(const PointType2D & p2D,
587 const ImageType * image,
592 index3D[0] = index3D[1] = 0;
593 index3D[2] = image->GetLargestPossibleRegion().GetIndex()[2]+slice;
594 image->TransformIndexToPhysicalPoint(index3D, p3D);
597 // p3D[2] = p[2];//(image->GetLargestPossibleRegion().GetIndex()[2]+slice)*image->GetSpacing()[2]
598 // + image->GetOrigin()[2];
600 //--------------------------------------------------------------------
603 //--------------------------------------------------------------------
604 template<class ImageType>
606 PointsUtils<ImageType>::Convert2DMapTo3DList(const MapPoint2DType & map,
607 const ImageType * image,
608 VectorPoint3DType & list)
610 typename MapPoint2DType::const_iterator iter = map.begin();
611 while (iter != map.end()) {
613 Convert2DTo3D(iter->second, image, iter->first, p);
618 //--------------------------------------------------------------------
621 //--------------------------------------------------------------------
622 template<class ImageType>
624 PointsUtils<ImageType>::Convert2DListTo3DList(const VectorPoint2DType & p2D,
626 const ImageType * image,
627 VectorPoint3DType & list)
629 for(uint i=0; i<p2D.size(); i++) {
631 Convert2DTo3D(p2D[i], image, slice, p);
635 //--------------------------------------------------------------------
638 //--------------------------------------------------------------------
639 template<class ImageType>
641 WriteListOfLandmarks(std::vector<typename ImageType::PointType> points,
642 std::string filename)
645 openFileForWriting(os, filename);
646 os << "LANDMARKS1" << std::endl;
647 for(uint i=0; i<points.size(); i++) {
648 const typename ImageType::PointType & p = points[i];
649 // Write it in the file
650 os << i << " " << p[0] << " " << p[1] << " " << p[2] << " 0 0 " << std::endl;
654 //--------------------------------------------------------------------
657 //--------------------------------------------------------------------
658 template<class ImageType>
659 typename ImageType::Pointer
660 Dilate(const ImageType * image, double radiusInMM,
661 typename ImageType::PixelType BG,
662 typename ImageType::PixelType FG,
665 typename ImageType::SizeType r;
666 for(uint i=0; i<ImageType::ImageDimension; i++)
667 r[i] = (uint)lrint(radiusInMM/image->GetSpacing()[i]);
668 return Dilate<ImageType>(image, r, BG, FG, extendSupport);
670 //--------------------------------------------------------------------
673 //--------------------------------------------------------------------
674 template<class ImageType>
675 typename ImageType::Pointer
676 Dilate(const ImageType * image, typename ImageType::PointType radiusInMM,
677 typename ImageType::PixelType BG,
678 typename ImageType::PixelType FG,
681 typename ImageType::SizeType r;
682 for(uint i=0; i<ImageType::ImageDimension; i++)
683 r[i] = (uint)lrint(radiusInMM[i]/image->GetSpacing()[i]);
684 return Dilate<ImageType>(image, r, BG, FG, extendSupport);
686 //--------------------------------------------------------------------
689 //--------------------------------------------------------------------
690 template<class ImageType>
691 typename ImageType::Pointer
692 Dilate(const ImageType * image, typename ImageType::SizeType radius,
693 typename ImageType::PixelType BG,
694 typename ImageType::PixelType FG,
697 // Create kernel for dilatation
698 typedef itk::BinaryBallStructuringElement<typename ImageType::PixelType,
699 ImageType::ImageDimension> KernelType;
700 KernelType structuringElement;
701 structuringElement.SetRadius(radius);
702 structuringElement.CreateStructuringElement();
704 typename ImageType::Pointer output;
706 typedef itk::ConstantPadImageFilter<ImageType, ImageType> PadFilterType;
707 typename PadFilterType::Pointer padFilter = PadFilterType::New();
708 padFilter->SetInput(image);
709 typename ImageType::SizeType lower;
710 typename ImageType::SizeType upper;
711 for(uint i=0; i<3; i++) {
712 lower[i] = upper[i] = 2*(radius[i]+1);
714 padFilter->SetPadLowerBound(lower);
715 padFilter->SetPadUpperBound(upper);
717 output = padFilter->GetOutput();
721 typedef itk::BinaryDilateImageFilter<ImageType, ImageType , KernelType> DilateFilterType;
722 typename DilateFilterType::Pointer dilateFilter = DilateFilterType::New();
723 dilateFilter->SetBackgroundValue(BG);
724 dilateFilter->SetForegroundValue(FG);
725 dilateFilter->SetBoundaryToForeground(false);
726 dilateFilter->SetKernel(structuringElement);
727 dilateFilter->SetInput(output);
728 dilateFilter->Update();
729 return dilateFilter->GetOutput();
731 //--------------------------------------------------------------------
734 //--------------------------------------------------------------------
735 template<class ImageType>
736 typename ImageType::Pointer
737 Opening(const ImageType * image, typename ImageType::SizeType radius,
738 typename ImageType::PixelType BG,
739 typename ImageType::PixelType FG)
742 typedef itk::BinaryBallStructuringElement<typename ImageType::PixelType,
743 ImageType::ImageDimension> KernelType;
744 KernelType structuringElement;
745 structuringElement.SetRadius(radius);
746 structuringElement.CreateStructuringElement();
749 typedef itk::BinaryMorphologicalOpeningImageFilter<ImageType, ImageType , KernelType> OpeningFilterType;
750 typename OpeningFilterType::Pointer open = OpeningFilterType::New();
751 open->SetInput(image);
752 open->SetBackgroundValue(BG);
753 open->SetForegroundValue(FG);
754 open->SetKernel(structuringElement);
756 return open->GetOutput();
758 //--------------------------------------------------------------------
762 //--------------------------------------------------------------------
763 template<class ValueType, class VectorType>
764 void ConvertOption(std::string optionName, uint given,
765 ValueType * values, VectorType & p,
766 uint dim, bool required)
768 if (required && (given == 0)) {
769 clitkExceptionMacro("The option --" << optionName << " must be set and have 1 or "
770 << dim << " values.");
773 for(uint i=0; i<dim; i++) p[i] = values[0];
777 for(uint i=0; i<dim; i++) p[i] = values[i];
780 if (given == 0) return;
781 clitkExceptionMacro("The option --" << optionName << " must have 1 or "
782 << dim << " values.");
784 //--------------------------------------------------------------------
787 //--------------------------------------------------------------------
789 http://www.gamedev.net/community/forums/topic.asp?topic_id=542870
790 Assuming the points are (Ax,Ay) (Bx,By) and (Cx,Cy), you need to compute:
791 (Bx - Ax) * (Cy - Ay) - (By - Ay) * (Cx - Ax)
792 This will equal zero if the point C is on the line formed by
793 points A and B, and will have a different sign depending on the
794 side. Which side this is depends on the orientation of your (x,y)
795 coordinates, but you can plug test values for A,B and C into this
796 formula to determine whether negative values are to the left or to
798 => to accelerate, start with formula, when change sign -> stop and fill
800 offsetToKeep = is used to determine which side of the line we
801 keep. The point along the mainDirection but 'offsetToKeep' mm away
805 template<class ImageType>
807 SliceBySliceSetBackgroundFromLineSeparation(ImageType * input,
808 std::vector<typename ImageType::PointType> & lA,
809 std::vector<typename ImageType::PointType> & lB,
810 typename ImageType::PixelType BG,
814 typedef itk::ImageSliceIteratorWithIndex<ImageType> SliceIteratorType;
815 SliceIteratorType siter = SliceIteratorType(input,
816 input->GetLargestPossibleRegion());
817 siter.SetFirstDirection(0);
818 siter.SetSecondDirection(1);
821 typename ImageType::PointType A;
822 typename ImageType::PointType B;
823 typename ImageType::PointType C;
824 assert(lA.size() == lB.size());
825 while ((i<lA.size()) && (!siter.IsAtEnd())) {
826 // Check that the current slice correspond to the current point
827 input->TransformIndexToPhysicalPoint(siter.GetIndex(), C);
831 if ((fabs(C[2] - lA[i][2]))>0.01) { // is !equal with a tolerance of 0.01 mm
834 // Define A,B,C points
842 // Check that the line is not a point (A=B)
843 bool p = (A[0] == B[0]) && (A[1] == B[1]);
846 C[mainDirection] += offsetToKeep; // I know I must keep this point
847 double s = (B[0] - A[0]) * (C[1] - A[1]) - (B[1] - A[1]) * (C[0] - A[0]);
848 bool isPositive = s<0;
849 while (!siter.IsAtEndOfSlice()) {
850 while (!siter.IsAtEndOfLine()) {
851 // Very slow, I know ... but image should be very small
852 input->TransformIndexToPhysicalPoint(siter.GetIndex(), C);
853 double s = (B[0] - A[0]) * (C[1] - A[1]) - (B[1] - A[1]) * (C[0] - A[0]);
854 if (s == 0) siter.Set(BG); // on the line, we decide to remove
856 if (s > 0) siter.Set(BG);
859 if (s < 0) siter.Set(BG);
868 } // End of current slice
872 //--------------------------------------------------------------------
874 //--------------------------------------------------------------------
875 template<class ImageType>
877 AndNot(ImageType * input,
878 const ImageType * object,
879 typename ImageType::PixelType BG)
881 typename ImageType::Pointer o;
883 if (!clitk::HaveSameSizeAndSpacing<ImageType, ImageType>(input, object)) {
884 o = clitk::ResizeImageLike<ImageType>(object, input, BG);
888 typedef clitk::BooleanOperatorLabelImageFilter<ImageType> BoolFilterType;
889 typename BoolFilterType::Pointer boolFilter = BoolFilterType::New();
890 boolFilter->InPlaceOn();
891 boolFilter->SetInput1(input);
892 if (resized) boolFilter->SetInput2(o);
893 else boolFilter->SetInput2(object);
894 boolFilter->SetBackgroundValue1(BG);
895 boolFilter->SetBackgroundValue2(BG);
896 boolFilter->SetOperationType(BoolFilterType::AndNot);
897 boolFilter->Update();
899 //--------------------------------------------------------------------
902 //--------------------------------------------------------------------
903 template<class ImageType>
904 typename ImageType::Pointer
905 Binarize(const ImageType * input,
906 typename ImageType::PixelType lower,
907 typename ImageType::PixelType upper,
908 typename ImageType::PixelType BG,
909 typename ImageType::PixelType FG)
911 typedef itk::BinaryThresholdImageFilter<ImageType, ImageType> BinaryThresholdFilterType;
912 typename BinaryThresholdFilterType::Pointer binarizeFilter = BinaryThresholdFilterType::New();
913 binarizeFilter->SetInput(input);
914 binarizeFilter->InPlaceOff();
915 binarizeFilter->SetLowerThreshold(lower);
916 binarizeFilter->SetUpperThreshold(upper);
917 binarizeFilter->SetInsideValue(FG);
918 binarizeFilter->SetOutsideValue(BG);
919 binarizeFilter->Update();
920 return binarizeFilter->GetOutput();
922 //--------------------------------------------------------------------
925 //--------------------------------------------------------------------
926 template<class ImageType>
928 GetMinMaxPointPosition(const ImageType * input,
929 typename ImageType::PointType & min,
930 typename ImageType::PointType & max)
932 typename ImageType::IndexType index = input->GetLargestPossibleRegion().GetIndex();
933 input->TransformIndexToPhysicalPoint(index, min);
934 index = index+input->GetLargestPossibleRegion().GetSize();
935 input->TransformIndexToPhysicalPoint(index, max);
937 //--------------------------------------------------------------------
940 //--------------------------------------------------------------------
941 template<class ImageType>
942 typename ImageType::PointType
943 FindExtremaPointInAGivenLine(const ImageType * input,
946 typename ImageType::PointType p,
947 typename ImageType::PixelType BG,
950 // Which direction ? Increasing or decreasing.
954 // Transform to pixel index
955 typename ImageType::IndexType index;
956 input->TransformPhysicalPointToIndex(p, index);
958 // Loop while inside the mask;
959 while (input->GetPixel(index) != BG) {
960 index[dimension] += d;
963 // Transform back to Physical Units
964 typename ImageType::PointType result;
965 input->TransformIndexToPhysicalPoint(index, result);
967 // Check that is is not too far away
968 double distance = p.EuclideanDistanceTo(result);
969 if (distance > distanceMax) {
970 result = p; // Get back to initial value
975 //--------------------------------------------------------------------
978 //--------------------------------------------------------------------
979 template<class PointType>
981 IsOnTheSameLineSide(PointType C, PointType A, PointType B, PointType like)
983 // Look at the position of point 'like' according to the AB line
984 double s = (B[0] - A[0]) * (like[1] - A[1]) - (B[1] - A[1]) * (like[0] - A[0]);
987 // Look the C position
988 s = (B[0] - A[0]) * (C[1] - A[1]) - (B[1] - A[1]) * (C[0] - A[0]);
990 if (negative && (s<=0)) return true;
991 if (!negative && (s>=0)) return true;
994 //--------------------------------------------------------------------
997 //--------------------------------------------------------------------
998 /* Consider an input object, for each slice, find the extrema
999 position according to a given direction and build a line segment
1000 passing throught this point in a given direction. Output is a
1001 vector of line (from point A to B), for each slice;
1003 template<class ImageType>
1005 SliceBySliceBuildLineSegmentAccordingToExtremaPosition(const ImageType * input,
1006 typename ImageType::PixelType BG,
1008 int extremaDirection,
1009 bool extremaOppositeFlag,
1012 std::vector<typename ImageType::PointType> & A,
1013 std::vector<typename ImageType::PointType> & B)
1016 typedef typename itk::Image<typename ImageType::PixelType, ImageType::ImageDimension-1> SliceType;
1018 // Build the list of slices
1019 std::vector<typename SliceType::Pointer> slices;
1020 clitk::ExtractSlices<ImageType>(input, sliceDimension, slices);
1022 // Build the list of 2D points
1023 std::map<int, typename SliceType::PointType> position2D;
1024 for(uint i=0; i<slices.size(); i++) {
1025 typename SliceType::PointType p;
1027 clitk::FindExtremaPointInAGivenDirection<SliceType>(slices[i], BG,
1028 extremaDirection, extremaOppositeFlag, p);
1034 // Convert 2D points in slice into 3D points
1035 clitk::PointsUtils<ImageType>::Convert2DMapTo3DList(position2D, input, A);
1037 // Create additional point just right to the previous ones, on the
1038 // given lineDirection, in order to create a horizontal/vertical line.
1039 for(uint i=0; i<A.size(); i++) {
1040 typename ImageType::PointType p = A[i];
1041 p[lineDirection] += 10;
1049 //--------------------------------------------------------------------
1052 //--------------------------------------------------------------------
1053 template<class ImageType>
1054 typename ImageType::Pointer
1055 SliceBySliceKeepMainCCL(const ImageType * input,
1056 typename ImageType::PixelType BG,
1057 typename ImageType::PixelType FG) {
1060 const int d = ImageType::ImageDimension-1;
1061 typedef typename itk::Image<typename ImageType::PixelType, d> SliceType;
1062 std::vector<typename SliceType::Pointer> slices;
1063 clitk::ExtractSlices<ImageType>(input, d, slices);
1066 // Labelize and keep the main one
1067 std::vector<typename SliceType::Pointer> o;
1068 for(uint i=0; i<slices.size(); i++) {
1070 o.push_back(clitk::Labelize<SliceType>(slices[i], BG, false, 1));
1071 o[i] = clitk::KeepLabels<SliceType>(o[i], BG, FG, 1, 1, true);
1076 typename ImageType::Pointer output;
1077 output = clitk::JoinSlices<ImageType>(o, input, d);
1081 //--------------------------------------------------------------------
1084 } // end of namespace