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://oncora1.lyon.fnclcc.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(typename ImageType::Pointer 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(typename ImageType::Pointer 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++) {
114 // DD((maxs[i] - mins[i])/image->GetSpacing()[i]);
115 regionSize[i] = lrint((maxs[i] - mins[i])/image->GetSpacing()[i]);
116 // DD(regionSize[i]);
120 region.SetIndex(regionStart);
121 region.SetSize(regionSize);
123 //--------------------------------------------------------------------
125 //--------------------------------------------------------------------
126 template<class ImageType, class TMaskImageType>
127 typename ImageType::Pointer
128 SetBackground(const ImageType * input,
129 const TMaskImageType * mask,
130 typename TMaskImageType::PixelType maskBG,
131 typename ImageType::PixelType outValue,
133 typedef SetBackgroundImageFilter<ImageType, TMaskImageType, ImageType>
134 SetBackgroundImageFilterType;
135 typename SetBackgroundImageFilterType::Pointer setBackgroundFilter
136 = SetBackgroundImageFilterType::New();
137 // if (inPlace) setBackgroundFilter->ReleaseDataFlagOn(); // No seg fault
138 setBackgroundFilter->SetInPlace(inPlace); // This is important to keep memory low
139 setBackgroundFilter->SetInput(input);
140 setBackgroundFilter->SetInput2(mask);
141 setBackgroundFilter->SetMaskValue(maskBG);
142 setBackgroundFilter->SetOutsideValue(outValue);
143 setBackgroundFilter->Update();
144 return setBackgroundFilter->GetOutput();
146 //--------------------------------------------------------------------
149 //--------------------------------------------------------------------
150 template<class ImageType>
151 int GetNumberOfConnectedComponentLabels(typename ImageType::Pointer input,
152 typename ImageType::PixelType BG,
153 bool isFullyConnected) {
154 // Connected Component label
155 typedef itk::ConnectedComponentImageFilter<ImageType, ImageType> ConnectFilterType;
156 typename ConnectFilterType::Pointer connectFilter = ConnectFilterType::New();
157 connectFilter->SetInput(input);
158 connectFilter->SetBackgroundValue(BG);
159 connectFilter->SetFullyConnected(isFullyConnected);
160 connectFilter->Update();
163 return connectFilter->GetObjectCount();
165 //--------------------------------------------------------------------
167 //--------------------------------------------------------------------
169 Warning : in this cas, we consider outputType like inputType, not
170 InternalImageType. Be sure it fits.
172 template<class ImageType>
173 typename ImageType::Pointer
174 Labelize(const ImageType * input,
175 typename ImageType::PixelType BG,
176 bool isFullyConnected,
177 int minimalComponentSize) {
178 // InternalImageType for storing large number of component
179 typedef itk::Image<int, ImageType::ImageDimension> InternalImageType;
181 // Connected Component label
182 typedef itk::ConnectedComponentImageFilter<ImageType, InternalImageType> ConnectFilterType;
183 typename ConnectFilterType::Pointer connectFilter = ConnectFilterType::New();
184 // connectFilter->ReleaseDataFlagOn();
185 connectFilter->SetInput(input);
186 connectFilter->SetBackgroundValue(BG);
187 connectFilter->SetFullyConnected(isFullyConnected);
189 // Sort by size and remove too small area.
190 typedef itk::RelabelComponentImageFilter<InternalImageType, ImageType> RelabelFilterType;
191 typename RelabelFilterType::Pointer relabelFilter = RelabelFilterType::New();
192 // relabelFilter->ReleaseDataFlagOn(); // if yes, fail when ExplosionControlledThresholdConnectedImageFilter ???
193 relabelFilter->SetInput(connectFilter->GetOutput());
194 relabelFilter->SetMinimumObjectSize(minimalComponentSize);
195 relabelFilter->Update();
197 // DD(relabelFilter->GetNumberOfObjects());
198 // DD(relabelFilter->GetOriginalNumberOfObjects());
199 // DD(relabelFilter->GetSizeOfObjectsInPhysicalUnits()[0]);
202 typename ImageType::Pointer output = relabelFilter->GetOutput();
205 //--------------------------------------------------------------------
208 //--------------------------------------------------------------------
210 Warning : in this cas, we consider outputType like inputType, not
211 InternalImageType. Be sure it fits.
213 template<class ImageType>
214 typename ImageType::Pointer
215 LabelizeAndCountNumberOfObjects(const ImageType * input,
216 typename ImageType::PixelType BG,
217 bool isFullyConnected,
218 int minimalComponentSize,
220 // InternalImageType for storing large number of component
221 typedef itk::Image<int, ImageType::ImageDimension> InternalImageType;
223 // Connected Component label
224 typedef itk::ConnectedComponentImageFilter<ImageType, InternalImageType> ConnectFilterType;
225 typename ConnectFilterType::Pointer connectFilter = ConnectFilterType::New();
226 // connectFilter->ReleaseDataFlagOn();
227 connectFilter->SetInput(input);
228 connectFilter->SetBackgroundValue(BG);
229 connectFilter->SetFullyConnected(isFullyConnected);
231 // Sort by size and remove too small area.
232 typedef itk::RelabelComponentImageFilter<InternalImageType, ImageType> RelabelFilterType;
233 typename RelabelFilterType::Pointer relabelFilter = RelabelFilterType::New();
234 // relabelFilter->ReleaseDataFlagOn(); // if yes, fail when ExplosionControlledThresholdConnectedImageFilter ???
235 relabelFilter->SetInput(connectFilter->GetOutput());
236 relabelFilter->SetMinimumObjectSize(minimalComponentSize);
237 relabelFilter->Update();
239 nb = relabelFilter->GetNumberOfObjects();
240 // DD(relabelFilter->GetNumberOfObjects());
241 // DD(relabelFilter->GetOriginalNumberOfObjects());
242 // DD(relabelFilter->GetSizeOfObjectsInPhysicalUnits()[0]);
245 typename ImageType::Pointer output = relabelFilter->GetOutput();
248 //--------------------------------------------------------------------
252 //--------------------------------------------------------------------
253 template<class ImageType>
254 typename ImageType::Pointer
255 RemoveLabels(typename ImageType::Pointer input,
256 typename ImageType::PixelType BG,
257 std::vector<typename ImageType::PixelType> & labelsToRemove) {
258 typename ImageType::Pointer working_image = input;
259 for (unsigned int i=0; i <labelsToRemove.size(); i++) {
260 typedef SetBackgroundImageFilter<ImageType, ImageType> SetBackgroundImageFilterType;
261 typename SetBackgroundImageFilterType::Pointer setBackgroundFilter = SetBackgroundImageFilterType::New();
262 setBackgroundFilter->SetInput(input);
263 setBackgroundFilter->SetInput2(input);
264 setBackgroundFilter->SetMaskValue(labelsToRemove[i]);
265 setBackgroundFilter->SetOutsideValue(BG);
266 setBackgroundFilter->Update();
267 working_image = setBackgroundFilter->GetOutput();
269 return working_image;
271 //--------------------------------------------------------------------
274 //--------------------------------------------------------------------
275 template<class ImageType>
276 typename ImageType::Pointer
277 KeepLabels(const ImageType * input,
278 typename ImageType::PixelType BG,
279 typename ImageType::PixelType FG,
280 typename ImageType::PixelType firstKeep,
281 typename ImageType::PixelType lastKeep,
283 typedef itk::BinaryThresholdImageFilter<ImageType, ImageType> BinarizeFilterType;
284 typename BinarizeFilterType::Pointer binarizeFilter = BinarizeFilterType::New();
285 binarizeFilter->SetInput(input);
286 binarizeFilter->SetLowerThreshold(firstKeep);
287 if (useLastKeep) binarizeFilter->SetUpperThreshold(lastKeep);
288 binarizeFilter->SetInsideValue(FG);
289 binarizeFilter->SetOutsideValue(BG);
290 binarizeFilter->Update();
291 return binarizeFilter->GetOutput();
293 //--------------------------------------------------------------------
296 //--------------------------------------------------------------------
297 template<class ImageType>
298 typename ImageType::Pointer
299 LabelizeAndSelectLabels(typename ImageType::Pointer input,
300 typename ImageType::PixelType BG,
301 typename ImageType::PixelType FG,
302 bool isFullyConnected,
303 int minimalComponentSize,
304 LabelizeParameters<typename ImageType::PixelType> * param)
306 typename ImageType::Pointer working_image;
307 working_image = Labelize<ImageType>(input, BG, isFullyConnected, minimalComponentSize);
308 working_image = RemoveLabels<ImageType>(working_image, BG, param->GetLabelsToRemove());
309 working_image = KeepLabels<ImageType>(working_image,
311 param->GetFirstKeep(),
312 param->GetLastKeep(),
313 param->GetUseLastKeep());
314 return working_image;
316 //--------------------------------------------------------------------
319 //--------------------------------------------------------------------
320 template<class ImageType>
321 typename ImageType::Pointer
322 ResizeImageLike(typename ImageType::Pointer input,
323 typename ImageType::Pointer like,
324 typename ImageType::PixelType backgroundValue)
326 typedef CropLikeImageFilter<ImageType> CropFilterType;
327 typename CropFilterType::Pointer cropFilter = CropFilterType::New();
328 cropFilter->SetInput(input);
329 cropFilter->SetCropLikeImage(like);
330 cropFilter->SetBackgroundValue(backgroundValue);
331 cropFilter->Update();
332 return cropFilter->GetOutput();
334 //--------------------------------------------------------------------
337 //--------------------------------------------------------------------
338 template<class MaskImageType>
339 typename MaskImageType::Pointer
340 SliceBySliceRelativePosition(const MaskImageType * input,
341 const MaskImageType * object,
344 std::string orientation,
345 bool uniqueConnectedComponent,
349 typedef SliceBySliceRelativePositionFilter<MaskImageType> SliceRelPosFilterType;
350 typename SliceRelPosFilterType::Pointer sliceRelPosFilter = SliceRelPosFilterType::New();
351 sliceRelPosFilter->VerboseStepFlagOff();
352 sliceRelPosFilter->WriteStepFlagOff();
353 sliceRelPosFilter->SetInput(input);
354 sliceRelPosFilter->SetInputObject(object);
355 sliceRelPosFilter->SetDirection(direction);
356 sliceRelPosFilter->SetFuzzyThreshold(threshold);
357 sliceRelPosFilter->AddOrientationTypeString(orientation);
358 sliceRelPosFilter->SetResampleBeforeRelativePositionFilter((spacing != -1));
359 sliceRelPosFilter->SetIntermediateSpacing(spacing);
360 sliceRelPosFilter->SetUniqueConnectedComponentBySlice(uniqueConnectedComponent);
361 sliceRelPosFilter->SetInverseOrientationFlag(inverseflag);
362 // sliceRelPosFilter->SetAutoCropFlag(true); ??
363 sliceRelPosFilter->Update();
364 return sliceRelPosFilter->GetOutput();
366 //--------------------------------------------------------------------
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 //--------------------------------------------------------------------
382 //--------------------------------------------------------------------
383 template<class ImageType>
385 FindExtremaPointInAGivenDirection(const ImageType * input,
386 typename ImageType::PixelType bg,
387 int direction, bool opposite,
388 typename ImageType::PointType refpoint,
390 typename ImageType::PointType & point)
393 loop over input pixels, store the index in the fg that is max
394 according to the given direction.
396 typedef itk::ImageRegionConstIteratorWithIndex<ImageType> IteratorType;
397 IteratorType iter(input, input->GetLargestPossibleRegion());
399 typename ImageType::IndexType max = input->GetLargestPossibleRegion().GetIndex();
400 if (opposite) max = max+input->GetLargestPossibleRegion().GetSize();
402 while (!iter.IsAtEnd()) {
403 if (iter.Get() != bg) {
404 bool test = iter.GetIndex()[direction] > max[direction];
405 if (opposite) test = !test;
407 typename ImageType::PointType p;
408 input->TransformIndexToPhysicalPoint(iter.GetIndex(), p);
409 if ((distanceMax==0) || (p.EuclideanDistanceTo(refpoint) < distanceMax)) {
410 max = iter.GetIndex();
417 if (!found) return false;
418 input->TransformIndexToPhysicalPoint(max, point);
421 //--------------------------------------------------------------------
424 //--------------------------------------------------------------------
425 template<class ImageType>
426 typename ImageType::Pointer
427 CropImageAbove(typename ImageType::Pointer image,
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 CropImageBelow(typename ImageType::Pointer image,
446 typename ImageType::PixelType BG)
448 typename ImageType::PointType p;
449 image->TransformIndexToPhysicalPoint(image->GetLargestPossibleRegion().GetIndex()+
450 image->GetLargestPossibleRegion().GetSize(), p);
451 return CropImageAlongOneAxis<ImageType>(image, dim, max, p[dim], autoCrop, BG);
453 //--------------------------------------------------------------------
456 //--------------------------------------------------------------------
457 template<class ImageType>
458 typename ImageType::Pointer
459 CropImageAlongOneAxis(typename ImageType::Pointer image,
460 int dim, double min, double max,
462 typename ImageType::PixelType BG)
464 // Compute region size
465 typename ImageType::RegionType region;
466 typename ImageType::SizeType size = image->GetLargestPossibleRegion().GetSize();
467 typename ImageType::PointType p = image->GetOrigin();
469 typename ImageType::IndexType start;
470 image->TransformPhysicalPointToIndex(p, start);
472 typename ImageType::IndexType end;
473 image->TransformPhysicalPointToIndex(p, end);
474 size[dim] = fabs(end[dim]-start[dim]);
475 region.SetIndex(start);
476 region.SetSize(size);
479 typedef itk::RegionOfInterestImageFilter<ImageType, ImageType> CropFilterType;
480 typename CropFilterType::Pointer cropFilter = CropFilterType::New();
481 cropFilter->SetInput(image);
482 cropFilter->SetRegionOfInterest(region);
483 cropFilter->Update();
484 typename ImageType::Pointer result = cropFilter->GetOutput();
488 result = AutoCrop<ImageType>(result, BG);
492 //--------------------------------------------------------------------
495 //--------------------------------------------------------------------
496 template<class ImageType>
498 ComputeCentroids(typename ImageType::Pointer image,
499 typename ImageType::PixelType BG,
500 std::vector<typename ImageType::PointType> & centroids)
502 typedef long LabelType;
503 static const unsigned int Dim = ImageType::ImageDimension;
504 typedef itk::ShapeLabelObject< LabelType, Dim > LabelObjectType;
505 typedef itk::LabelMap< LabelObjectType > LabelMapType;
506 typedef itk::LabelImageToLabelMapFilter<ImageType, LabelMapType> ImageToMapFilterType;
507 typename ImageToMapFilterType::Pointer imageToLabelFilter = ImageToMapFilterType::New();
508 typedef itk::ShapeLabelMapFilter<LabelMapType, ImageType> ShapeFilterType;
509 typename ShapeFilterType::Pointer statFilter = ShapeFilterType::New();
510 imageToLabelFilter->SetBackgroundValue(BG);
511 imageToLabelFilter->SetInput(image);
512 statFilter->SetInput(imageToLabelFilter->GetOutput());
513 statFilter->Update();
514 typename LabelMapType::Pointer labelMap = statFilter->GetOutput();
517 typename ImageType::PointType dummy;
518 centroids.push_back(dummy); // label 0 -> no centroid, use dummy point
519 for(uint i=1; i<labelMap->GetNumberOfLabelObjects()+1; i++) {
520 centroids.push_back(labelMap->GetLabelObject(i)->GetCentroid());
523 //--------------------------------------------------------------------
526 //--------------------------------------------------------------------
527 template<class ImageType>
529 ExtractSlices(typename ImageType::Pointer image,
531 std::vector<typename itk::Image<typename ImageType::PixelType,
532 ImageType::ImageDimension-1>::Pointer > & slices)
534 typedef ExtractSliceFilter<ImageType> ExtractSliceFilterType;
535 typedef typename ExtractSliceFilterType::SliceType SliceType;
536 typename ExtractSliceFilterType::Pointer
537 extractSliceFilter = ExtractSliceFilterType::New();
538 extractSliceFilter->SetInput(image);
539 extractSliceFilter->SetDirection(direction);
540 extractSliceFilter->Update();
541 extractSliceFilter->GetOutputSlices(slices);
543 //--------------------------------------------------------------------
546 //--------------------------------------------------------------------
547 template<class ImageType>
548 typename ImageType::Pointer
549 JoinSlices(std::vector<typename itk::Image<typename ImageType::PixelType,
550 ImageType::ImageDimension-1>::Pointer > & slices,
551 typename ImageType::Pointer input,
553 typedef typename itk::Image<typename ImageType::PixelType, ImageType::ImageDimension-1> SliceType;
554 typedef itk::JoinSeriesImageFilter<SliceType, ImageType> JoinSeriesFilterType;
555 typename JoinSeriesFilterType::Pointer joinFilter = JoinSeriesFilterType::New();
556 joinFilter->SetOrigin(input->GetOrigin()[direction]);
557 joinFilter->SetSpacing(input->GetSpacing()[direction]);
558 for(unsigned int i=0; i<slices.size(); i++) {
559 joinFilter->PushBackInput(slices[i]);
561 joinFilter->Update();
562 return joinFilter->GetOutput();
564 //--------------------------------------------------------------------
567 //--------------------------------------------------------------------
568 template<class ImageType>
570 PointsUtils<ImageType>::Convert2DTo3D(const PointType2D & p,
577 p3D[2] = (image->GetLargestPossibleRegion().GetIndex()[2]+slice)*image->GetSpacing()[2]
578 + image->GetOrigin()[2];
580 //--------------------------------------------------------------------
583 //--------------------------------------------------------------------
584 template<class ImageType>
586 PointsUtils<ImageType>::Convert2DTo3DList(const MapPoint2DType & map,
588 VectorPoint3DType & list)
590 typename MapPoint2DType::const_iterator iter = map.begin();
591 while (iter != map.end()) {
593 Convert2DTo3D(iter->second, image, iter->first, p);
598 //--------------------------------------------------------------------
600 //--------------------------------------------------------------------
601 template<class ImageType>
603 WriteListOfLandmarks(std::vector<typename ImageType::PointType> points,
604 std::string filename)
607 openFileForWriting(os, filename);
608 os << "LANDMARKS1" << std::endl;
609 for(uint i=0; i<points.size(); i++) {
610 const typename ImageType::PointType & p = points[i];
611 // Write it in the file
612 os << i << " " << p[0] << " " << p[1] << " " << p[2] << " 0 0 " << std::endl;
616 //--------------------------------------------------------------------
619 //--------------------------------------------------------------------
620 template<class ImageType>
621 typename ImageType::Pointer
622 Dilate(typename ImageType::Pointer image,
624 typename ImageType::PixelType BG,
625 typename ImageType::PixelType FG,
628 typename ImageType::SizeType r;
629 for(uint i=0; i<ImageType::ImageDimension; i++)
630 r[i] = (uint)lrint(radiusInMM/image->GetSpacing()[i]);
631 return Dilate<ImageType>(image, r, BG, FG, extendSupport);
633 //--------------------------------------------------------------------
636 //--------------------------------------------------------------------
637 template<class ImageType>
638 typename ImageType::Pointer
639 Dilate(typename ImageType::Pointer image,
640 typename ImageType::PointType radiusInMM,
641 typename ImageType::PixelType BG,
642 typename ImageType::PixelType FG,
645 typename ImageType::SizeType r;
646 for(uint i=0; i<ImageType::ImageDimension; i++)
647 r[i] = (uint)lrint(radiusInMM[i]/image->GetSpacing()[i]);
648 return Dilate<ImageType>(image, r, BG, FG, extendSupport);
650 //--------------------------------------------------------------------
653 //--------------------------------------------------------------------
654 template<class ImageType>
655 typename ImageType::Pointer
656 Dilate(typename ImageType::Pointer image,
657 typename ImageType::SizeType radius,
658 typename ImageType::PixelType BG,
659 typename ImageType::PixelType FG,
662 // Create kernel for dilatation
663 typedef itk::BinaryBallStructuringElement<typename ImageType::PixelType,
664 ImageType::ImageDimension> KernelType;
665 KernelType structuringElement;
666 structuringElement.SetRadius(radius);
667 structuringElement.CreateStructuringElement();
670 typedef itk::ConstantPadImageFilter<ImageType, ImageType> PadFilterType;
671 typename PadFilterType::Pointer padFilter = PadFilterType::New();
672 padFilter->SetInput(image);
673 typename ImageType::SizeType lower;
674 typename ImageType::SizeType upper;
675 for(uint i=0; i<3; i++) {
676 lower[i] = upper[i] = 2*(radius[i]+1);
678 padFilter->SetPadLowerBound(lower);
679 padFilter->SetPadUpperBound(upper);
681 image = padFilter->GetOutput();
685 typedef itk::BinaryDilateImageFilter<ImageType, ImageType , KernelType> DilateFilterType;
686 typename DilateFilterType::Pointer dilateFilter = DilateFilterType::New();
687 dilateFilter->SetBackgroundValue(BG);
688 dilateFilter->SetForegroundValue(FG);
689 dilateFilter->SetBoundaryToForeground(false);
690 dilateFilter->SetKernel(structuringElement);
691 dilateFilter->SetInput(image);
692 dilateFilter->Update();
693 return image = dilateFilter->GetOutput();
695 //--------------------------------------------------------------------
698 //--------------------------------------------------------------------
699 template<class ValueType, class VectorType>
700 void ConvertOption(std::string optionName, uint given,
701 ValueType * values, VectorType & p,
702 uint dim, bool required)
704 if (required && (given == 0)) {
705 clitkExceptionMacro("The option --" << optionName << " must be set and have 1 or "
706 << dim << " values.");
709 for(uint i=0; i<dim; i++) p[i] = values[0];
713 for(uint i=0; i<dim; i++) p[i] = values[i];
716 if (given == 0) return;
717 clitkExceptionMacro("The option --" << optionName << " must have 1 or "
718 << dim << " values.");
720 //--------------------------------------------------------------------
723 //--------------------------------------------------------------------
725 http://www.gamedev.net/community/forums/topic.asp?topic_id=542870
726 Assuming the points are (Ax,Ay) (Bx,By) and (Cx,Cy), you need to compute:
727 (Bx - Ax) * (Cy - Ay) - (By - Ay) * (Cx - Ax)
728 This will equal zero if the point C is on the line formed by
729 points A and B, and will have a different sign depending on the
730 side. Which side this is depends on the orientation of your (x,y)
731 coordinates, but you can plug test values for A,B and C into this
732 formula to determine whether negative values are to the left or to
734 => to accelerate, start with formula, when change sign -> stop and fill
736 template<class ImageType>
738 SliceBySliceSetBackgroundFromLineSeparation(typename ImageType::Pointer input,
739 std::vector<typename ImageType::PointType> & lA,
740 std::vector<typename ImageType::PointType> & lB,
741 typename ImageType::PixelType BG,
746 typedef itk::ImageSliceIteratorWithIndex<ImageType> SliceIteratorType;
747 SliceIteratorType siter = SliceIteratorType(input,
748 input->GetLargestPossibleRegion());
749 siter.SetFirstDirection(0);
750 siter.SetSecondDirection(1);
753 typename ImageType::PointType A;
754 typename ImageType::PointType B;
755 typename ImageType::PointType C;
756 while (!siter.IsAtEnd()) {
757 // Check that the current slice correspond to the current point
758 input->TransformIndexToPhysicalPoint(siter.GetIndex(), C);
759 if (C[2] != lA[i][2]) {
764 // Define A,B,C points
768 C[mainDirection] += offsetToKeep; // I know I must keep this point
769 double s = (B[0] - A[0]) * (C[1] - A[1]) - (B[1] - A[1]) * (C[0] - A[0]);
770 bool isPositive = s<0;
771 while (!siter.IsAtEndOfSlice()) {
772 while (!siter.IsAtEndOfLine()) {
773 // Very slow, I know ... but image should be very small
774 input->TransformIndexToPhysicalPoint(siter.GetIndex(), C);
775 double s = (B[0] - A[0]) * (C[1] - A[1]) - (B[1] - A[1]) * (C[0] - A[0]);
776 if (s == 0) siter.Set(BG); // on the line, we decide to remove
778 if (s > 0) siter.Set(BG);
781 if (s < 0) siter.Set(BG);
792 //--------------------------------------------------------------------