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>
38 //--------------------------------------------------------------------
39 template<class ImageType>
40 void ComputeBBFromImageRegion(const ImageType * image,
41 typename ImageType::RegionType region,
42 typename itk::BoundingBox<unsigned long,
43 ImageType::ImageDimension>::Pointer bb) {
44 typedef typename ImageType::IndexType IndexType;
47 for(unsigned int i=0; i<image->GetImageDimension(); i++) {
48 firstIndex[i] = region.GetIndex()[i];
49 lastIndex[i] = firstIndex[i]+region.GetSize()[i];
52 typedef itk::BoundingBox<unsigned long,
53 ImageType::ImageDimension> BBType;
54 typedef typename BBType::PointType PointType;
57 image->TransformIndexToPhysicalPoint(firstIndex, firstPoint);
58 image->TransformIndexToPhysicalPoint(lastIndex, lastPoint);
60 bb->SetMaximum(lastPoint);
61 bb->SetMinimum(firstPoint);
63 //--------------------------------------------------------------------
66 //--------------------------------------------------------------------
67 template<int Dimension>
68 void ComputeBBIntersection(typename itk::BoundingBox<unsigned long, Dimension>::Pointer bbo,
69 typename itk::BoundingBox<unsigned long, Dimension>::Pointer bbi1,
70 typename itk::BoundingBox<unsigned long, Dimension>::Pointer bbi2) {
72 typedef itk::BoundingBox<unsigned long, Dimension> BBType;
73 typedef typename BBType::PointType PointType;
77 for(unsigned int i=0; i<Dimension; i++) {
78 firstPoint[i] = std::max(bbi1->GetMinimum()[i],
79 bbi2->GetMinimum()[i]);
80 lastPoint[i] = std::min(bbi1->GetMaximum()[i],
81 bbi2->GetMaximum()[i]);
84 bbo->SetMaximum(lastPoint);
85 bbo->SetMinimum(firstPoint);
87 //--------------------------------------------------------------------
90 //--------------------------------------------------------------------
91 template<class ImageType>
92 void ComputeRegionFromBB(const ImageType * image,
93 const typename itk::BoundingBox<unsigned long,
94 ImageType::ImageDimension>::Pointer bb,
95 typename ImageType::RegionType & region) {
97 typedef typename ImageType::IndexType IndexType;
98 typedef typename ImageType::PointType PointType;
99 typedef typename ImageType::RegionType RegionType;
100 typedef typename ImageType::SizeType SizeType;
102 // Region starting point
103 IndexType regionStart;
104 PointType start = bb->GetMinimum();
105 image->TransformPhysicalPointToIndex(start, regionStart);
109 PointType maxs = bb->GetMaximum();
110 PointType mins = bb->GetMinimum();
111 for(unsigned int i=0; i<ImageType::ImageDimension; i++) {
112 regionSize[i] = lrint((maxs[i] - mins[i])/image->GetSpacing()[i]);
116 region.SetIndex(regionStart);
117 region.SetSize(regionSize);
119 //--------------------------------------------------------------------
121 //--------------------------------------------------------------------
122 template<class ImageType, class TMaskImageType>
123 typename ImageType::Pointer
124 SetBackground(const ImageType * input,
125 const TMaskImageType * mask,
126 typename TMaskImageType::PixelType maskBG,
127 typename ImageType::PixelType outValue,
129 typedef SetBackgroundImageFilter<ImageType, TMaskImageType, ImageType>
130 SetBackgroundImageFilterType;
131 typename SetBackgroundImageFilterType::Pointer setBackgroundFilter
132 = SetBackgroundImageFilterType::New();
133 // if (inPlace) setBackgroundFilter->ReleaseDataFlagOn(); // No seg fault
134 setBackgroundFilter->SetInPlace(inPlace); // This is important to keep memory low
135 setBackgroundFilter->SetInput(input);
136 setBackgroundFilter->SetInput2(mask);
137 setBackgroundFilter->SetMaskValue(maskBG);
138 setBackgroundFilter->SetOutsideValue(outValue);
139 setBackgroundFilter->Update();
140 return setBackgroundFilter->GetOutput();
142 //--------------------------------------------------------------------
145 //--------------------------------------------------------------------
146 template<class ImageType>
147 int GetNumberOfConnectedComponentLabels(const ImageType * input,
148 typename ImageType::PixelType BG,
149 bool isFullyConnected) {
150 // Connected Component label
151 typedef itk::ConnectedComponentImageFilter<ImageType, ImageType> ConnectFilterType;
152 typename ConnectFilterType::Pointer connectFilter = ConnectFilterType::New();
153 connectFilter->SetInput(input);
154 connectFilter->SetBackgroundValue(BG);
155 connectFilter->SetFullyConnected(isFullyConnected);
156 connectFilter->Update();
159 return connectFilter->GetObjectCount();
161 //--------------------------------------------------------------------
163 //--------------------------------------------------------------------
165 Warning : in this cas, we consider outputType like inputType, not
166 InternalImageType. Be sure it fits.
168 template<class ImageType>
169 typename ImageType::Pointer
170 Labelize(const ImageType * input,
171 typename ImageType::PixelType BG,
172 bool isFullyConnected,
173 int minimalComponentSize) {
174 // InternalImageType for storing large number of component
175 typedef itk::Image<int, ImageType::ImageDimension> InternalImageType;
177 // Connected Component label
178 typedef itk::ConnectedComponentImageFilter<ImageType, InternalImageType> ConnectFilterType;
179 typename ConnectFilterType::Pointer connectFilter = ConnectFilterType::New();
180 // connectFilter->ReleaseDataFlagOn();
181 connectFilter->SetInput(input);
182 connectFilter->SetBackgroundValue(BG);
183 connectFilter->SetFullyConnected(isFullyConnected);
185 // Sort by size and remove too small area.
186 typedef itk::RelabelComponentImageFilter<InternalImageType, ImageType> RelabelFilterType;
187 typename RelabelFilterType::Pointer relabelFilter = RelabelFilterType::New();
188 // relabelFilter->ReleaseDataFlagOn(); // if yes, fail when ExplosionControlledThresholdConnectedImageFilter ???
189 relabelFilter->SetInput(connectFilter->GetOutput());
190 relabelFilter->SetMinimumObjectSize(minimalComponentSize);
191 relabelFilter->Update();
194 typename ImageType::Pointer output = relabelFilter->GetOutput();
197 //--------------------------------------------------------------------
200 //--------------------------------------------------------------------
202 Warning : in this cas, we consider outputType like inputType, not
203 InternalImageType. Be sure it fits.
205 template<class ImageType>
206 typename ImageType::Pointer
207 LabelizeAndCountNumberOfObjects(const ImageType * input,
208 typename ImageType::PixelType BG,
209 bool isFullyConnected,
210 int minimalComponentSize,
212 // InternalImageType for storing large number of component
213 typedef itk::Image<int, ImageType::ImageDimension> InternalImageType;
215 // Connected Component label
216 typedef itk::ConnectedComponentImageFilter<ImageType, InternalImageType> ConnectFilterType;
217 typename ConnectFilterType::Pointer connectFilter = ConnectFilterType::New();
218 // connectFilter->ReleaseDataFlagOn();
219 connectFilter->SetInput(input);
220 connectFilter->SetBackgroundValue(BG);
221 connectFilter->SetFullyConnected(isFullyConnected);
223 // Sort by size and remove too small area.
224 typedef itk::RelabelComponentImageFilter<InternalImageType, ImageType> RelabelFilterType;
225 typename RelabelFilterType::Pointer relabelFilter = RelabelFilterType::New();
226 // relabelFilter->ReleaseDataFlagOn(); // if yes, fail when ExplosionControlledThresholdConnectedImageFilter ???
227 relabelFilter->SetInput(connectFilter->GetOutput());
228 relabelFilter->SetMinimumObjectSize(minimalComponentSize);
229 relabelFilter->Update();
231 nb = relabelFilter->GetNumberOfObjects();
232 // DD(relabelFilter->GetOriginalNumberOfObjects());
233 // DD(relabelFilter->GetSizeOfObjectsInPhysicalUnits()[0]);
236 typename ImageType::Pointer output = relabelFilter->GetOutput();
239 //--------------------------------------------------------------------
243 //--------------------------------------------------------------------
244 template<class ImageType>
245 typename ImageType::Pointer
246 RemoveLabels(const ImageType * input,
247 typename ImageType::PixelType BG,
248 std::vector<typename ImageType::PixelType> & labelsToRemove) {
249 assert(labelsToRemove.size() != 0);
250 typename ImageType::Pointer working_image;// = input;
251 for (unsigned int i=0; i <labelsToRemove.size(); i++) {
252 typedef SetBackgroundImageFilter<ImageType, ImageType> SetBackgroundImageFilterType;
253 typename SetBackgroundImageFilterType::Pointer setBackgroundFilter = SetBackgroundImageFilterType::New();
254 setBackgroundFilter->SetInput(input);
255 setBackgroundFilter->SetInput2(input);
256 setBackgroundFilter->SetMaskValue(labelsToRemove[i]);
257 setBackgroundFilter->SetOutsideValue(BG);
258 setBackgroundFilter->Update();
259 working_image = setBackgroundFilter->GetOutput();
261 return working_image;
263 //--------------------------------------------------------------------
266 //--------------------------------------------------------------------
267 template<class ImageType>
268 typename ImageType::Pointer
269 KeepLabels(const ImageType * input,
270 typename ImageType::PixelType BG,
271 typename ImageType::PixelType FG,
272 typename ImageType::PixelType firstKeep,
273 typename ImageType::PixelType lastKeep,
275 typedef itk::BinaryThresholdImageFilter<ImageType, ImageType> BinarizeFilterType;
276 typename BinarizeFilterType::Pointer binarizeFilter = BinarizeFilterType::New();
277 binarizeFilter->SetInput(input);
278 binarizeFilter->SetLowerThreshold(firstKeep);
279 if (useLastKeep) binarizeFilter->SetUpperThreshold(lastKeep);
280 binarizeFilter->SetInsideValue(FG);
281 binarizeFilter->SetOutsideValue(BG);
282 binarizeFilter->Update();
283 return binarizeFilter->GetOutput();
285 //--------------------------------------------------------------------
288 //--------------------------------------------------------------------
289 template<class ImageType>
290 typename ImageType::Pointer
291 LabelizeAndSelectLabels(const ImageType * input,
292 typename ImageType::PixelType BG,
293 typename ImageType::PixelType FG,
294 bool isFullyConnected,
295 int minimalComponentSize,
296 LabelizeParameters<typename ImageType::PixelType> * param)
298 typename ImageType::Pointer working_image;
299 working_image = Labelize<ImageType>(input, BG, isFullyConnected, minimalComponentSize);
300 if (param->GetLabelsToRemove().size() != 0)
301 working_image = RemoveLabels<ImageType>(working_image, BG, param->GetLabelsToRemove());
302 working_image = KeepLabels<ImageType>(working_image,
304 param->GetFirstKeep(),
305 param->GetLastKeep(),
306 param->GetUseLastKeep());
307 return working_image;
309 //--------------------------------------------------------------------
312 //--------------------------------------------------------------------
313 template<class ImageType>
314 typename ImageType::Pointer
315 ResizeImageLike(const ImageType * input,
316 const itk::ImageBase<ImageType::ImageDimension> * like,
317 typename ImageType::PixelType backgroundValue)
319 typedef CropLikeImageFilter<ImageType> CropFilterType;
320 typename CropFilterType::Pointer cropFilter = CropFilterType::New();
321 cropFilter->SetInput(input);
322 cropFilter->SetCropLikeImage(like);
323 cropFilter->SetBackgroundValue(backgroundValue);
324 cropFilter->Update();
325 return cropFilter->GetOutput();
327 //--------------------------------------------------------------------
330 //--------------------------------------------------------------------
331 template<class MaskImageType>
332 typename MaskImageType::Pointer
333 SliceBySliceRelativePosition(const MaskImageType * input,
334 const MaskImageType * object,
337 std::string orientation,
338 bool uniqueConnectedComponent,
341 bool singleObjectCCL)
343 typedef SliceBySliceRelativePositionFilter<MaskImageType> SliceRelPosFilterType;
344 typename SliceRelPosFilterType::Pointer sliceRelPosFilter = SliceRelPosFilterType::New();
345 sliceRelPosFilter->VerboseStepFlagOff();
346 sliceRelPosFilter->WriteStepFlagOff();
347 sliceRelPosFilter->SetInput(input);
348 sliceRelPosFilter->SetInputObject(object);
349 sliceRelPosFilter->SetDirection(direction);
350 sliceRelPosFilter->SetFuzzyThreshold(threshold);
351 sliceRelPosFilter->AddOrientationTypeString(orientation);
352 sliceRelPosFilter->SetIntermediateSpacingFlag((spacing != -1));
353 sliceRelPosFilter->SetIntermediateSpacing(spacing);
354 sliceRelPosFilter->SetUniqueConnectedComponentBySlice(uniqueConnectedComponent);
355 sliceRelPosFilter->SetUseASingleObjectConnectedComponentBySliceFlag(singleObjectCCL);
356 // sliceRelPosFilter->SetInverseOrientationFlag(inverseflag);
357 sliceRelPosFilter->SetAutoCropFlag(autocropFlag);
358 sliceRelPosFilter->IgnoreEmptySliceObjectFlagOn();
359 sliceRelPosFilter->Update();
360 return sliceRelPosFilter->GetOutput();
362 //--------------------------------------------------------------------
364 //--------------------------------------------------------------------
365 template<class ImageType>
367 FindExtremaPointInAGivenDirection(const ImageType * input,
368 typename ImageType::PixelType bg,
369 int direction, bool opposite,
370 typename ImageType::PointType & point)
372 typename ImageType::PointType dummy;
373 return FindExtremaPointInAGivenDirection(input, bg, direction,
374 opposite, dummy, 0, point);
376 //--------------------------------------------------------------------
378 //--------------------------------------------------------------------
379 template<class ImageType>
381 FindExtremaPointInAGivenDirection(const ImageType * input,
382 typename ImageType::PixelType bg,
383 int direction, bool opposite,
384 typename ImageType::PointType refpoint,
386 typename ImageType::PointType & point)
389 loop over input pixels, store the index in the fg that is max
390 according to the given direction.
392 typedef itk::ImageRegionConstIteratorWithIndex<ImageType> IteratorType;
393 IteratorType iter(input, input->GetLargestPossibleRegion());
395 typename ImageType::IndexType max = input->GetLargestPossibleRegion().GetIndex();
396 if (opposite) max = max+input->GetLargestPossibleRegion().GetSize();
398 while (!iter.IsAtEnd()) {
399 if (iter.Get() != bg) {
400 bool test = iter.GetIndex()[direction] > max[direction];
401 if (opposite) test = !test;
403 typename ImageType::PointType p;
404 input->TransformIndexToPhysicalPoint(iter.GetIndex(), p);
405 if ((distanceMax==0) || (p.EuclideanDistanceTo(refpoint) < distanceMax)) {
406 max = iter.GetIndex();
413 if (!found) return false;
414 input->TransformIndexToPhysicalPoint(max, point);
417 //--------------------------------------------------------------------
420 //--------------------------------------------------------------------
421 template<class ImageType>
422 typename ImageType::Pointer
423 CropImageAbove(const ImageType * image,
424 int dim, double min, bool autoCrop,
425 typename ImageType::PixelType BG)
427 return CropImageAlongOneAxis<ImageType>(image, dim,
428 image->GetOrigin()[dim],
432 //--------------------------------------------------------------------
435 //--------------------------------------------------------------------
436 template<class ImageType>
437 typename ImageType::Pointer
438 CropImageBelow(const ImageType * image,
439 int dim, double max, bool autoCrop,
440 typename ImageType::PixelType BG)
442 typename ImageType::PointType p;
443 image->TransformIndexToPhysicalPoint(image->GetLargestPossibleRegion().GetIndex()+
444 image->GetLargestPossibleRegion().GetSize(), p);
445 return CropImageAlongOneAxis<ImageType>(image, dim, max, p[dim], autoCrop, BG);
447 //--------------------------------------------------------------------
450 //--------------------------------------------------------------------
451 template<class ImageType>
452 typename ImageType::Pointer
453 CropImageAlongOneAxis(const ImageType * image,
454 int dim, double min, double max,
455 bool autoCrop, typename ImageType::PixelType BG)
457 // Compute region size
458 typename ImageType::RegionType region;
459 typename ImageType::SizeType size = image->GetLargestPossibleRegion().GetSize();
460 typename ImageType::PointType p = image->GetOrigin();
462 typename ImageType::IndexType start;
463 image->TransformPhysicalPointToIndex(p, start);
465 typename ImageType::IndexType end;
466 image->TransformPhysicalPointToIndex(p, end);
467 size[dim] = fabs(end[dim]-start[dim]);
468 region.SetIndex(start);
469 region.SetSize(size);
472 typedef itk::RegionOfInterestImageFilter<ImageType, ImageType> CropFilterType;
473 typename CropFilterType::Pointer cropFilter = CropFilterType::New();
474 cropFilter->SetInput(image);
475 cropFilter->SetRegionOfInterest(region);
476 cropFilter->Update();
477 typename ImageType::Pointer result = cropFilter->GetOutput();
481 result = AutoCrop<ImageType>(result, BG);
485 //--------------------------------------------------------------------
488 //--------------------------------------------------------------------
489 template<class ImageType>
491 ComputeCentroids(const ImageType * image,
492 typename ImageType::PixelType BG,
493 std::vector<typename ImageType::PointType> & centroids)
495 typedef long LabelType;
496 static const unsigned int Dim = ImageType::ImageDimension;
497 typedef itk::ShapeLabelObject< LabelType, Dim > LabelObjectType;
498 typedef itk::LabelMap< LabelObjectType > LabelMapType;
499 typedef itk::LabelImageToLabelMapFilter<ImageType, LabelMapType> ImageToMapFilterType;
500 typename ImageToMapFilterType::Pointer imageToLabelFilter = ImageToMapFilterType::New();
501 typedef itk::ShapeLabelMapFilter<LabelMapType, ImageType> ShapeFilterType;
502 typename ShapeFilterType::Pointer statFilter = ShapeFilterType::New();
503 imageToLabelFilter->SetBackgroundValue(BG);
504 imageToLabelFilter->SetInput(image);
505 statFilter->SetInput(imageToLabelFilter->GetOutput());
506 statFilter->Update();
507 typename LabelMapType::Pointer labelMap = statFilter->GetOutput();
510 typename ImageType::PointType dummy;
511 centroids.push_back(dummy); // label 0 -> no centroid, use dummy point
512 for(uint i=1; i<labelMap->GetNumberOfLabelObjects()+1; i++) {
513 centroids.push_back(labelMap->GetLabelObject(i)->GetCentroid());
516 //--------------------------------------------------------------------
519 //--------------------------------------------------------------------
520 template<class ImageType>
522 ComputeCentroids2(const ImageType * image,
523 typename ImageType::PixelType BG,
524 std::vector<typename ImageType::PointType> & centroids)
526 typedef long LabelType;
527 static const unsigned int Dim = ImageType::ImageDimension;
528 typedef itk::ShapeLabelObject< LabelType, Dim > LabelObjectType;
529 typedef itk::LabelMap< LabelObjectType > LabelMapType;
530 typedef itk::LabelImageToLabelMapFilter<ImageType, LabelMapType> ImageToMapFilterType;
531 typename ImageToMapFilterType::Pointer imageToLabelFilter = ImageToMapFilterType::New();
532 typedef itk::ShapeLabelMapFilter<LabelMapType, ImageType> ShapeFilterType;
533 typename ShapeFilterType::Pointer statFilter = ShapeFilterType::New();
534 imageToLabelFilter->SetBackgroundValue(BG);
535 imageToLabelFilter->SetInput(image);
536 statFilter->SetInput(imageToLabelFilter->GetOutput());
537 statFilter->Update();
538 typename LabelMapType::Pointer labelMap = statFilter->GetOutput();
541 typename ImageType::PointType dummy;
542 centroids.push_back(dummy); // label 0 -> no centroid, use dummy point
543 for(uint i=1; i<labelMap->GetNumberOfLabelObjects()+1; i++) {
544 centroids.push_back(labelMap->GetLabelObject(i)->GetCentroid());
547 for(uint i=1; i<labelMap->GetNumberOfLabelObjects()+1; i++) {
548 DD(labelMap->GetLabelObject(i)->GetBinaryPrincipalAxes());
549 DD(labelMap->GetLabelObject(i)->GetBinaryFlatness());
550 DD(labelMap->GetLabelObject(i)->GetRoundness ());
552 // search for the point on the boundary alog PA
557 //--------------------------------------------------------------------
560 //--------------------------------------------------------------------
561 template<class ImageType>
563 ExtractSlices(const ImageType * image, int direction,
564 std::vector<typename itk::Image<typename ImageType::PixelType,
565 ImageType::ImageDimension-1>::Pointer > & slices)
567 typedef ExtractSliceFilter<ImageType> ExtractSliceFilterType;
568 typedef typename ExtractSliceFilterType::SliceType SliceType;
569 typename ExtractSliceFilterType::Pointer
570 extractSliceFilter = ExtractSliceFilterType::New();
571 extractSliceFilter->SetInput(image);
572 extractSliceFilter->SetDirection(direction);
573 extractSliceFilter->Update();
574 extractSliceFilter->GetOutputSlices(slices);
576 //--------------------------------------------------------------------
579 //--------------------------------------------------------------------
580 template<class ImageType>
582 PointsUtils<ImageType>::Convert2DTo3D(const PointType2D & p2D,
583 const ImageType * image,
588 index3D[0] = index3D[1] = 0;
589 index3D[2] = image->GetLargestPossibleRegion().GetIndex()[2]+slice;
590 image->TransformIndexToPhysicalPoint(index3D, p3D);
593 // p3D[2] = p[2];//(image->GetLargestPossibleRegion().GetIndex()[2]+slice)*image->GetSpacing()[2]
594 // + image->GetOrigin()[2];
596 //--------------------------------------------------------------------
599 //--------------------------------------------------------------------
600 template<class ImageType>
602 PointsUtils<ImageType>::Convert2DMapTo3DList(const MapPoint2DType & map,
603 const ImageType * image,
604 VectorPoint3DType & list)
606 typename MapPoint2DType::const_iterator iter = map.begin();
607 while (iter != map.end()) {
609 Convert2DTo3D(iter->second, image, iter->first, p);
614 //--------------------------------------------------------------------
617 //--------------------------------------------------------------------
618 template<class ImageType>
620 PointsUtils<ImageType>::Convert2DListTo3DList(const VectorPoint2DType & p2D,
622 const ImageType * image,
623 VectorPoint3DType & list)
625 for(uint i=0; i<p2D.size(); i++) {
627 Convert2DTo3D(p2D[i], image, slice, p);
631 //--------------------------------------------------------------------
634 //--------------------------------------------------------------------
635 template<class ImageType>
637 WriteListOfLandmarks(std::vector<typename ImageType::PointType> points,
638 std::string filename)
641 openFileForWriting(os, filename);
642 os << "LANDMARKS1" << std::endl;
643 for(uint i=0; i<points.size(); i++) {
644 const typename ImageType::PointType & p = points[i];
645 // Write it in the file
646 os << i << " " << p[0] << " " << p[1] << " " << p[2] << " 0 0 " << std::endl;
650 //--------------------------------------------------------------------
653 //--------------------------------------------------------------------
654 template<class ImageType>
655 typename ImageType::Pointer
656 Dilate(const ImageType * image, double radiusInMM,
657 typename ImageType::PixelType BG,
658 typename ImageType::PixelType FG,
661 typename ImageType::SizeType r;
662 for(uint i=0; i<ImageType::ImageDimension; i++)
663 r[i] = (uint)lrint(radiusInMM/image->GetSpacing()[i]);
664 return Dilate<ImageType>(image, r, BG, FG, extendSupport);
666 //--------------------------------------------------------------------
669 //--------------------------------------------------------------------
670 template<class ImageType>
671 typename ImageType::Pointer
672 Dilate(const ImageType * image, typename ImageType::PointType radiusInMM,
673 typename ImageType::PixelType BG,
674 typename ImageType::PixelType FG,
677 typename ImageType::SizeType r;
678 for(uint i=0; i<ImageType::ImageDimension; i++)
679 r[i] = (uint)lrint(radiusInMM[i]/image->GetSpacing()[i]);
680 return Dilate<ImageType>(image, r, BG, FG, extendSupport);
682 //--------------------------------------------------------------------
685 //--------------------------------------------------------------------
686 template<class ImageType>
687 typename ImageType::Pointer
688 Dilate(const ImageType * image, typename ImageType::SizeType radius,
689 typename ImageType::PixelType BG,
690 typename ImageType::PixelType FG,
693 // Create kernel for dilatation
694 typedef itk::BinaryBallStructuringElement<typename ImageType::PixelType,
695 ImageType::ImageDimension> KernelType;
696 KernelType structuringElement;
697 structuringElement.SetRadius(radius);
698 structuringElement.CreateStructuringElement();
700 typename ImageType::Pointer output;
702 typedef itk::ConstantPadImageFilter<ImageType, ImageType> PadFilterType;
703 typename PadFilterType::Pointer padFilter = PadFilterType::New();
704 padFilter->SetInput(image);
705 typename ImageType::SizeType lower;
706 typename ImageType::SizeType upper;
707 for(uint i=0; i<3; i++) {
708 lower[i] = upper[i] = 2*(radius[i]+1);
710 padFilter->SetPadLowerBound(lower);
711 padFilter->SetPadUpperBound(upper);
713 output = padFilter->GetOutput();
717 typedef itk::BinaryDilateImageFilter<ImageType, ImageType , KernelType> DilateFilterType;
718 typename DilateFilterType::Pointer dilateFilter = DilateFilterType::New();
719 dilateFilter->SetBackgroundValue(BG);
720 dilateFilter->SetForegroundValue(FG);
721 dilateFilter->SetBoundaryToForeground(false);
722 dilateFilter->SetKernel(structuringElement);
723 dilateFilter->SetInput(output);
724 dilateFilter->Update();
725 return dilateFilter->GetOutput();
727 //--------------------------------------------------------------------
730 //--------------------------------------------------------------------
731 template<class ValueType, class VectorType>
732 void ConvertOption(std::string optionName, uint given,
733 ValueType * values, VectorType & p,
734 uint dim, bool required)
736 if (required && (given == 0)) {
737 clitkExceptionMacro("The option --" << optionName << " must be set and have 1 or "
738 << dim << " values.");
741 for(uint i=0; i<dim; i++) p[i] = values[0];
745 for(uint i=0; i<dim; i++) p[i] = values[i];
748 if (given == 0) return;
749 clitkExceptionMacro("The option --" << optionName << " must have 1 or "
750 << dim << " values.");
752 //--------------------------------------------------------------------
755 //--------------------------------------------------------------------
757 http://www.gamedev.net/community/forums/topic.asp?topic_id=542870
758 Assuming the points are (Ax,Ay) (Bx,By) and (Cx,Cy), you need to compute:
759 (Bx - Ax) * (Cy - Ay) - (By - Ay) * (Cx - Ax)
760 This will equal zero if the point C is on the line formed by
761 points A and B, and will have a different sign depending on the
762 side. Which side this is depends on the orientation of your (x,y)
763 coordinates, but you can plug test values for A,B and C into this
764 formula to determine whether negative values are to the left or to
766 => to accelerate, start with formula, when change sign -> stop and fill
768 offsetToKeep = is used to determine which side of the line we
769 keep. The point along the mainDirection but 'offsetToKeep' mm away
773 template<class ImageType>
775 SliceBySliceSetBackgroundFromLineSeparation(ImageType * input,
776 std::vector<typename ImageType::PointType> & lA,
777 std::vector<typename ImageType::PointType> & lB,
778 typename ImageType::PixelType BG,
782 typedef itk::ImageSliceIteratorWithIndex<ImageType> SliceIteratorType;
783 SliceIteratorType siter = SliceIteratorType(input,
784 input->GetLargestPossibleRegion());
785 siter.SetFirstDirection(0);
786 siter.SetSecondDirection(1);
789 typename ImageType::PointType A;
790 typename ImageType::PointType B;
791 typename ImageType::PointType C;
792 assert(lA.size() == lB.size());
793 while ((i<lA.size()) && (!siter.IsAtEnd())) {
794 // Check that the current slice correspond to the current point
795 input->TransformIndexToPhysicalPoint(siter.GetIndex(), C);
799 if ((fabs(C[2] - lA[i][2]))>0.01) { // is !equal with a tolerance of 0.01 mm
802 // Define A,B,C points
810 // Check that the line is not a point (A=B)
811 bool p = (A[0] == B[0]) && (A[1] == B[1]);
814 C[mainDirection] += offsetToKeep; // I know I must keep this point
815 double s = (B[0] - A[0]) * (C[1] - A[1]) - (B[1] - A[1]) * (C[0] - A[0]);
816 bool isPositive = s<0;
817 while (!siter.IsAtEndOfSlice()) {
818 while (!siter.IsAtEndOfLine()) {
819 // Very slow, I know ... but image should be very small
820 input->TransformIndexToPhysicalPoint(siter.GetIndex(), C);
821 double s = (B[0] - A[0]) * (C[1] - A[1]) - (B[1] - A[1]) * (C[0] - A[0]);
822 if (s == 0) siter.Set(BG); // on the line, we decide to remove
824 if (s > 0) siter.Set(BG);
827 if (s < 0) siter.Set(BG);
836 } // End of current slice
840 //--------------------------------------------------------------------
842 //--------------------------------------------------------------------
843 template<class ImageType>
845 AndNot(ImageType * input,
846 const ImageType * object,
847 typename ImageType::PixelType BG)
849 typename ImageType::Pointer o;
851 if (!clitk::HaveSameSizeAndSpacing<ImageType, ImageType>(input, object)) {
852 o = clitk::ResizeImageLike<ImageType>(object, input, BG);
856 typedef clitk::BooleanOperatorLabelImageFilter<ImageType> BoolFilterType;
857 typename BoolFilterType::Pointer boolFilter = BoolFilterType::New();
858 boolFilter->InPlaceOn();
859 boolFilter->SetInput1(input);
860 if (resized) boolFilter->SetInput2(o);
861 else boolFilter->SetInput2(object);
862 boolFilter->SetBackgroundValue1(BG);
863 boolFilter->SetBackgroundValue2(BG);
864 boolFilter->SetOperationType(BoolFilterType::AndNot);
865 boolFilter->Update();
867 //--------------------------------------------------------------------
870 //--------------------------------------------------------------------
871 template<class ImageType>
872 typename ImageType::Pointer
873 Binarize(const ImageType * input,
874 typename ImageType::PixelType lower,
875 typename ImageType::PixelType upper,
876 typename ImageType::PixelType BG,
877 typename ImageType::PixelType FG)
879 typedef itk::BinaryThresholdImageFilter<ImageType, ImageType> BinaryThresholdFilterType;
880 typename BinaryThresholdFilterType::Pointer binarizeFilter = BinaryThresholdFilterType::New();
881 binarizeFilter->SetInput(input);
882 binarizeFilter->SetLowerThreshold(lower);
883 binarizeFilter->SetUpperThreshold(upper);
884 binarizeFilter->SetInsideValue(FG);
885 binarizeFilter->SetOutsideValue(BG);
886 binarizeFilter->Update();
887 return binarizeFilter->GetOutput();
889 //--------------------------------------------------------------------
892 //--------------------------------------------------------------------
893 template<class ImageType>
895 GetMinMaxPointPosition(const ImageType * input,
896 typename ImageType::PointType & min,
897 typename ImageType::PointType & max)
899 typename ImageType::IndexType index = input->GetLargestPossibleRegion().GetIndex();
900 input->TransformIndexToPhysicalPoint(index, min);
901 index = index+input->GetLargestPossibleRegion().GetSize();
902 input->TransformIndexToPhysicalPoint(index, max);
904 //--------------------------------------------------------------------
907 //--------------------------------------------------------------------
908 template<class ImageType>
909 typename ImageType::PointType
910 FindExtremaPointInAGivenLine(const ImageType * input,
913 typename ImageType::PointType p,
914 typename ImageType::PixelType BG,
917 // Which direction ? Increasing or decreasing.
921 // Transform to pixel index
922 typename ImageType::IndexType index;
923 input->TransformPhysicalPointToIndex(p, index);
925 // Loop while inside the mask;
926 while (input->GetPixel(index) != BG) {
927 index[dimension] += d;
930 // Transform back to Physical Units
931 typename ImageType::PointType result;
932 input->TransformIndexToPhysicalPoint(index, result);
934 // Check that is is not too far away
935 double distance = p.EuclideanDistanceTo(result);
936 if (distance > distanceMax) {
937 result = p; // Get back to initial value
942 //--------------------------------------------------------------------
945 //--------------------------------------------------------------------
946 template<class PointType>
948 IsOnTheSameLineSide(PointType C, PointType A, PointType B, PointType like)
950 // Look at the position of point 'like' according to the AB line
951 double s = (B[0] - A[0]) * (like[1] - A[1]) - (B[1] - A[1]) * (like[0] - A[0]);
954 // Look the C position
955 s = (B[0] - A[0]) * (C[1] - A[1]) - (B[1] - A[1]) * (C[0] - A[0]);
957 if (negative && (s<=0)) return true;
958 if (!negative && (s>=0)) return true;
961 //--------------------------------------------------------------------
964 //--------------------------------------------------------------------
965 /* Consider an input object, for each slice, find the extrema
966 position according to a given direction and build a line segment
967 passing throught this point in a given direction. Output is a
968 vector of line (from point A to B), for each slice;
970 template<class ImageType>
972 SliceBySliceBuildLineSegmentAccordingToExtremaPosition(const ImageType * input,
973 typename ImageType::PixelType BG,
975 int extremaDirection,
976 bool extremaOppositeFlag,
979 std::vector<typename ImageType::PointType> & A,
980 std::vector<typename ImageType::PointType> & B)
983 typedef typename itk::Image<typename ImageType::PixelType, ImageType::ImageDimension-1> SliceType;
985 // Build the list of slices
986 std::vector<typename SliceType::Pointer> slices;
987 clitk::ExtractSlices<ImageType>(input, sliceDimension, slices);
989 // Build the list of 2D points
990 std::map<int, typename SliceType::PointType> position2D;
991 for(uint i=0; i<slices.size(); i++) {
992 typename SliceType::PointType p;
994 clitk::FindExtremaPointInAGivenDirection<SliceType>(slices[i], BG,
995 extremaDirection, extremaOppositeFlag, p);
1001 // Convert 2D points in slice into 3D points
1002 clitk::PointsUtils<ImageType>::Convert2DMapTo3DList(position2D, input, A);
1004 // Create additional point just right to the previous ones, on the
1005 // given lineDirection, in order to create a horizontal/vertical line.
1006 for(uint i=0; i<A.size(); i++) {
1007 typename ImageType::PointType p = A[i];
1008 p[lineDirection] += 10;
1016 //--------------------------------------------------------------------
1019 } // end of namespace