From 9a4f7a1a65fcaa160878e50a12a61bdacf78607a Mon Sep 17 00:00:00 2001 From: dsarrut Date: Thu, 3 Mar 2011 10:11:38 +0000 Subject: [PATCH] replace SmartPointer with raw pointer for function (recommended) --- itk/clitkResampleImageWithOptionsFilter.txx | 15 +- itk/clitkSegmentationUtils.h | 136 +- itk/clitkSegmentationUtils.txx | 1594 ++++++++++--------- 3 files changed, 957 insertions(+), 788 deletions(-) diff --git a/itk/clitkResampleImageWithOptionsFilter.txx b/itk/clitkResampleImageWithOptionsFilter.txx index 7c05dc5..0807fc0 100644 --- a/itk/clitkResampleImageWithOptionsFilter.txx +++ b/itk/clitkResampleImageWithOptionsFilter.txx @@ -55,6 +55,7 @@ ResampleImageWithOptionsFilter():itk::ImageToImageFilterCopyInformation(input); outputImage->SetLargestPossibleRegion(m_OutputRegion); outputImage->SetSpacing(m_OutputSpacing); + outputImage->FillBuffer(GetDefaultPixelValue()); // Init Gaussian sigma if (m_GaussianSigma[0] != -1) { // Gaussian filter set by user @@ -188,19 +191,11 @@ GenerateData() InputImagePointer input = dynamic_cast(itk::ProcessObject::GetInput(0)); static const unsigned int dim = InputImageType::ImageDimension; - // Set regions and allocate - //DD(this->GetOutput()->GetLargestPossibleRegion()); - //this->GetOutput()->SetRegions(m_OutputRegion); - //this->GetOutput()->Allocate(); - // this->GetOutput()->FillBuffer(m_DefaultPixelValue); - // Create main Resample Image Filter typedef itk::ResampleImageFilter FilterType; typename FilterType::Pointer filter = FilterType::New(); filter->GraftOutput(this->GetOutput()); - // this->GetOutput()->Print(std::cout); this->GetOutput()->SetBufferedRegion(this->GetOutput()->GetLargestPossibleRegion()); - // this->GetOutput()->Print(std::cout); // Print options if needed if (m_VerboseOptions) { diff --git a/itk/clitkSegmentationUtils.h b/itk/clitkSegmentationUtils.h index aafe722..6fa67aa 100644 --- a/itk/clitkSegmentationUtils.h +++ b/itk/clitkSegmentationUtils.h @@ -27,11 +27,17 @@ // itk #include +/* + According to + http://answerpot.com/showthread.php?357451-Itk::SmartPointer%20-%20problem%20making%20code%20const-correct + it is better to take raw pointer as argument instead of SmartPointer. +*/ + namespace clitk { //-------------------------------------------------------------------- template - void ComputeBBFromImageRegion(typename ImageType::Pointer image, + void ComputeBBFromImageRegion(const ImageType * image, typename ImageType::RegionType region, typename itk::BoundingBox::Pointer bb); @@ -44,9 +50,9 @@ namespace clitk { //-------------------------------------------------------------------- template - void ComputeRegionFromBB(typename ImageType::Pointer image, + void ComputeRegionFromBB(const ImageType * image, const typename itk::BoundingBox::Pointer bb, + ImageType::ImageDimension>::Pointer bb, typename ImageType::RegionType & region); //-------------------------------------------------------------------- template @@ -61,7 +67,7 @@ namespace clitk { //-------------------------------------------------------------------- template - int GetNumberOfConnectedComponentLabels(typename ImageType::Pointer input, + int GetNumberOfConnectedComponentLabels(const ImageType * input, typename ImageType::PixelType BG, bool isFullyConnected); //-------------------------------------------------------------------- @@ -70,10 +76,8 @@ namespace clitk { //-------------------------------------------------------------------- template typename TImageType::Pointer - Labelize(const TImageType * input, - typename TImageType::PixelType BG, - bool isFullyConnected, - int minimalComponentSize); + Labelize(const TImageType * input, typename TImageType::PixelType BG, + bool isFullyConnected, int minimalComponentSize); template typename TImageType::Pointer LabelizeAndCountNumberOfObjects(const TImageType * input, @@ -87,7 +91,7 @@ namespace clitk { //-------------------------------------------------------------------- template typename ImageType::Pointer - RemoveLabels(typename ImageType::Pointer input, + RemoveLabels(const ImageType * input, typename ImageType::PixelType BG, std::vector & labelsToRemove); //-------------------------------------------------------------------- @@ -96,7 +100,7 @@ namespace clitk { //-------------------------------------------------------------------- template typename ImageType::Pointer - AutoCrop(typename ImageType::Pointer input, + AutoCrop(const ImageType * input, typename ImageType::PixelType BG) { typedef clitk::AutoCropFilter AutoCropFilterType; typename AutoCropFilterType::Pointer autoCropFilter = AutoCropFilterType::New(); @@ -123,7 +127,7 @@ namespace clitk { //-------------------------------------------------------------------- template typename TImageType::Pointer - LabelizeAndSelectLabels(typename TImageType::Pointer input, + LabelizeAndSelectLabels(const TImageType * input, typename TImageType::PixelType BG, typename TImageType::PixelType FG, bool isFullyConnected, @@ -133,8 +137,8 @@ namespace clitk { //-------------------------------------------------------------------- template typename ImageType::Pointer - ResizeImageLike(typename ImageType::Pointer input, - typename ImageType::Pointer like, + ResizeImageLike(const ImageType * input, + const itk::ImageBase * like, typename ImageType::PixelType BG); @@ -180,21 +184,19 @@ namespace clitk { //-------------------------------------------------------------------- template typename ImageType::Pointer - CropImageAlongOneAxis(typename ImageType::Pointer image, + CropImageAlongOneAxis(const ImageType * image, int dim, double min, double max, bool autoCrop = false, typename ImageType::PixelType BG=0); template typename ImageType::Pointer - CropImageAbove(typename ImageType::Pointer image, - int dim, double min, - bool autoCrop = false, + CropImageAbove(const ImageType * image, + int dim, double min, bool autoCrop = false, typename ImageType::PixelType BG=0); template typename ImageType::Pointer - CropImageBelow(typename ImageType::Pointer image, - int dim, double max, - bool autoCrop = false, + CropImageBelow(const ImageType * image, + int dim, double max,bool autoCrop = false, typename ImageType::PixelType BG=0); //-------------------------------------------------------------------- @@ -202,7 +204,7 @@ namespace clitk { //-------------------------------------------------------------------- template void - ComputeCentroids(typename ImageType::Pointer image, + ComputeCentroids(const ImageType * image, typename ImageType::PixelType BG, std::vector & centroids); //-------------------------------------------------------------------- @@ -211,10 +213,9 @@ namespace clitk { //-------------------------------------------------------------------- template void - ExtractSlices(typename ImageType::Pointer image, - int dim, + ExtractSlices(const ImageType * image, int dim, std::vector< typename itk::Image::Pointer > & slices); + ImageType::ImageDimension-1>::Pointer > & slices); //-------------------------------------------------------------------- @@ -222,9 +223,8 @@ namespace clitk { template typename ImageType::Pointer JoinSlices(std::vector::Pointer > & slices, - typename ImageType::Pointer input, - int dim); + ImageType::ImageDimension-1>::Pointer > & slices, + const ImageType * input, int dim); //-------------------------------------------------------------------- @@ -234,21 +234,23 @@ namespace clitk { class PointsUtils { typedef typename ImageType::PointType PointType3D; + typedef typename ImageType::IndexType IndexType3D; typedef typename ImageType::PixelType PixelType; typedef typename ImageType::Pointer ImagePointer; typedef typename ImageType::ConstPointer ImageConstPointer; typedef itk::Image SliceType; typedef typename SliceType::PointType PointType2D; + typedef typename SliceType::IndexType IndexType2D; typedef std::map MapPoint2DType; typedef std::vector VectorPoint3DType; public: static void Convert2DTo3D(const PointType2D & p2D, - ImagePointer image, + const ImageType * image, const int slice, PointType3D & p3D); static void Convert2DTo3DList(const MapPoint2DType & map, - ImagePointer image, + const ImageType * image, VectorPoint3DType & list); }; @@ -263,25 +265,22 @@ namespace clitk { //-------------------------------------------------------------------- template typename ImageType::Pointer - Dilate(typename ImageType::Pointer image, - double radiusInMM, - typename ImageType::PixelType BG, - typename ImageType::PixelType FG, - bool extendSupport); + Dilate(const ImageType * image, double radiusInMM, + typename ImageType::PixelType BG, + typename ImageType::PixelType FG, + bool extendSupport); template typename ImageType::Pointer - Dilate(typename ImageType::Pointer image, - typename ImageType::SizeType radius, - typename ImageType::PixelType BG, - typename ImageType::PixelType FG, - bool extendSupport); + Dilate(const ImageType * image, typename ImageType::SizeType radius, + typename ImageType::PixelType BG, + typename ImageType::PixelType FG, + bool extendSupport); template typename ImageType::Pointer - Dilate(typename ImageType::Pointer image, - typename ImageType::PointType radiusInMM, - typename ImageType::PixelType BG, - typename ImageType::PixelType FG, - bool extendSupport); + Dilate(const ImageType * image, typename ImageType::PointType radiusInMM, + typename ImageType::PixelType BG, + typename ImageType::PixelType FG, + bool extendSupport); //-------------------------------------------------------------------- //-------------------------------------------------------------------- @@ -289,20 +288,65 @@ namespace clitk { void ConvertOption(std::string optionName, uint given, ValueType * values, VectorType & p, uint dim, bool required); -#define ConvertOptionMacro(OPTIONNAME, VAR, DIM, REQUIRED) \ +#define ConvertOptionMacro(OPTIONNAME, VAR, DIM, REQUIRED) \ ConvertOption(#OPTIONNAME, OPTIONNAME##_given, OPTIONNAME##_arg, VAR, DIM, REQUIRED); //-------------------------------------------------------------------- //-------------------------------------------------------------------- template void - SliceBySliceSetBackgroundFromLineSeparation(typename ImageType::Pointer input, + SliceBySliceSetBackgroundFromLineSeparation(ImageType * input, std::vector & lA, std::vector & lB, typename ImageType::PixelType BG, int mainDirection, double offsetToKeep); + //-------------------------------------------------------------------- + + + //-------------------------------------------------------------------- + template + void AndNot(ImageType * input, + const ImageType * object, + typename ImageType::PixelType BG=0); + //-------------------------------------------------------------------- + + + //-------------------------------------------------------------------- + template + typename ImageType::Pointer + Binarize(const ImageType * input, + typename ImageType::PixelType lower, + typename ImageType::PixelType upper, + typename ImageType::PixelType BG=0, + typename ImageType::PixelType FG=1); + //-------------------------------------------------------------------- + + //-------------------------------------------------------------------- + template + void + GetMinMaxPointPosition(const ImageType * input, + typename ImageType::PointType & min, + typename ImageType::PointType & max); + //-------------------------------------------------------------------- + + //-------------------------------------------------------------------- + template + typename ImageType::PointType + FindExtremaPointInAGivenLine(const ImageType * input, + int dimension, bool inverse, + typename ImageType::PointType p, + typename ImageType::PixelType BG, + double distanceMax); + //-------------------------------------------------------------------- + + + //-------------------------------------------------------------------- + template + bool + IsOnTheSameLineSide(PointType C, PointType A, PointType B, PointType like); + //-------------------------------------------------------------------- } diff --git a/itk/clitkSegmentationUtils.txx b/itk/clitkSegmentationUtils.txx index 1760395..9b226b2 100644 --- a/itk/clitkSegmentationUtils.txx +++ b/itk/clitkSegmentationUtils.txx @@ -35,759 +35,889 @@ namespace clitk { -//-------------------------------------------------------------------- -template -void ComputeBBFromImageRegion(typename ImageType::Pointer image, - typename ImageType::RegionType region, - typename itk::BoundingBox::Pointer bb) { - typedef typename ImageType::IndexType IndexType; - IndexType firstIndex; - IndexType lastIndex; - for(unsigned int i=0; iGetImageDimension(); i++) { - firstIndex[i] = region.GetIndex()[i]; - lastIndex[i] = firstIndex[i]+region.GetSize()[i]; - } - - typedef itk::BoundingBox BBType; - typedef typename BBType::PointType PointType; - PointType lastPoint; - PointType firstPoint; - image->TransformIndexToPhysicalPoint(firstIndex, firstPoint); - image->TransformIndexToPhysicalPoint(lastIndex, lastPoint); - - bb->SetMaximum(lastPoint); - bb->SetMinimum(firstPoint); -} -//-------------------------------------------------------------------- - - -//-------------------------------------------------------------------- -template -void ComputeBBIntersection(typename itk::BoundingBox::Pointer bbo, - typename itk::BoundingBox::Pointer bbi1, - typename itk::BoundingBox::Pointer bbi2) { - - typedef itk::BoundingBox BBType; - typedef typename BBType::PointType PointType; - PointType lastPoint; - PointType firstPoint; - - for(unsigned int i=0; iGetMinimum()[i], - bbi2->GetMinimum()[i]); - lastPoint[i] = std::min(bbi1->GetMaximum()[i], - bbi2->GetMaximum()[i]); - } - - bbo->SetMaximum(lastPoint); - bbo->SetMinimum(firstPoint); -} -//-------------------------------------------------------------------- - - -//-------------------------------------------------------------------- -template -void ComputeRegionFromBB(typename ImageType::Pointer image, - const typename itk::BoundingBox::Pointer bb, - typename ImageType::RegionType & region) { - // Types - typedef typename ImageType::IndexType IndexType; - typedef typename ImageType::PointType PointType; - typedef typename ImageType::RegionType RegionType; - typedef typename ImageType::SizeType SizeType; - - // Region starting point - IndexType regionStart; - PointType start = bb->GetMinimum(); - image->TransformPhysicalPointToIndex(start, regionStart); - - // Region size - SizeType regionSize; - PointType maxs = bb->GetMaximum(); - PointType mins = bb->GetMinimum(); - for(unsigned int i=0; iGetSpacing()[i]); - regionSize[i] = lrint((maxs[i] - mins[i])/image->GetSpacing()[i]); - // DD(regionSize[i]); + //-------------------------------------------------------------------- + template + void ComputeBBFromImageRegion(const ImageType * image, + typename ImageType::RegionType region, + typename itk::BoundingBox::Pointer bb) { + typedef typename ImageType::IndexType IndexType; + IndexType firstIndex; + IndexType lastIndex; + for(unsigned int i=0; iGetImageDimension(); i++) { + firstIndex[i] = region.GetIndex()[i]; + lastIndex[i] = firstIndex[i]+region.GetSize()[i]; + } + + typedef itk::BoundingBox BBType; + typedef typename BBType::PointType PointType; + PointType lastPoint; + PointType firstPoint; + image->TransformIndexToPhysicalPoint(firstIndex, firstPoint); + image->TransformIndexToPhysicalPoint(lastIndex, lastPoint); + + bb->SetMaximum(lastPoint); + bb->SetMinimum(firstPoint); } + //-------------------------------------------------------------------- + + + //-------------------------------------------------------------------- + template + void ComputeBBIntersection(typename itk::BoundingBox::Pointer bbo, + typename itk::BoundingBox::Pointer bbi1, + typename itk::BoundingBox::Pointer bbi2) { + + typedef itk::BoundingBox BBType; + typedef typename BBType::PointType PointType; + PointType lastPoint; + PointType firstPoint; + + for(unsigned int i=0; iGetMinimum()[i], + bbi2->GetMinimum()[i]); + lastPoint[i] = std::min(bbi1->GetMaximum()[i], + bbi2->GetMaximum()[i]); + } + + bbo->SetMaximum(lastPoint); + bbo->SetMinimum(firstPoint); + } + //-------------------------------------------------------------------- + + + //-------------------------------------------------------------------- + template + void ComputeRegionFromBB(const ImageType * image, + const typename itk::BoundingBox::Pointer bb, + typename ImageType::RegionType & region) { + // Types + typedef typename ImageType::IndexType IndexType; + typedef typename ImageType::PointType PointType; + typedef typename ImageType::RegionType RegionType; + typedef typename ImageType::SizeType SizeType; + + // Region starting point + IndexType regionStart; + PointType start = bb->GetMinimum(); + image->TransformPhysicalPointToIndex(start, regionStart); + + // Region size + SizeType regionSize; + PointType maxs = bb->GetMaximum(); + PointType mins = bb->GetMinimum(); + for(unsigned int i=0; iGetSpacing()[i]); + } - // Create region - region.SetIndex(regionStart); - region.SetSize(regionSize); -} -//-------------------------------------------------------------------- - -//-------------------------------------------------------------------- -template -typename ImageType::Pointer -SetBackground(const ImageType * input, - const TMaskImageType * mask, - typename TMaskImageType::PixelType maskBG, - typename ImageType::PixelType outValue, - bool inPlace) { - typedef SetBackgroundImageFilter - SetBackgroundImageFilterType; - typename SetBackgroundImageFilterType::Pointer setBackgroundFilter - = SetBackgroundImageFilterType::New(); - // if (inPlace) setBackgroundFilter->ReleaseDataFlagOn(); // No seg fault - setBackgroundFilter->SetInPlace(inPlace); // This is important to keep memory low - setBackgroundFilter->SetInput(input); - setBackgroundFilter->SetInput2(mask); - setBackgroundFilter->SetMaskValue(maskBG); - setBackgroundFilter->SetOutsideValue(outValue); - setBackgroundFilter->Update(); - return setBackgroundFilter->GetOutput(); -} -//-------------------------------------------------------------------- - - -//-------------------------------------------------------------------- -template -int GetNumberOfConnectedComponentLabels(typename ImageType::Pointer input, - typename ImageType::PixelType BG, - bool isFullyConnected) { - // Connected Component label - typedef itk::ConnectedComponentImageFilter ConnectFilterType; - typename ConnectFilterType::Pointer connectFilter = ConnectFilterType::New(); - connectFilter->SetInput(input); - connectFilter->SetBackgroundValue(BG); - connectFilter->SetFullyConnected(isFullyConnected); - connectFilter->Update(); - - // Return result - return connectFilter->GetObjectCount(); -} -//-------------------------------------------------------------------- - -//-------------------------------------------------------------------- -/* - Warning : in this cas, we consider outputType like inputType, not - InternalImageType. Be sure it fits. - */ -template -typename ImageType::Pointer -Labelize(const ImageType * input, - typename ImageType::PixelType BG, - bool isFullyConnected, - int minimalComponentSize) { - // InternalImageType for storing large number of component - typedef itk::Image InternalImageType; - - // Connected Component label - typedef itk::ConnectedComponentImageFilter ConnectFilterType; - typename ConnectFilterType::Pointer connectFilter = ConnectFilterType::New(); - // connectFilter->ReleaseDataFlagOn(); - connectFilter->SetInput(input); - connectFilter->SetBackgroundValue(BG); - connectFilter->SetFullyConnected(isFullyConnected); + // Create region + region.SetIndex(regionStart); + region.SetSize(regionSize); + } + //-------------------------------------------------------------------- + + //-------------------------------------------------------------------- + template + typename ImageType::Pointer + SetBackground(const ImageType * input, + const TMaskImageType * mask, + typename TMaskImageType::PixelType maskBG, + typename ImageType::PixelType outValue, + bool inPlace) { + typedef SetBackgroundImageFilter + SetBackgroundImageFilterType; + typename SetBackgroundImageFilterType::Pointer setBackgroundFilter + = SetBackgroundImageFilterType::New(); + // if (inPlace) setBackgroundFilter->ReleaseDataFlagOn(); // No seg fault + setBackgroundFilter->SetInPlace(inPlace); // This is important to keep memory low + setBackgroundFilter->SetInput(input); + setBackgroundFilter->SetInput2(mask); + setBackgroundFilter->SetMaskValue(maskBG); + setBackgroundFilter->SetOutsideValue(outValue); + setBackgroundFilter->Update(); + return setBackgroundFilter->GetOutput(); + } + //-------------------------------------------------------------------- + + + //-------------------------------------------------------------------- + template + int GetNumberOfConnectedComponentLabels(const ImageType * input, + typename ImageType::PixelType BG, + bool isFullyConnected) { + // Connected Component label + typedef itk::ConnectedComponentImageFilter ConnectFilterType; + typename ConnectFilterType::Pointer connectFilter = ConnectFilterType::New(); + connectFilter->SetInput(input); + connectFilter->SetBackgroundValue(BG); + connectFilter->SetFullyConnected(isFullyConnected); + connectFilter->Update(); - // Sort by size and remove too small area. - typedef itk::RelabelComponentImageFilter RelabelFilterType; - typename RelabelFilterType::Pointer relabelFilter = RelabelFilterType::New(); - // relabelFilter->ReleaseDataFlagOn(); // if yes, fail when ExplosionControlledThresholdConnectedImageFilter ??? - relabelFilter->SetInput(connectFilter->GetOutput()); - relabelFilter->SetMinimumObjectSize(minimalComponentSize); - relabelFilter->Update(); - - // DD(relabelFilter->GetNumberOfObjects()); - // DD(relabelFilter->GetOriginalNumberOfObjects()); - // DD(relabelFilter->GetSizeOfObjectsInPhysicalUnits()[0]); - - // Return result - typename ImageType::Pointer output = relabelFilter->GetOutput(); - return output; -} -//-------------------------------------------------------------------- - - -//-------------------------------------------------------------------- -/* - Warning : in this cas, we consider outputType like inputType, not - InternalImageType. Be sure it fits. - */ -template -typename ImageType::Pointer -LabelizeAndCountNumberOfObjects(const ImageType * input, - typename ImageType::PixelType BG, - bool isFullyConnected, - int minimalComponentSize, - int & nb) { - // InternalImageType for storing large number of component - typedef itk::Image InternalImageType; + // Return result + return connectFilter->GetObjectCount(); + } + //-------------------------------------------------------------------- + + //-------------------------------------------------------------------- + /* + Warning : in this cas, we consider outputType like inputType, not + InternalImageType. Be sure it fits. + */ + template + typename ImageType::Pointer + Labelize(const ImageType * input, + typename ImageType::PixelType BG, + bool isFullyConnected, + int minimalComponentSize) { + // InternalImageType for storing large number of component + typedef itk::Image InternalImageType; - // Connected Component label - typedef itk::ConnectedComponentImageFilter ConnectFilterType; - typename ConnectFilterType::Pointer connectFilter = ConnectFilterType::New(); - // connectFilter->ReleaseDataFlagOn(); - connectFilter->SetInput(input); - connectFilter->SetBackgroundValue(BG); - connectFilter->SetFullyConnected(isFullyConnected); + // Connected Component label + typedef itk::ConnectedComponentImageFilter ConnectFilterType; + typename ConnectFilterType::Pointer connectFilter = ConnectFilterType::New(); + // connectFilter->ReleaseDataFlagOn(); + connectFilter->SetInput(input); + connectFilter->SetBackgroundValue(BG); + connectFilter->SetFullyConnected(isFullyConnected); - // Sort by size and remove too small area. - typedef itk::RelabelComponentImageFilter RelabelFilterType; - typename RelabelFilterType::Pointer relabelFilter = RelabelFilterType::New(); - // relabelFilter->ReleaseDataFlagOn(); // if yes, fail when ExplosionControlledThresholdConnectedImageFilter ??? - relabelFilter->SetInput(connectFilter->GetOutput()); - relabelFilter->SetMinimumObjectSize(minimalComponentSize); - relabelFilter->Update(); - - nb = relabelFilter->GetNumberOfObjects(); - // DD(relabelFilter->GetNumberOfObjects()); - // DD(relabelFilter->GetOriginalNumberOfObjects()); - // DD(relabelFilter->GetSizeOfObjectsInPhysicalUnits()[0]); - - // Return result - typename ImageType::Pointer output = relabelFilter->GetOutput(); - return output; -} -//-------------------------------------------------------------------- - - - -//-------------------------------------------------------------------- -template -typename ImageType::Pointer -RemoveLabels(typename ImageType::Pointer input, - typename ImageType::PixelType BG, - std::vector & labelsToRemove) { - typename ImageType::Pointer working_image = input; - for (unsigned int i=0; i SetBackgroundImageFilterType; - typename SetBackgroundImageFilterType::Pointer setBackgroundFilter = SetBackgroundImageFilterType::New(); - setBackgroundFilter->SetInput(input); - setBackgroundFilter->SetInput2(input); - setBackgroundFilter->SetMaskValue(labelsToRemove[i]); - setBackgroundFilter->SetOutsideValue(BG); - setBackgroundFilter->Update(); - working_image = setBackgroundFilter->GetOutput(); + // Sort by size and remove too small area. + typedef itk::RelabelComponentImageFilter RelabelFilterType; + typename RelabelFilterType::Pointer relabelFilter = RelabelFilterType::New(); + // relabelFilter->ReleaseDataFlagOn(); // if yes, fail when ExplosionControlledThresholdConnectedImageFilter ??? + relabelFilter->SetInput(connectFilter->GetOutput()); + relabelFilter->SetMinimumObjectSize(minimalComponentSize); + relabelFilter->Update(); + + // Return result + typename ImageType::Pointer output = relabelFilter->GetOutput(); + return output; } - return working_image; -} -//-------------------------------------------------------------------- + //-------------------------------------------------------------------- -//-------------------------------------------------------------------- -template -typename ImageType::Pointer -KeepLabels(const ImageType * input, - typename ImageType::PixelType BG, - typename ImageType::PixelType FG, - typename ImageType::PixelType firstKeep, - typename ImageType::PixelType lastKeep, - bool useLastKeep) { - typedef itk::BinaryThresholdImageFilter BinarizeFilterType; - typename BinarizeFilterType::Pointer binarizeFilter = BinarizeFilterType::New(); - binarizeFilter->SetInput(input); - binarizeFilter->SetLowerThreshold(firstKeep); - if (useLastKeep) binarizeFilter->SetUpperThreshold(lastKeep); - binarizeFilter->SetInsideValue(FG); - binarizeFilter->SetOutsideValue(BG); - binarizeFilter->Update(); - return binarizeFilter->GetOutput(); -} -//-------------------------------------------------------------------- - - -//-------------------------------------------------------------------- -template -typename ImageType::Pointer -LabelizeAndSelectLabels(typename ImageType::Pointer input, - typename ImageType::PixelType BG, - typename ImageType::PixelType FG, - bool isFullyConnected, - int minimalComponentSize, - LabelizeParameters * param) -{ - typename ImageType::Pointer working_image; - working_image = Labelize(input, BG, isFullyConnected, minimalComponentSize); - working_image = RemoveLabels(working_image, BG, param->GetLabelsToRemove()); - working_image = KeepLabels(working_image, - BG, FG, - param->GetFirstKeep(), - param->GetLastKeep(), - param->GetUseLastKeep()); - return working_image; -} -//-------------------------------------------------------------------- - - -//-------------------------------------------------------------------- -template -typename ImageType::Pointer -ResizeImageLike(typename ImageType::Pointer input, - typename ImageType::Pointer like, - typename ImageType::PixelType backgroundValue) -{ - typedef CropLikeImageFilter CropFilterType; - typename CropFilterType::Pointer cropFilter = CropFilterType::New(); - cropFilter->SetInput(input); - cropFilter->SetCropLikeImage(like); - cropFilter->SetBackgroundValue(backgroundValue); - cropFilter->Update(); - return cropFilter->GetOutput(); -} -//-------------------------------------------------------------------- - - -//-------------------------------------------------------------------- -template -typename MaskImageType::Pointer -SliceBySliceRelativePosition(const MaskImageType * input, - const MaskImageType * object, - int direction, - double threshold, - std::string orientation, - bool uniqueConnectedComponent, - double spacing, - bool inverseflag) -{ - typedef SliceBySliceRelativePositionFilter SliceRelPosFilterType; - typename SliceRelPosFilterType::Pointer sliceRelPosFilter = SliceRelPosFilterType::New(); - sliceRelPosFilter->VerboseStepFlagOff(); - sliceRelPosFilter->WriteStepFlagOff(); - sliceRelPosFilter->SetInput(input); - sliceRelPosFilter->SetInputObject(object); - sliceRelPosFilter->SetDirection(direction); - sliceRelPosFilter->SetFuzzyThreshold(threshold); - sliceRelPosFilter->AddOrientationTypeString(orientation); - sliceRelPosFilter->SetResampleBeforeRelativePositionFilter((spacing != -1)); - sliceRelPosFilter->SetIntermediateSpacing(spacing); - sliceRelPosFilter->SetUniqueConnectedComponentBySlice(uniqueConnectedComponent); - sliceRelPosFilter->SetInverseOrientationFlag(inverseflag); - // sliceRelPosFilter->SetAutoCropFlag(true); ?? - sliceRelPosFilter->Update(); - return sliceRelPosFilter->GetOutput(); -} -//-------------------------------------------------------------------- - -//-------------------------------------------------------------------- -template -bool -FindExtremaPointInAGivenDirection(const ImageType * input, - typename ImageType::PixelType bg, - int direction, bool opposite, - typename ImageType::PointType & point) -{ - typename ImageType::PointType dummy; - return FindExtremaPointInAGivenDirection(input, bg, direction, - opposite, dummy, 0, point); -} -//-------------------------------------------------------------------- - -//-------------------------------------------------------------------- -template -bool -FindExtremaPointInAGivenDirection(const ImageType * input, - typename ImageType::PixelType bg, - int direction, bool opposite, - typename ImageType::PointType refpoint, - double distanceMax, - typename ImageType::PointType & point) -{ + //-------------------------------------------------------------------- /* - loop over input pixels, store the index in the fg that is max - according to the given direction. + Warning : in this cas, we consider outputType like inputType, not + InternalImageType. Be sure it fits. */ - typedef itk::ImageRegionConstIteratorWithIndex IteratorType; - IteratorType iter(input, input->GetLargestPossibleRegion()); - iter.GoToBegin(); - typename ImageType::IndexType max = input->GetLargestPossibleRegion().GetIndex(); - if (opposite) max = max+input->GetLargestPossibleRegion().GetSize(); - bool found=false; - while (!iter.IsAtEnd()) { - if (iter.Get() != bg) { - bool test = iter.GetIndex()[direction] > max[direction]; - if (opposite) test = !test; - if (test) { - typename ImageType::PointType p; - input->TransformIndexToPhysicalPoint(iter.GetIndex(), p); - if ((distanceMax==0) || (p.EuclideanDistanceTo(refpoint) < distanceMax)) { - max = iter.GetIndex(); - found = true; + template + typename ImageType::Pointer + LabelizeAndCountNumberOfObjects(const ImageType * input, + typename ImageType::PixelType BG, + bool isFullyConnected, + int minimalComponentSize, + int & nb) { + // InternalImageType for storing large number of component + typedef itk::Image InternalImageType; + + // Connected Component label + typedef itk::ConnectedComponentImageFilter ConnectFilterType; + typename ConnectFilterType::Pointer connectFilter = ConnectFilterType::New(); + // connectFilter->ReleaseDataFlagOn(); + connectFilter->SetInput(input); + connectFilter->SetBackgroundValue(BG); + connectFilter->SetFullyConnected(isFullyConnected); + + // Sort by size and remove too small area. + typedef itk::RelabelComponentImageFilter RelabelFilterType; + typename RelabelFilterType::Pointer relabelFilter = RelabelFilterType::New(); + // relabelFilter->ReleaseDataFlagOn(); // if yes, fail when ExplosionControlledThresholdConnectedImageFilter ??? + relabelFilter->SetInput(connectFilter->GetOutput()); + relabelFilter->SetMinimumObjectSize(minimalComponentSize); + relabelFilter->Update(); + + nb = relabelFilter->GetNumberOfObjects(); + // DD(relabelFilter->GetOriginalNumberOfObjects()); + // DD(relabelFilter->GetSizeOfObjectsInPhysicalUnits()[0]); + + // Return result + typename ImageType::Pointer output = relabelFilter->GetOutput(); + return output; + } + //-------------------------------------------------------------------- + + + + //-------------------------------------------------------------------- + template + typename ImageType::Pointer + RemoveLabels(const ImageType * input, + typename ImageType::PixelType BG, + std::vector & labelsToRemove) { + assert(labelsToRemove.size() != 0); + typename ImageType::Pointer working_image;// = input; + for (unsigned int i=0; i SetBackgroundImageFilterType; + typename SetBackgroundImageFilterType::Pointer setBackgroundFilter = SetBackgroundImageFilterType::New(); + setBackgroundFilter->SetInput(input); + setBackgroundFilter->SetInput2(input); + setBackgroundFilter->SetMaskValue(labelsToRemove[i]); + setBackgroundFilter->SetOutsideValue(BG); + setBackgroundFilter->Update(); + working_image = setBackgroundFilter->GetOutput(); + } + return working_image; + } + //-------------------------------------------------------------------- + + + //-------------------------------------------------------------------- + template + typename ImageType::Pointer + KeepLabels(const ImageType * input, + typename ImageType::PixelType BG, + typename ImageType::PixelType FG, + typename ImageType::PixelType firstKeep, + typename ImageType::PixelType lastKeep, + bool useLastKeep) { + typedef itk::BinaryThresholdImageFilter BinarizeFilterType; + typename BinarizeFilterType::Pointer binarizeFilter = BinarizeFilterType::New(); + binarizeFilter->SetInput(input); + binarizeFilter->SetLowerThreshold(firstKeep); + if (useLastKeep) binarizeFilter->SetUpperThreshold(lastKeep); + binarizeFilter->SetInsideValue(FG); + binarizeFilter->SetOutsideValue(BG); + binarizeFilter->Update(); + return binarizeFilter->GetOutput(); + } + //-------------------------------------------------------------------- + + + //-------------------------------------------------------------------- + template + typename ImageType::Pointer + LabelizeAndSelectLabels(const ImageType * input, + typename ImageType::PixelType BG, + typename ImageType::PixelType FG, + bool isFullyConnected, + int minimalComponentSize, + LabelizeParameters * param) + { + typename ImageType::Pointer working_image; + working_image = Labelize(input, BG, isFullyConnected, minimalComponentSize); + working_image = RemoveLabels(working_image, BG, param->GetLabelsToRemove()); + working_image = KeepLabels(working_image, + BG, FG, + param->GetFirstKeep(), + param->GetLastKeep(), + param->GetUseLastKeep()); + return working_image; + } + //-------------------------------------------------------------------- + + + //-------------------------------------------------------------------- + template + typename ImageType::Pointer + ResizeImageLike(const ImageType * input, + const itk::ImageBase * like, + typename ImageType::PixelType backgroundValue) + { + typedef CropLikeImageFilter CropFilterType; + typename CropFilterType::Pointer cropFilter = CropFilterType::New(); + cropFilter->SetInput(input); + cropFilter->SetCropLikeImage(like); + cropFilter->SetBackgroundValue(backgroundValue); + cropFilter->Update(); + return cropFilter->GetOutput(); + } + //-------------------------------------------------------------------- + + + //-------------------------------------------------------------------- + template + typename MaskImageType::Pointer + SliceBySliceRelativePosition(const MaskImageType * input, + const MaskImageType * object, + int direction, + double threshold, + std::string orientation, + bool uniqueConnectedComponent, + double spacing, + bool inverseflag) + { + typedef SliceBySliceRelativePositionFilter SliceRelPosFilterType; + typename SliceRelPosFilterType::Pointer sliceRelPosFilter = SliceRelPosFilterType::New(); + sliceRelPosFilter->VerboseStepFlagOff(); + sliceRelPosFilter->WriteStepFlagOff(); + sliceRelPosFilter->SetInput(input); + sliceRelPosFilter->SetInputObject(object); + sliceRelPosFilter->SetDirection(direction); + sliceRelPosFilter->SetFuzzyThreshold(threshold); + sliceRelPosFilter->AddOrientationTypeString(orientation); + sliceRelPosFilter->SetIntermediateSpacingFlag((spacing != -1)); + sliceRelPosFilter->SetIntermediateSpacing(spacing); + sliceRelPosFilter->SetUniqueConnectedComponentBySlice(uniqueConnectedComponent); + sliceRelPosFilter->SetInverseOrientationFlag(inverseflag); + // sliceRelPosFilter->SetAutoCropFlag(true); ?? + sliceRelPosFilter->Update(); + return sliceRelPosFilter->GetOutput(); + } + //-------------------------------------------------------------------- + + //-------------------------------------------------------------------- + template + bool + FindExtremaPointInAGivenDirection(const ImageType * input, + typename ImageType::PixelType bg, + int direction, bool opposite, + typename ImageType::PointType & point) + { + typename ImageType::PointType dummy; + return FindExtremaPointInAGivenDirection(input, bg, direction, + opposite, dummy, 0, point); + } + //-------------------------------------------------------------------- + + //-------------------------------------------------------------------- + template + bool + FindExtremaPointInAGivenDirection(const ImageType * input, + typename ImageType::PixelType bg, + int direction, bool opposite, + typename ImageType::PointType refpoint, + double distanceMax, + typename ImageType::PointType & point) + { + /* + loop over input pixels, store the index in the fg that is max + according to the given direction. + */ + typedef itk::ImageRegionConstIteratorWithIndex IteratorType; + IteratorType iter(input, input->GetLargestPossibleRegion()); + iter.GoToBegin(); + typename ImageType::IndexType max = input->GetLargestPossibleRegion().GetIndex(); + if (opposite) max = max+input->GetLargestPossibleRegion().GetSize(); + bool found=false; + while (!iter.IsAtEnd()) { + if (iter.Get() != bg) { + bool test = iter.GetIndex()[direction] > max[direction]; + if (opposite) test = !test; + if (test) { + typename ImageType::PointType p; + input->TransformIndexToPhysicalPoint(iter.GetIndex(), p); + if ((distanceMax==0) || (p.EuclideanDistanceTo(refpoint) < distanceMax)) { + max = iter.GetIndex(); + found = true; + } } } + ++iter; } - ++iter; - } - if (!found) return false; - input->TransformIndexToPhysicalPoint(max, point); - return true; -} -//-------------------------------------------------------------------- - - -//-------------------------------------------------------------------- -template -typename ImageType::Pointer -CropImageAbove(typename ImageType::Pointer image, - int dim, double min, - bool autoCrop, - typename ImageType::PixelType BG) -{ - return CropImageAlongOneAxis(image, dim, - image->GetOrigin()[dim], - min, - autoCrop, BG); -} -//-------------------------------------------------------------------- - - -//-------------------------------------------------------------------- -template -typename ImageType::Pointer -CropImageBelow(typename ImageType::Pointer image, - int dim, double max, - bool autoCrop, - typename ImageType::PixelType BG) -{ - typename ImageType::PointType p; - image->TransformIndexToPhysicalPoint(image->GetLargestPossibleRegion().GetIndex()+ - image->GetLargestPossibleRegion().GetSize(), p); - return CropImageAlongOneAxis(image, dim, max, p[dim], autoCrop, BG); -} -//-------------------------------------------------------------------- - - -//-------------------------------------------------------------------- -template -typename ImageType::Pointer -CropImageAlongOneAxis(typename ImageType::Pointer image, - int dim, double min, double max, - bool autoCrop, - typename ImageType::PixelType BG) -{ - // Compute region size - typename ImageType::RegionType region; - typename ImageType::SizeType size = image->GetLargestPossibleRegion().GetSize(); - typename ImageType::PointType p = image->GetOrigin(); - p[dim] = min; - typename ImageType::IndexType start; - image->TransformPhysicalPointToIndex(p, start); - p[dim] = max; - typename ImageType::IndexType end; - image->TransformPhysicalPointToIndex(p, end); - size[dim] = fabs(end[dim]-start[dim]); - region.SetIndex(start); - region.SetSize(size); + if (!found) return false; + input->TransformIndexToPhysicalPoint(max, point); + return true; + } + //-------------------------------------------------------------------- + + + //-------------------------------------------------------------------- + template + typename ImageType::Pointer + CropImageAbove(const ImageType * image, + int dim, double min, bool autoCrop, + typename ImageType::PixelType BG) + { + return CropImageAlongOneAxis(image, dim, + image->GetOrigin()[dim], + min, + autoCrop, BG); + } + //-------------------------------------------------------------------- + + + //-------------------------------------------------------------------- + template + typename ImageType::Pointer + CropImageBelow(const ImageType * image, + int dim, double max, bool autoCrop, + typename ImageType::PixelType BG) + { + typename ImageType::PointType p; + image->TransformIndexToPhysicalPoint(image->GetLargestPossibleRegion().GetIndex()+ + image->GetLargestPossibleRegion().GetSize(), p); + return CropImageAlongOneAxis(image, dim, max, p[dim], autoCrop, BG); + } + //-------------------------------------------------------------------- + + + //-------------------------------------------------------------------- + template + typename ImageType::Pointer + CropImageAlongOneAxis(const ImageType * image, + int dim, double min, double max, + bool autoCrop, typename ImageType::PixelType BG) + { + // Compute region size + typename ImageType::RegionType region; + typename ImageType::SizeType size = image->GetLargestPossibleRegion().GetSize(); + typename ImageType::PointType p = image->GetOrigin(); + p[dim] = min; + typename ImageType::IndexType start; + image->TransformPhysicalPointToIndex(p, start); + p[dim] = max; + typename ImageType::IndexType end; + image->TransformPhysicalPointToIndex(p, end); + size[dim] = fabs(end[dim]-start[dim]); + region.SetIndex(start); + region.SetSize(size); - // Perform Crop - typedef itk::RegionOfInterestImageFilter CropFilterType; - typename CropFilterType::Pointer cropFilter = CropFilterType::New(); - cropFilter->SetInput(image); - cropFilter->SetRegionOfInterest(region); - cropFilter->Update(); - typename ImageType::Pointer result = cropFilter->GetOutput(); + // Perform Crop + typedef itk::RegionOfInterestImageFilter CropFilterType; + typename CropFilterType::Pointer cropFilter = CropFilterType::New(); + cropFilter->SetInput(image); + cropFilter->SetRegionOfInterest(region); + cropFilter->Update(); + typename ImageType::Pointer result = cropFilter->GetOutput(); - // Auto Crop - if (autoCrop) { - result = AutoCrop(result, BG); - } - return result; -} -//-------------------------------------------------------------------- - - -//-------------------------------------------------------------------- -template -void -ComputeCentroids(typename ImageType::Pointer image, - typename ImageType::PixelType BG, - std::vector & centroids) -{ - typedef long LabelType; - static const unsigned int Dim = ImageType::ImageDimension; - typedef itk::ShapeLabelObject< LabelType, Dim > LabelObjectType; - typedef itk::LabelMap< LabelObjectType > LabelMapType; - typedef itk::LabelImageToLabelMapFilter ImageToMapFilterType; - typename ImageToMapFilterType::Pointer imageToLabelFilter = ImageToMapFilterType::New(); - typedef itk::ShapeLabelMapFilter ShapeFilterType; - typename ShapeFilterType::Pointer statFilter = ShapeFilterType::New(); - imageToLabelFilter->SetBackgroundValue(BG); - imageToLabelFilter->SetInput(image); - statFilter->SetInput(imageToLabelFilter->GetOutput()); - statFilter->Update(); - typename LabelMapType::Pointer labelMap = statFilter->GetOutput(); - - centroids.clear(); - typename ImageType::PointType dummy; - centroids.push_back(dummy); // label 0 -> no centroid, use dummy point - for(uint i=1; iGetNumberOfLabelObjects()+1; i++) { - centroids.push_back(labelMap->GetLabelObject(i)->GetCentroid()); - } -} -//-------------------------------------------------------------------- - - -//-------------------------------------------------------------------- -template -void -ExtractSlices(typename ImageType::Pointer image, - int direction, - std::vector::Pointer > & slices) -{ - typedef ExtractSliceFilter ExtractSliceFilterType; - typedef typename ExtractSliceFilterType::SliceType SliceType; - typename ExtractSliceFilterType::Pointer - extractSliceFilter = ExtractSliceFilterType::New(); - extractSliceFilter->SetInput(image); - extractSliceFilter->SetDirection(direction); - extractSliceFilter->Update(); - extractSliceFilter->GetOutputSlices(slices); -} -//-------------------------------------------------------------------- - - -//-------------------------------------------------------------------- -template -typename ImageType::Pointer -JoinSlices(std::vector::Pointer > & slices, - typename ImageType::Pointer input, - int direction) { - typedef typename itk::Image SliceType; - typedef itk::JoinSeriesImageFilter JoinSeriesFilterType; - typename JoinSeriesFilterType::Pointer joinFilter = JoinSeriesFilterType::New(); - joinFilter->SetOrigin(input->GetOrigin()[direction]); - joinFilter->SetSpacing(input->GetSpacing()[direction]); - for(unsigned int i=0; iPushBackInput(slices[i]); - } - joinFilter->Update(); - return joinFilter->GetOutput(); -} -//-------------------------------------------------------------------- - - -//-------------------------------------------------------------------- -template -void -PointsUtils::Convert2DTo3D(const PointType2D & p, - ImagePointer image, - const int slice, - PointType3D & p3D) -{ - p3D[0] = p[0]; - p3D[1] = p[1]; - p3D[2] = (image->GetLargestPossibleRegion().GetIndex()[2]+slice)*image->GetSpacing()[2] - + image->GetOrigin()[2]; -} -//-------------------------------------------------------------------- - - -//-------------------------------------------------------------------- -template -void -PointsUtils::Convert2DTo3DList(const MapPoint2DType & map, - ImagePointer image, - VectorPoint3DType & list) -{ - typename MapPoint2DType::const_iterator iter = map.begin(); - while (iter != map.end()) { - PointType3D p; - Convert2DTo3D(iter->second, image, iter->first, p); - list.push_back(p); - ++iter; - } -} -//-------------------------------------------------------------------- - -//-------------------------------------------------------------------- -template -void -WriteListOfLandmarks(std::vector points, - std::string filename) -{ - std::ofstream os; - openFileForWriting(os, filename); - os << "LANDMARKS1" << std::endl; - for(uint i=0; i -typename ImageType::Pointer -Dilate(typename ImageType::Pointer image, - double radiusInMM, - typename ImageType::PixelType BG, - typename ImageType::PixelType FG, - bool extendSupport) -{ - typename ImageType::SizeType r; - for(uint i=0; iGetSpacing()[i]); - return Dilate(image, r, BG, FG, extendSupport); -} -//-------------------------------------------------------------------- - - -//-------------------------------------------------------------------- -template -typename ImageType::Pointer -Dilate(typename ImageType::Pointer image, - typename ImageType::PointType radiusInMM, - typename ImageType::PixelType BG, - typename ImageType::PixelType FG, - bool extendSupport) -{ - typename ImageType::SizeType r; - for(uint i=0; iGetSpacing()[i]); - return Dilate(image, r, BG, FG, extendSupport); -} -//-------------------------------------------------------------------- - - -//-------------------------------------------------------------------- -template -typename ImageType::Pointer -Dilate(typename ImageType::Pointer image, - typename ImageType::SizeType radius, - typename ImageType::PixelType BG, - typename ImageType::PixelType FG, - bool extendSupport) -{ - // Create kernel for dilatation - typedef itk::BinaryBallStructuringElement KernelType; - KernelType structuringElement; - structuringElement.SetRadius(radius); - structuringElement.CreateStructuringElement(); - - if (extendSupport) { - typedef itk::ConstantPadImageFilter PadFilterType; - typename PadFilterType::Pointer padFilter = PadFilterType::New(); - padFilter->SetInput(image); - typename ImageType::SizeType lower; - typename ImageType::SizeType upper; - for(uint i=0; i<3; i++) { - lower[i] = upper[i] = 2*(radius[i]+1); + // Auto Crop + if (autoCrop) { + result = AutoCrop(result, BG); + } + return result; + } + //-------------------------------------------------------------------- + + + //-------------------------------------------------------------------- + template + void + ComputeCentroids(const ImageType * image, + typename ImageType::PixelType BG, + std::vector & centroids) + { + typedef long LabelType; + static const unsigned int Dim = ImageType::ImageDimension; + typedef itk::ShapeLabelObject< LabelType, Dim > LabelObjectType; + typedef itk::LabelMap< LabelObjectType > LabelMapType; + typedef itk::LabelImageToLabelMapFilter ImageToMapFilterType; + typename ImageToMapFilterType::Pointer imageToLabelFilter = ImageToMapFilterType::New(); + typedef itk::ShapeLabelMapFilter ShapeFilterType; + typename ShapeFilterType::Pointer statFilter = ShapeFilterType::New(); + imageToLabelFilter->SetBackgroundValue(BG); + imageToLabelFilter->SetInput(image); + statFilter->SetInput(imageToLabelFilter->GetOutput()); + statFilter->Update(); + typename LabelMapType::Pointer labelMap = statFilter->GetOutput(); + + centroids.clear(); + typename ImageType::PointType dummy; + centroids.push_back(dummy); // label 0 -> no centroid, use dummy point + for(uint i=1; iGetNumberOfLabelObjects()+1; i++) { + centroids.push_back(labelMap->GetLabelObject(i)->GetCentroid()); + } + } + //-------------------------------------------------------------------- + + + //-------------------------------------------------------------------- + template + void + ExtractSlices(const ImageType * image, int direction, + std::vector::Pointer > & slices) + { + typedef ExtractSliceFilter ExtractSliceFilterType; + typedef typename ExtractSliceFilterType::SliceType SliceType; + typename ExtractSliceFilterType::Pointer + extractSliceFilter = ExtractSliceFilterType::New(); + extractSliceFilter->SetInput(image); + extractSliceFilter->SetDirection(direction); + extractSliceFilter->Update(); + extractSliceFilter->GetOutputSlices(slices); + } + //-------------------------------------------------------------------- + + + //-------------------------------------------------------------------- + template + typename ImageType::Pointer + JoinSlices(std::vector::Pointer > & slices, + const ImageType * input, + int direction) { + typedef typename itk::Image SliceType; + typedef itk::JoinSeriesImageFilter JoinSeriesFilterType; + typename JoinSeriesFilterType::Pointer joinFilter = JoinSeriesFilterType::New(); + joinFilter->SetOrigin(input->GetOrigin()[direction]); + joinFilter->SetSpacing(input->GetSpacing()[direction]); + for(unsigned int i=0; iPushBackInput(slices[i]); + } + joinFilter->Update(); + return joinFilter->GetOutput(); + } + //-------------------------------------------------------------------- + + + //-------------------------------------------------------------------- + template + void + PointsUtils::Convert2DTo3D(const PointType2D & p2D, + const ImageType * image, + const int slice, + PointType3D & p3D) + { + IndexType3D index3D; + index3D[0] = index3D[1] = 0; + index3D[2] = image->GetLargestPossibleRegion().GetIndex()[2]+slice; + image->TransformIndexToPhysicalPoint(index3D, p3D); + p3D[0] = p2D[0]; + p3D[1] = p2D[1]; + // p3D[2] = p[2];//(image->GetLargestPossibleRegion().GetIndex()[2]+slice)*image->GetSpacing()[2] + // + image->GetOrigin()[2]; + } + //-------------------------------------------------------------------- + + + //-------------------------------------------------------------------- + template + void + PointsUtils::Convert2DTo3DList(const MapPoint2DType & map, + const ImageType * image, + VectorPoint3DType & list) + { + typename MapPoint2DType::const_iterator iter = map.begin(); + while (iter != map.end()) { + PointType3D p; + Convert2DTo3D(iter->second, image, iter->first, p); + list.push_back(p); + ++iter; } - padFilter->SetPadLowerBound(lower); - padFilter->SetPadUpperBound(upper); - padFilter->Update(); - image = padFilter->GetOutput(); - } - - // Dilate filter - typedef itk::BinaryDilateImageFilter DilateFilterType; - typename DilateFilterType::Pointer dilateFilter = DilateFilterType::New(); - dilateFilter->SetBackgroundValue(BG); - dilateFilter->SetForegroundValue(FG); - dilateFilter->SetBoundaryToForeground(false); - dilateFilter->SetKernel(structuringElement); - dilateFilter->SetInput(image); - dilateFilter->Update(); - return image = dilateFilter->GetOutput(); -} -//-------------------------------------------------------------------- - - -//-------------------------------------------------------------------- -template -void ConvertOption(std::string optionName, uint given, - ValueType * values, VectorType & p, - uint dim, bool required) -{ - if (required && (given == 0)) { - clitkExceptionMacro("The option --" << optionName << " must be set and have 1 or " + } + //-------------------------------------------------------------------- + + //-------------------------------------------------------------------- + template + void + WriteListOfLandmarks(std::vector points, + std::string filename) + { + std::ofstream os; + openFileForWriting(os, filename); + os << "LANDMARKS1" << std::endl; + for(uint i=0; i + typename ImageType::Pointer + Dilate(const ImageType * image, double radiusInMM, + typename ImageType::PixelType BG, + typename ImageType::PixelType FG, + bool extendSupport) + { + typename ImageType::SizeType r; + for(uint i=0; iGetSpacing()[i]); + return Dilate(image, r, BG, FG, extendSupport); + } + //-------------------------------------------------------------------- + + + //-------------------------------------------------------------------- + template + typename ImageType::Pointer + Dilate(const ImageType * image, typename ImageType::PointType radiusInMM, + typename ImageType::PixelType BG, + typename ImageType::PixelType FG, + bool extendSupport) + { + typename ImageType::SizeType r; + for(uint i=0; iGetSpacing()[i]); + return Dilate(image, r, BG, FG, extendSupport); + } + //-------------------------------------------------------------------- + + + //-------------------------------------------------------------------- + template + typename ImageType::Pointer + Dilate(const ImageType * image, typename ImageType::SizeType radius, + typename ImageType::PixelType BG, + typename ImageType::PixelType FG, + bool extendSupport) + { + // Create kernel for dilatation + typedef itk::BinaryBallStructuringElement KernelType; + KernelType structuringElement; + structuringElement.SetRadius(radius); + structuringElement.CreateStructuringElement(); + + typename ImageType::Pointer output; + if (extendSupport) { + typedef itk::ConstantPadImageFilter PadFilterType; + typename PadFilterType::Pointer padFilter = PadFilterType::New(); + padFilter->SetInput(image); + typename ImageType::SizeType lower; + typename ImageType::SizeType upper; + for(uint i=0; i<3; i++) { + lower[i] = upper[i] = 2*(radius[i]+1); + } + padFilter->SetPadLowerBound(lower); + padFilter->SetPadUpperBound(upper); + padFilter->Update(); + output = padFilter->GetOutput(); + } + + // Dilate filter + typedef itk::BinaryDilateImageFilter DilateFilterType; + typename DilateFilterType::Pointer dilateFilter = DilateFilterType::New(); + dilateFilter->SetBackgroundValue(BG); + dilateFilter->SetForegroundValue(FG); + dilateFilter->SetBoundaryToForeground(false); + dilateFilter->SetKernel(structuringElement); + dilateFilter->SetInput(output); + dilateFilter->Update(); + return dilateFilter->GetOutput(); + } + //-------------------------------------------------------------------- + + + //-------------------------------------------------------------------- + template + void ConvertOption(std::string optionName, uint given, + ValueType * values, VectorType & p, + uint dim, bool required) + { + if (required && (given == 0)) { + clitkExceptionMacro("The option --" << optionName << " must be set and have 1 or " + << dim << " values."); + } + if (given == 1) { + for(uint i=0; i to accelerate, start with formula, when change sign -> stop and fill -*/ -template -void -SliceBySliceSetBackgroundFromLineSeparation(typename ImageType::Pointer input, - std::vector & lA, - std::vector & lB, - typename ImageType::PixelType BG, - int mainDirection, - double offsetToKeep) -{ + //-------------------------------------------------------------------- + + + //-------------------------------------------------------------------- + /* + http://www.gamedev.net/community/forums/topic.asp?topic_id=542870 + Assuming the points are (Ax,Ay) (Bx,By) and (Cx,Cy), you need to compute: + (Bx - Ax) * (Cy - Ay) - (By - Ay) * (Cx - Ax) + This will equal zero if the point C is on the line formed by + points A and B, and will have a different sign depending on the + side. Which side this is depends on the orientation of your (x,y) + coordinates, but you can plug test values for A,B and C into this + formula to determine whether negative values are to the left or to + the right. + => to accelerate, start with formula, when change sign -> stop and fill + + offsetToKeep = is used to determine which side of the line we + keep. The point along the mainDirection but 'offsetToKeep' mm away + is kept. - typedef itk::ImageSliceIteratorWithIndex SliceIteratorType; - SliceIteratorType siter = SliceIteratorType(input, - input->GetLargestPossibleRegion()); - siter.SetFirstDirection(0); - siter.SetSecondDirection(1); - siter.GoToBegin(); - int i=0; - typename ImageType::PointType A; - typename ImageType::PointType B; - typename ImageType::PointType C; - while (!siter.IsAtEnd()) { - // Check that the current slice correspond to the current point - input->TransformIndexToPhysicalPoint(siter.GetIndex(), C); - if (C[2] != lA[i][2]) { + */ + template + void + SliceBySliceSetBackgroundFromLineSeparation(ImageType * input, + std::vector & lA, + std::vector & lB, + typename ImageType::PixelType BG, + int mainDirection, + double offsetToKeep) + { + typedef itk::ImageSliceIteratorWithIndex SliceIteratorType; + SliceIteratorType siter = SliceIteratorType(input, + input->GetLargestPossibleRegion()); + siter.SetFirstDirection(0); + siter.SetSecondDirection(1); + siter.GoToBegin(); + uint i=0; + typename ImageType::PointType A; + typename ImageType::PointType B; + typename ImageType::PointType C; + assert(lA.size() == lB.size()); + while ((iTransformIndexToPhysicalPoint(siter.GetIndex(), C); // DD(C); + // DD(i); // DD(lA[i]); - } - else { - // Define A,B,C points - A = lA[i]; - B = lB[i]; - C = A; - C[mainDirection] += offsetToKeep; // I know I must keep this point - double s = (B[0] - A[0]) * (C[1] - A[1]) - (B[1] - A[1]) * (C[0] - A[0]); - bool isPositive = s<0; - while (!siter.IsAtEndOfSlice()) { - while (!siter.IsAtEndOfLine()) { - // Very slow, I know ... but image should be very small - input->TransformIndexToPhysicalPoint(siter.GetIndex(), C); - double s = (B[0] - A[0]) * (C[1] - A[1]) - (B[1] - A[1]) * (C[0] - A[0]); - if (s == 0) siter.Set(BG); // on the line, we decide to remove - if (isPositive) { - if (s > 0) siter.Set(BG); - } - else { - if (s < 0) siter.Set(BG); - } - ++siter; - } - siter.NextLine(); + if ((fabs(C[2] - lA[i][2]))>0.01) { // is !equal with a tolerance of 0.01 mm } - ++i; + else { + // Define A,B,C points + A = lA[i]; + B = lB[i]; + C = A; + // DD(A); + // DD(B); + // DD(C); + + // Check that the line is not a point (A=B) + bool p = (A[0] == B[0]) && (A[1] == B[1]); + + if (!p) { + C[mainDirection] += offsetToKeep; // I know I must keep this point + double s = (B[0] - A[0]) * (C[1] - A[1]) - (B[1] - A[1]) * (C[0] - A[0]); + bool isPositive = s<0; + while (!siter.IsAtEndOfSlice()) { + while (!siter.IsAtEndOfLine()) { + // Very slow, I know ... but image should be very small + input->TransformIndexToPhysicalPoint(siter.GetIndex(), C); + double s = (B[0] - A[0]) * (C[1] - A[1]) - (B[1] - A[1]) * (C[0] - A[0]); + if (s == 0) siter.Set(BG); // on the line, we decide to remove + if (isPositive) { + if (s > 0) siter.Set(BG); + } + else { + if (s < 0) siter.Set(BG); + } + ++siter; + } + siter.NextLine(); + } // end loop slice + } + + ++i; + } // End of current slice + siter.NextSlice(); + } + } + //-------------------------------------------------------------------- + + //-------------------------------------------------------------------- + template + void + AndNot(ImageType * input, + const ImageType * object, + typename ImageType::PixelType BG) + { + typename ImageType::Pointer o; + bool resized=false; + if (!clitk::HaveSameSizeAndSpacing(input, object)) { + o = clitk::ResizeImageLike(object, input, BG); + resized = true; + } + + typedef clitk::BooleanOperatorLabelImageFilter BoolFilterType; + typename BoolFilterType::Pointer boolFilter = BoolFilterType::New(); + boolFilter->InPlaceOn(); + boolFilter->SetInput1(input); + if (resized) boolFilter->SetInput2(o); + else boolFilter->SetInput2(object); + boolFilter->SetBackgroundValue1(BG); + boolFilter->SetBackgroundValue2(BG); + boolFilter->SetOperationType(BoolFilterType::AndNot); + boolFilter->Update(); + } + //-------------------------------------------------------------------- + + + //-------------------------------------------------------------------- + template + typename ImageType::Pointer + Binarize(const ImageType * input, + typename ImageType::PixelType lower, + typename ImageType::PixelType upper, + typename ImageType::PixelType BG=0, + typename ImageType::PixelType FG=1) + { + typedef itk::BinaryThresholdImageFilter BinaryThresholdFilterType; + typename BinaryThresholdFilterType::Pointer binarizeFilter = BinaryThresholdFilterType::New(); + binarizeFilter->SetInput(input); + binarizeFilter->SetLowerThreshold(lower); + binarizeFilter->SetUpperThreshold(upper); + binarizeFilter->SetInsideValue(FG); + binarizeFilter->SetOutsideValue(BG); + binarizeFilter->Update(); + return binarizeFilter->GetOutput(); + } + //-------------------------------------------------------------------- + + + //-------------------------------------------------------------------- + template + void + GetMinMaxPointPosition(const ImageType * input, + typename ImageType::PointType & min, + typename ImageType::PointType & max) + { + typename ImageType::IndexType index = input->GetLargestPossibleRegion().GetIndex(); + input->TransformIndexToPhysicalPoint(index, min); + index = index+input->GetLargestPossibleRegion().GetSize(); + input->TransformIndexToPhysicalPoint(index, max); + } + //-------------------------------------------------------------------- + + + //-------------------------------------------------------------------- + template + typename ImageType::PointType + FindExtremaPointInAGivenLine(const ImageType * input, + int dimension, + bool inverse, + typename ImageType::PointType p, + typename ImageType::PixelType BG, + double distanceMax) + { + // Which direction ? Increasing or decreasing. + int d=1; + if (inverse) d=-1; + + // Transform to pixel index + typename ImageType::IndexType index; + input->TransformPhysicalPointToIndex(p, index); + + // Loop while inside the mask; + while (input->GetPixel(index) != BG) { + index[dimension] += d; } - siter.NextSlice(); + + // Transform back to Physical Units + typename ImageType::PointType result; + input->TransformIndexToPhysicalPoint(index, result); + + // Check that is is not too far away + double distance = p.EuclideanDistanceTo(result); + if (distance > distanceMax) { + result = p; // Get back to initial value + } + + return result; + } + //-------------------------------------------------------------------- + + + //-------------------------------------------------------------------- + template + bool + IsOnTheSameLineSide(PointType C, PointType A, PointType B, PointType like) + { + // Look at the position of point 'like' according to the AB line + double s = (B[0] - A[0]) * (like[1] - A[1]) - (B[1] - A[1]) * (like[0] - A[0]); + bool negative = s<0; + + // Look the C position + s = (B[0] - A[0]) * (C[1] - A[1]) - (B[1] - A[1]) * (C[0] - A[0]); + + if (negative && (s<=0)) return true; + if (!negative && (s>=0)) return true; + return false; } -} -//-------------------------------------------------------------------- -} + //-------------------------------------------------------------------- + + +} // end of namespace + -- 2.45.1