diff --git a/Modules/Filtering/ParabolicMorphology/CMakeLists.txt b/Modules/Filtering/ParabolicMorphology/CMakeLists.txt new file mode 100644 index 00000000000..e55581e3de3 --- /dev/null +++ b/Modules/Filtering/ParabolicMorphology/CMakeLists.txt @@ -0,0 +1,3 @@ +project(ParabolicMorphology) + +itk_module_impl() diff --git a/Modules/Filtering/ParabolicMorphology/README.md b/Modules/Filtering/ParabolicMorphology/README.md new file mode 100644 index 00000000000..a19c2514fbd --- /dev/null +++ b/Modules/Filtering/ParabolicMorphology/README.md @@ -0,0 +1,29 @@ +# ParabolicMorphology + +In-tree ITK module providing fast mathematical morphology using +parabolic structuring functions. Implements distance transforms, +binary erosions, dilations, openings, closings by spheres, grayscale +operations, and sharpenings in linear time per pixel via +separable parabolic filtering. + +## Origin + +Ingested from the standalone remote module +[**InsightSoftwareConsortium/ITKParabolicMorphology**](https://github.com/InsightSoftwareConsortium/ITKParabolicMorphology) +on 2026-05-12, following the v4 ingestion guidelines defined in +InsightSoftwareConsortium/ITK#6204. The upstream repository will be +archived read-only after this PR merges; it remains reachable at the +URL above for historical reference (notably the `doc/`, `examples/`, +and `oldtests/` directories, which were intentionally left in the +upstream archive). + +## References + +- Beare R. *Morphology with parabolic structuring elements.* The + Insight Journal. January-June. 2008. + +- Beare R., Jackway P. *Parabolic Morphology in ITK.* The Insight + Journal. 2012. +- Felzenszwalb P., Huttenlocher D. *Distance Transforms of Sampled + Functions.* Theory of Computing 8(19), 415-428, 2012. + diff --git a/Modules/Filtering/ParabolicMorphology/include/itkBinaryCloseParabolicImageFilter.h b/Modules/Filtering/ParabolicMorphology/include/itkBinaryCloseParabolicImageFilter.h new file mode 100644 index 00000000000..38805016b84 --- /dev/null +++ b/Modules/Filtering/ParabolicMorphology/include/itkBinaryCloseParabolicImageFilter.h @@ -0,0 +1,186 @@ +/*========================================================================= + * + * Copyright NumFOCUS + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * https://www.apache.org/licenses/LICENSE-2.0.txt + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + *=========================================================================*/ +#ifndef itkBinaryCloseParabolicImageFilter_h +#define itkBinaryCloseParabolicImageFilter_h + +#include "itkParabolicErodeImageFilter.h" +#include "itkParabolicDilateImageFilter.h" +#include "itkGreaterEqualValImageFilter.h" +#include "itkBinaryThresholdImageFilter.h" + +namespace itk +{ +/** + * \class BinaryCloseParabolicImageFilter + * \brief Class for binary morphological opening operation. + * + * This class uses the parabolic morphology operations to do very + * efficient erosions by circles/spheres. The operations are efficient + * because the underlying parabolic operations are separable and + * the operations are implicitly short circuited in comparison to a + * full distance transform approach. + * + * The basic idea is that a binary erosion or dilation by a circle/sphere + * can be carried out by thresholding a distance transform. By using + * the parabolic filters we can avoid computing the entire distance + * transform and instead only compute the subset we are interested + * in. + * + * Note that the circles and spheres may not be quite what you + * expect, because this class doesn't explicitly use Bresenham circles + * as most of the others do. A voxel's centre needs to be less than or + * equal to the circle radius, rather than any part of the voxel + * inside the circle. + * + * Also note that the inputs must be 0/1 not 0/max for pixel type. + * + * This filter was developed as a result of discussions with + * M.Starring on the ITK mailing list. + * + * Core methods described in the InsightJournal article: + * "Morphology with parabolic structuring elements" + * + * https://hdl.handle.net/1926/1370 + * + * \sa itkParabolicErodeImageFilter + * + * \ingroup ParabolicMorphology + * + * \author Richard Beare, Department of Medicine, Monash University, + * Australia. + **/ + +template +class ITK_TEMPLATE_EXPORT BinaryCloseParabolicImageFilter : public ImageToImageFilter + +{ +public: + ITK_DISALLOW_COPY_AND_MOVE(BinaryCloseParabolicImageFilter); + + /** Standard class type alias. */ + using Self = BinaryCloseParabolicImageFilter; + using Superclass = ImageToImageFilter; + using Pointer = SmartPointer; + using ConstPointer = SmartPointer; + + /** Method for creation through the object factory. */ + itkNewMacro(Self); + + /** Runtime information support. */ + itkTypeMacro(BinaryCloseParabolicImageFilter, ImageToImageFilter); + + /** Pixel Type of the input image */ + using InputImageType = TInputImage; + using OutputImageType = TOutputImage; + using PixelType = typename TInputImage::PixelType; + using RealType = typename NumericTraits::RealType; + using ScalarRealType = typename NumericTraits::ScalarRealType; + using OutputPixelType = typename TOutputImage::PixelType; + + /** Smart pointer type alias support. */ + using InputImagePointer = typename TInputImage::Pointer; + using InputImageConstPointer = typename TInputImage::ConstPointer; + + using InternalRealType = typename NumericTraits::FloatType; + // perhaps a bit dodgy, change to int if you want to do enormous + // binary operations + using InternalIntType = short; + + /** Image dimension. */ + static constexpr unsigned int ImageDimension = TInputImage::ImageDimension; + + using RadiusType = typename itk::FixedArray; + + void + SetRadius(ScalarRealType radius); + + itkSetMacro(Radius, RadiusType); + itkGetConstReferenceMacro(Radius, RadiusType); + + void + SetUseImageSpacing(bool g) + { + m_RectErode->SetUseImageSpacing(g); + m_RectDilate->SetUseImageSpacing(g); + m_CircErode->SetUseImageSpacing(g); + m_CircDilate->SetUseImageSpacing(g); + } + + /** + * Set/Get whether the erosion is circular/rectangular - + * default is true (circular) + */ + itkSetMacro(Circular, bool); + itkGetConstReferenceMacro(Circular, bool); + itkBooleanMacro(Circular); + + /** A safe border is added to input image to avoid borders effects + * and remove it once the closing is done */ + itkSetMacro(SafeBorder, bool); + itkGetConstReferenceMacro(SafeBorder, bool); + itkBooleanMacro(SafeBorder); + + /** Image related type alias. */ + + /* add in the traits here */ +protected: + void + GenerateData() override; + + BinaryCloseParabolicImageFilter(); + ~BinaryCloseParabolicImageFilter() override = default; + void + PrintSelf(std::ostream & os, Indent indent) const override; + + using InternalRealImageType = typename itk::Image; + using InternalIntImageType = typename itk::Image; + using CircErodeType = typename itk::ParabolicErodeImageFilter; + using RectErodeType = typename itk::ParabolicErodeImageFilter; + using CircDilateType = typename itk::ParabolicDilateImageFilter; + using RectDilateType = typename itk::ParabolicDilateImageFilter; + + using CCastTypeA = typename itk::GreaterEqualValImageFilter; + using CCastTypeB = typename itk::BinaryThresholdImageFilter; + + using RCastTypeA = typename itk::GreaterEqualValImageFilter; + using RCastTypeB = typename itk::BinaryThresholdImageFilter; + +private: + RadiusType m_Radius; + bool m_Circular; + bool m_SafeBorder; + + typename CircErodeType::Pointer m_CircErode; + typename CircDilateType::Pointer m_CircDilate; + + typename CCastTypeA::Pointer m_CircCastA; + typename CCastTypeB::Pointer m_CircCastB; + + typename RectErodeType::Pointer m_RectErode; + typename RectDilateType::Pointer m_RectDilate; + + typename RCastTypeA::Pointer m_RectCastA; + typename RCastTypeB::Pointer m_RectCastB; +}; +} // end namespace itk + +#ifndef ITK_MANUAL_INSTANTIATION +# include "itkBinaryCloseParabolicImageFilter.hxx" +#endif + +#endif //__itkBinaryCloseParabolicImageFilter_h diff --git a/Modules/Filtering/ParabolicMorphology/include/itkBinaryCloseParabolicImageFilter.hxx b/Modules/Filtering/ParabolicMorphology/include/itkBinaryCloseParabolicImageFilter.hxx new file mode 100644 index 00000000000..db0c621d5bc --- /dev/null +++ b/Modules/Filtering/ParabolicMorphology/include/itkBinaryCloseParabolicImageFilter.hxx @@ -0,0 +1,235 @@ +/*========================================================================= + * + * Copyright NumFOCUS + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * https://www.apache.org/licenses/LICENSE-2.0.txt + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + *=========================================================================*/ +#ifndef itkBinaryCloseParabolicImageFilter_hxx +#define itkBinaryCloseParabolicImageFilter_hxx + +#include "itkProgressAccumulator.h" +#include "itkParabolicErodeImageFilter.h" +#include "itkCropImageFilter.h" +#include "itkConstantPadImageFilter.h" +#include "itkMath.h" + +namespace itk +{ +template +BinaryCloseParabolicImageFilter::BinaryCloseParabolicImageFilter() +{ + this->SetNumberOfRequiredOutputs(1); + this->SetNumberOfRequiredInputs(1); + this->m_CircErode = CircErodeType::New(); + this->m_CircDilate = CircDilateType::New(); + this->m_CircCastA = CCastTypeA::New(); + this->m_CircCastB = CCastTypeB::New(); + + this->m_RectErode = RectErodeType::New(); + this->m_RectDilate = RectDilateType::New(); + this->m_RectCastA = RCastTypeA::New(); + this->m_RectCastB = RCastTypeB::New(); + this->m_Circular = true; + // Need to call this after filters are created + this->SetUseImageSpacing(false); + this->SetSafeBorder(true); +} + +template +void +BinaryCloseParabolicImageFilter::SetRadius(ScalarRealType radius) +{ + RadiusType s; + + s.Fill(radius); + this->SetRadius(s); +} + +template +void +BinaryCloseParabolicImageFilter::GenerateData() +{ + // Allocate the output + this->AllocateOutputs(); + // set up the scaling before we pass control over to superclass + typename TInputImage::SizeType Pad; + // ScalarRealType margin = 0.0; + + // ScalarRealType mxRad = (ScalarRealType)(*std::max_element(m_Radius.Begin(), + // m_Radius.End())); + // this needs to be examined more closely + // margin = 1.0/(pow(mxRad, TInputImage::ImageDimension) * 10); + // margin = std::min(margin, 0.00001); + // std::cout << "Margin = " << margin << std::endl; + + if (this->m_RectErode->GetUseImageSpacing()) + { + // radius is in mm + RadiusType R; + for (unsigned P = 0; P < InputImageType::ImageDimension; P++) + { + typename TInputImage::SpacingValueType tsp = this->GetInput()->GetSpacing()[P]; + R[P] = 0.5 * (m_Radius[P] * m_Radius[P]) + tsp * tsp; + Pad[P] = (typename TInputImage::SizeType::SizeValueType)(itk::Math::rnd_halfinttoeven(m_Radius[P] / tsp + 1) + 1); + } + m_RectErode->SetScale(R); + m_CircErode->SetScale(R); + m_RectDilate->SetScale(R); + m_CircDilate->SetScale(R); + } + else + { + // radius is in pixels + RadiusType R; + // this gives us a little bit of a margin + for (unsigned P = 0; P < InputImageType::ImageDimension; P++) + { + R[P] = (0.5 * m_Radius[P] * m_Radius[P] + 1); + Pad[P] = (typename TInputImage::SizeType::SizeValueType)(m_Radius[P] + 1); + } + // std::cout << "no image spacing " << m_Radius << R << std::endl; + m_RectErode->SetScale(R); + m_CircErode->SetScale(R); + m_RectDilate->SetScale(R); + m_CircDilate->SetScale(R); + } + + // std::cout << "Padding " << Pad << std::endl; + + if (m_Circular) + { + ProgressAccumulator::Pointer progress = ProgressAccumulator::New(); + progress->SetMiniPipelineFilter(this); + InputImageConstPointer inputImage; + inputImage = this->GetInput(); + + progress->RegisterInternalFilter(m_CircErode, 0.4f); + progress->RegisterInternalFilter(m_CircCastA, 0.1f); + progress->RegisterInternalFilter(m_CircDilate, 0.4f); + progress->RegisterInternalFilter(m_CircCastB, 0.1f); + + m_CircCastB->SetInput(m_CircDilate->GetOutput()); + // m_CircCastB->SetUpperThreshold(margin); + m_CircCastB->SetUpperThreshold(0); + m_CircCastB->SetOutsideValue(1); + m_CircCastB->SetInsideValue(0); + + m_CircErode->SetInput(m_CircCastB->GetOutput()); + m_CircCastA->SetInput(m_CircErode->GetOutput()); + // m_CircCastA->SetVal(1.0-margin); + m_CircCastA->SetVal(1.0); + + if (m_SafeBorder) + { + using PadType = typename itk::ConstantPadImageFilter; + typename PadType::Pointer pad = PadType::New(); + pad->SetPadLowerBound(Pad); + pad->SetPadUpperBound(Pad); + pad->SetConstant(0); + pad->SetInput(inputImage); + + m_CircDilate->SetInput(pad->GetOutput()); + + // writeIm(m_CircCastB->GetOutput(), "dil.nii.gz"); + // writeIm(m_CircCastA->GetOutput(), "ero.nii.gz"); + // m_CircCastA->UpdateOutputInformation(); + using CropType = typename itk::CropImageFilter; + typename CropType::Pointer crop = CropType::New(); + crop->SetInput(m_CircCastA->GetOutput()); + crop->SetUpperBoundaryCropSize(Pad); + crop->SetLowerBoundaryCropSize(Pad); + + crop->GraftOutput(this->GetOutput()); + crop->Update(); + + this->GraftOutput(crop->GetOutput()); + } + else + { + m_CircDilate->SetInput(inputImage); + m_CircCastA->GraftOutput(this->GetOutput()); + m_CircCastA->Update(); + + this->GraftOutput(m_CircCastA->GetOutput()); + } + } + else + { + ProgressAccumulator::Pointer progress = ProgressAccumulator::New(); + progress->SetMiniPipelineFilter(this); + InputImageConstPointer inputImage; + inputImage = this->GetInput(); + + progress->RegisterInternalFilter(m_RectErode, 0.4f); + progress->RegisterInternalFilter(m_RectCastA, 0.1f); + progress->RegisterInternalFilter(m_RectDilate, 0.4f); + progress->RegisterInternalFilter(m_RectCastB, 0.1f); + + m_RectCastB->SetInput(m_RectDilate->GetOutput()); + m_RectCastB->SetUpperThreshold(0); + m_RectCastB->SetOutsideValue(1); + m_RectCastB->SetInsideValue(0); + + m_RectErode->SetInput(m_RectCastB->GetOutput()); + m_RectCastA->SetInput(m_RectErode->GetOutput()); + m_RectCastA->SetVal(1); + + if (m_SafeBorder) + { + using PadType = typename itk::ConstantPadImageFilter; + typename PadType::Pointer pad = PadType::New(); + pad->SetPadLowerBound(Pad); + pad->SetPadUpperBound(Pad); + pad->SetConstant(0); + pad->SetInput(inputImage); + + m_RectDilate->SetInput(pad->GetOutput()); + + using CropType = typename itk::CropImageFilter; + typename CropType::Pointer crop = CropType::New(); + crop->SetInput(m_RectCastA->GetOutput()); + crop->SetUpperBoundaryCropSize(Pad); + crop->SetLowerBoundaryCropSize(Pad); + + crop->GraftOutput(this->GetOutput()); + crop->Update(); + this->GraftOutput(crop->GetOutput()); + } + else + { + m_RectDilate->SetInput(inputImage); + m_RectCastA->GraftOutput(this->GetOutput()); + m_RectCastA->Update(); + this->GraftOutput(m_RectCastA->GetOutput()); + } + } +} + +template +void +BinaryCloseParabolicImageFilter::PrintSelf(std::ostream & os, Indent indent) const +{ + Superclass::PrintSelf(os, indent); + if (this->m_CircErode->GetUseImageSpacing()) + { + os << "Radius in world units: " << this->GetRadius() << std::endl; + } + else + { + os << "Radius in voxels: " << this->GetRadius() << std::endl; + } + os << "Safe border: " << this->GetSafeBorder() << std::endl; +} +} // namespace itk +#endif diff --git a/Modules/Filtering/ParabolicMorphology/include/itkBinaryDilateParabolicImageFilter.h b/Modules/Filtering/ParabolicMorphology/include/itkBinaryDilateParabolicImageFilter.h new file mode 100644 index 00000000000..3710b4c523a --- /dev/null +++ b/Modules/Filtering/ParabolicMorphology/include/itkBinaryDilateParabolicImageFilter.h @@ -0,0 +1,165 @@ +/*========================================================================= + * + * Copyright NumFOCUS + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * https://www.apache.org/licenses/LICENSE-2.0.txt + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + *=========================================================================*/ +#ifndef itkBinaryDilateParabolicImageFilter_h +#define itkBinaryDilateParabolicImageFilter_h + +#include "itkParabolicDilateImageFilter.h" +#include "itkBinaryThresholdImageFilter.h" + +namespace itk +{ +/** + * \class BinaryDilateParabolicImageFilter + * \brief Class for binary morphological dilation operation. + * + * This class uses the parabolic morphology operations to do very + * efficient erosions by circles/spheres. The operations are efficient + * because the underlying parabolic operations are separable and + * the operations are implicitly short circuited in comparison to a + * full distance transform approach. + * + * The basic idea is that a binary erosion or dilation by a circle/sphere + * can be carried out by thresholding a distance transform. By using + * the parabolic filters we can avoid computing the entire distance + * transform and instead only compute the subset we are interested + * in. + * + * Note that the circles and spheres may not be quite what you + * expect, because this class doesn't explicitly use Bresenham circles + * as most of the others do. A voxel's centre needs to be less than or + * equal to the circle radius, rather than any part of the voxel + * inside the circle. + * + * This filter was developed as a result of discussions with + * M.Starring on the ITK mailing list. + * + * Also note that the inputs must be 0/1 not 0/max for pixel type. + * + * Core methods described in the InsightJournal article: + * "Morphology with parabolic structuring elements" + * + * https://hdl.handle.net/1926/1370 + * + * \sa itkParabolicDilateImageFilter + * + * \ingroup ParabolicMorphology + * + * \author Richard Beare, Department of Medicine, Monash University, + * Australia. + **/ + +template +class ITK_TEMPLATE_EXPORT BinaryDilateParabolicImageFilter : public ImageToImageFilter + +{ +public: + ITK_DISALLOW_COPY_AND_MOVE(BinaryDilateParabolicImageFilter); + + /** Standard class type alias. */ + using Self = BinaryDilateParabolicImageFilter; + using Superclass = ImageToImageFilter; + using Pointer = SmartPointer; + using ConstPointer = SmartPointer; + + /** Method for creation through the object factory. */ + itkNewMacro(Self); + + /** Runtime information support. */ + itkTypeMacro(BinaryDilateParabolicImageFilter, ImageToImageFilter); + + /** Pixel Type of the input image */ + using InputImageType = TInputImage; + using OutputImageType = TOutputImage; + using PixelType = typename TInputImage::PixelType; + using RealType = typename NumericTraits::RealType; + using ScalarRealType = typename NumericTraits::ScalarRealType; + using OutputPixelType = typename TOutputImage::PixelType; + + /** Smart pointer type alias support. */ + using InputImagePointer = typename TInputImage::Pointer; + using InputImageConstPointer = typename TInputImage::ConstPointer; + + using InternalRealType = typename NumericTraits::FloatType; + // perhaps a bit dodgy, change to int if you want to do enormous + // binary operations + using InternalIntType = short; + + /** Image dimension. */ + static constexpr unsigned int ImageDimension = TInputImage::ImageDimension; + + using RadiusType = typename itk::FixedArray; + + void + SetRadius(ScalarRealType radius); + + itkSetMacro(Radius, RadiusType); + itkGetConstReferenceMacro(Radius, RadiusType); + + void + Modified() const override; + + void + SetUseImageSpacing(bool g) + { + m_RectPara->SetUseImageSpacing(g); + m_CircPara->SetUseImageSpacing(g); + } + + /** + * Set/Get whether the erosion is circular/rectangular - + * default is true (circular) + */ + itkSetMacro(Circular, bool); + itkGetConstReferenceMacro(Circular, bool); + itkBooleanMacro(Circular); + /** Image related type alias. */ + + /* add in the traits here */ +protected: + void + GenerateData() override; + + BinaryDilateParabolicImageFilter(); + ~BinaryDilateParabolicImageFilter() override = default; + void + PrintSelf(std::ostream & os, Indent indent) const override; + + using InternalRealImageType = typename itk::Image; + using InternalIntImageType = typename itk::Image; + using CircParabolicType = typename itk::ParabolicDilateImageFilter; + using RectParabolicType = typename itk::ParabolicDilateImageFilter; + using CCastType = typename itk::BinaryThresholdImageFilter; + using RCastType = typename itk::BinaryThresholdImageFilter; + +private: + RadiusType m_Radius; + bool m_Circular; + + typename CircParabolicType::Pointer m_CircPara; + typename CCastType::Pointer m_CircCast; + + typename RectParabolicType::Pointer m_RectPara; + typename RCastType::Pointer m_RectCast; +}; +} // end namespace itk + +#ifndef ITK_MANUAL_INSTANTIATION +# include "itkBinaryDilateParabolicImageFilter.hxx" +#endif + +#endif //__itkBinaryDilateParabolicImageFilter_h diff --git a/Modules/Filtering/ParabolicMorphology/include/itkBinaryDilateParabolicImageFilter.hxx b/Modules/Filtering/ParabolicMorphology/include/itkBinaryDilateParabolicImageFilter.hxx new file mode 100644 index 00000000000..ba64a971e5a --- /dev/null +++ b/Modules/Filtering/ParabolicMorphology/include/itkBinaryDilateParabolicImageFilter.hxx @@ -0,0 +1,155 @@ +/*========================================================================= + * + * Copyright NumFOCUS + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * https://www.apache.org/licenses/LICENSE-2.0.txt + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + *=========================================================================*/ +#ifndef itkBinaryDilateParabolicImageFilter_hxx +#define itkBinaryDilateParabolicImageFilter_hxx + +#include "itkProgressAccumulator.h" +#include "itkParabolicDilateImageFilter.h" + +namespace itk +{ +template +BinaryDilateParabolicImageFilter::BinaryDilateParabolicImageFilter() +{ + this->SetNumberOfRequiredOutputs(1); + this->SetNumberOfRequiredInputs(1); + this->m_CircPara = CircParabolicType::New(); + this->m_CircCast = CCastType::New(); + + this->m_RectPara = RectParabolicType::New(); + this->m_RectCast = RCastType::New(); + this->m_Circular = true; + // Need to call this after filters are created + this->SetUseImageSpacing(false); +} + +template +void +BinaryDilateParabolicImageFilter::SetRadius(ScalarRealType radius) +{ + RadiusType s; + + s.Fill(radius); + this->SetRadius(s); +} + +template +void +BinaryDilateParabolicImageFilter::GenerateData() +{ + // Allocate the output + this->AllocateOutputs(); + // set up the scaling before we pass control over to superclass + if (this->m_RectPara->GetUseImageSpacing()) + { + // radius is in mm + RadiusType R; + for (unsigned P = 0; P < InputImageType::ImageDimension; P++) + { + R[P] = 0.5 * m_Radius[P] * m_Radius[P]; + // this->SetScale(0.5*m_Radius[P] * m_Radius[P]); + } + m_RectPara->SetScale(R); + m_CircPara->SetScale(R); + } + else + { + // radius is in pixels + RadiusType R; + // this gives us a little bit of a margin + for (unsigned P = 0; P < InputImageType::ImageDimension; P++) + { + R[P] = (0.5 * m_Radius[P] * m_Radius[P] + 1); + } + // std::cout << "no image spacing " << m_Radius << R << std::endl; + m_RectPara->SetScale(R); + m_CircPara->SetScale(R); + } + + if (m_Circular) + { + ProgressAccumulator::Pointer progress = ProgressAccumulator::New(); + progress->SetMiniPipelineFilter(this); + InputImageConstPointer inputImage; + inputImage = this->GetInput(); + + progress->RegisterInternalFilter(m_CircPara, 0.8f); + progress->RegisterInternalFilter(m_CircCast, 0.2f); + + m_CircPara->SetInput(inputImage); + m_CircCast->SetInput(m_CircPara->GetOutput()); + // m_CircCast->SetInsideValue(0); + // m_CircCast->SetOutsideValue(1); + // setting the correct threshold value is a little tricky - needs would + // to produce a result matching a bresenham circle, but these + // circles are such that the voxel centres need to be less than radius + m_CircCast->SetUpperThreshold(0); + m_CircCast->SetOutsideValue(1); + m_CircCast->SetInsideValue(0); + m_CircCast->GraftOutput(this->GetOutput()); + m_CircCast->Update(); + this->GraftOutput(m_CircCast->GetOutput()); + } + else + { + ProgressAccumulator::Pointer progress = ProgressAccumulator::New(); + progress->SetMiniPipelineFilter(this); + InputImageConstPointer inputImage; + inputImage = this->GetInput(); + + progress->RegisterInternalFilter(m_RectPara, 0.8f); + progress->RegisterInternalFilter(m_RectCast, 0.2f); + + m_RectPara->SetInput(inputImage); + m_RectCast->SetInput(m_RectPara->GetOutput()); + m_RectCast->SetUpperThreshold(0); + m_RectCast->SetOutsideValue(1); + m_RectCast->SetInsideValue(0); + m_RectCast->GraftOutput(this->GetOutput()); + m_RectCast->Update(); + this->GraftOutput(m_RectCast->GetOutput()); + } +} + +template +void +BinaryDilateParabolicImageFilter::Modified() const +{ + Superclass::Modified(); + m_CircPara->Modified(); + m_CircCast->Modified(); + m_RectPara->Modified(); + m_RectCast->Modified(); +} + +template +void +BinaryDilateParabolicImageFilter::PrintSelf(std::ostream & os, Indent indent) const +{ + Superclass::PrintSelf(os, indent); + if (this->m_CircPara->GetUseImageSpacing()) + { + os << "Radius in world units: " << this->GetRadius() << std::endl; + } + else + { + os << "Radius in voxels: " << this->GetRadius() << std::endl; + } +} +} // namespace itk +#endif diff --git a/Modules/Filtering/ParabolicMorphology/include/itkBinaryErodeParabolicImageFilter.h b/Modules/Filtering/ParabolicMorphology/include/itkBinaryErodeParabolicImageFilter.h new file mode 100644 index 00000000000..fed4ed3b65e --- /dev/null +++ b/Modules/Filtering/ParabolicMorphology/include/itkBinaryErodeParabolicImageFilter.h @@ -0,0 +1,166 @@ +/*========================================================================= + * + * Copyright NumFOCUS + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * https://www.apache.org/licenses/LICENSE-2.0.txt + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + *=========================================================================*/ +#ifndef itkBinaryErodeParabolicImageFilter_h +#define itkBinaryErodeParabolicImageFilter_h + +#include "itkParabolicErodeImageFilter.h" +#include "itkGreaterEqualValImageFilter.h" + +namespace itk +{ +/** + * \class BinaryErodeParabolicImageFilter + * \brief Class for binary morphological erosion operation. + * + * This class uses the parabolic morphology operations to do very + * efficient erosions by circles/spheres. The operations are efficient + * because the underlying parabolic operations are separable and + * the operations are implicitly short circuited in comparison to a + * full distance transform approach. + * + * The basic idea is that a binary erosion or dilation by a circle/sphere + * can be carried out by thresholding a distance transform. By using + * the parabolic filters we can avoid computing the entire distance + * transform and instead only compute the subset we are interested + * in. + * + * Note that the circles and spheres may not be quite what you + * expect, because this class doesn't explicitly use Bresenham circles + * as most of the others do. A voxel's centre needs to be less than or + * equal to the circle radius, rather than any part of the voxel + * inside the circle. + * + * Also note that the inputs must be 0/1 not 0/max for pixel type. + * + * This filter was developed as a result of discussions with + * M.Starring on the ITK mailing list. + * + * Core methods described in the InsightJournal article: + * "Morphology with parabolic structuring elements" + * + * https://hdl.handle.net/1926/1370 + * + * \sa itkParabolicErodeImageFilter + * + * + * \ingroup ParabolicMorphology + * + * \author Richard Beare, Department of Medicine, Monash University, + * Australia. + **/ + +template +class ITK_TEMPLATE_EXPORT BinaryErodeParabolicImageFilter : public ImageToImageFilter + +{ +public: + ITK_DISALLOW_COPY_AND_MOVE(BinaryErodeParabolicImageFilter); + + /** Standard class type alias. */ + using Self = BinaryErodeParabolicImageFilter; + using Superclass = ImageToImageFilter; + using Pointer = SmartPointer; + using ConstPointer = SmartPointer; + + /** Method for creation through the object factory. */ + itkNewMacro(Self); + + /** Runtime information support. */ + itkTypeMacro(BinaryErodeParabolicImageFilter, ImageToImageFilter); + + /** Pixel Type of the input image */ + using InputImageType = TInputImage; + using OutputImageType = TOutputImage; + using PixelType = typename TInputImage::PixelType; + using RealType = typename NumericTraits::RealType; + using ScalarRealType = typename NumericTraits::ScalarRealType; + using OutputPixelType = typename TOutputImage::PixelType; + + /** Smart pointer type alias support. */ + using InputImagePointer = typename TInputImage::Pointer; + using InputImageConstPointer = typename TInputImage::ConstPointer; + + using InternalRealType = typename NumericTraits::FloatType; + // perhaps a bit dodgy, change to int if you want to do enormous + // binary operations + using InternalIntType = short; + + /** Image dimension. */ + static constexpr unsigned int ImageDimension = TInputImage::ImageDimension; + + using RadiusType = typename itk::FixedArray; + + void + SetRadius(ScalarRealType radius); + + itkSetMacro(Radius, RadiusType); + itkGetConstReferenceMacro(Radius, RadiusType); + + void + SetUseImageSpacing(bool g) + { + m_RectPara->SetUseImageSpacing(g); + m_CircPara->SetUseImageSpacing(g); + } + + /** + * Set/Get whether the erosion is circular/rectangular - + * default is true (circular) + */ + itkSetMacro(Circular, bool); + itkGetConstReferenceMacro(Circular, bool); + itkBooleanMacro(Circular); + /** Image related type alias. */ + + /* add in the traits here */ + void + Modified() const override; + +protected: + void + GenerateData() override; + + BinaryErodeParabolicImageFilter(); + ~BinaryErodeParabolicImageFilter() override = default; + void + PrintSelf(std::ostream & os, Indent indent) const override; + + using InternalRealImageType = typename itk::Image; + using InternalIntImageType = typename itk::Image; + using CircParabolicType = typename itk::ParabolicErodeImageFilter; + using RectParabolicType = typename itk::ParabolicErodeImageFilter; + using CCastType = typename itk::GreaterEqualValImageFilter; + using RCastType = typename itk::GreaterEqualValImageFilter; + +private: + RadiusType m_Radius; + bool m_Circular; + + typename CircParabolicType::Pointer m_CircPara; + typename CCastType::Pointer m_CircCast; + + typename RectParabolicType::Pointer m_RectPara; + typename RCastType::Pointer m_RectCast; +}; +} // end namespace itk + +#ifndef ITK_MANUAL_INSTANTIATION +# include "itkBinaryErodeParabolicImageFilter.hxx" +#endif + +#endif //__itkBinaryErodeParabolicImageFilter_h diff --git a/Modules/Filtering/ParabolicMorphology/include/itkBinaryErodeParabolicImageFilter.hxx b/Modules/Filtering/ParabolicMorphology/include/itkBinaryErodeParabolicImageFilter.hxx new file mode 100644 index 00000000000..6519659f858 --- /dev/null +++ b/Modules/Filtering/ParabolicMorphology/include/itkBinaryErodeParabolicImageFilter.hxx @@ -0,0 +1,155 @@ +/*========================================================================= + * + * Copyright NumFOCUS + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * https://www.apache.org/licenses/LICENSE-2.0.txt + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + *=========================================================================*/ +#ifndef itkBinaryErodeParabolicImageFilter_hxx +#define itkBinaryErodeParabolicImageFilter_hxx + +#include "itkProgressAccumulator.h" +#include "itkParabolicErodeImageFilter.h" + +namespace itk +{ +template +BinaryErodeParabolicImageFilter::BinaryErodeParabolicImageFilter() +{ + this->SetNumberOfRequiredOutputs(1); + this->SetNumberOfRequiredInputs(1); + this->m_CircPara = CircParabolicType::New(); + this->m_CircCast = CCastType::New(); + + this->m_RectPara = RectParabolicType::New(); + this->m_RectCast = RCastType::New(); + this->m_Circular = true; + // Need to call this after filters are created + this->SetUseImageSpacing(false); +} + +template +void +BinaryErodeParabolicImageFilter::SetRadius(ScalarRealType radius) +{ + RadiusType s; + + s.Fill(radius); + this->SetRadius(s); +} + +template +void +BinaryErodeParabolicImageFilter::GenerateData() +{ + // Allocate the output + this->AllocateOutputs(); + // set up the scaling before we pass control over to superclass + if (this->m_RectPara->GetUseImageSpacing()) + { + // radius is in mm + RadiusType R; + for (unsigned P = 0; P < InputImageType::ImageDimension; P++) + { + R[P] = 0.5 * m_Radius[P] * m_Radius[P]; + // this->SetScale(0.5*m_Radius[P] * m_Radius[P]); + } + m_RectPara->SetScale(R); + m_CircPara->SetScale(R); + } + else + { + // radius is in pixels + RadiusType R; + // this gives us a little bit of a margin + for (unsigned P = 0; P < InputImageType::ImageDimension; P++) + { + R[P] = (0.5 * m_Radius[P] * m_Radius[P] + 1); + } + // std::cout << "no image spacing " << m_Radius << R << std::endl; + m_RectPara->SetScale(R); + m_CircPara->SetScale(R); + } + + if (m_Circular) + { + ProgressAccumulator::Pointer progress = ProgressAccumulator::New(); + progress->SetMiniPipelineFilter(this); + InputImageConstPointer inputImage; + inputImage = this->GetInput(); + + progress->RegisterInternalFilter(m_CircPara, 0.8f); + progress->RegisterInternalFilter(m_CircCast, 0.2f); + + m_CircPara->SetInput(inputImage); + m_CircCast->SetInput(m_CircPara->GetOutput()); + m_CircCast->SetVal(1.0); + // m_CircCast->SetInsideValue(0); + // m_CircCast->SetOutsideValue(1); + // setting the correct threshold value is a little tricky - needs would + // to produce a result matching a bresenham circle, but these + // circles are such that the voxel centres need to be less than radius + // m_CircCast->SetUpperThreshold(1 - + // itk::NumericTraits::min()); + // m_CircCast->SetUpperThreshold(0.99); + + m_CircCast->GraftOutput(this->GetOutput()); + m_CircCast->Update(); + this->GraftOutput(m_CircCast->GetOutput()); + } + else + { + ProgressAccumulator::Pointer progress = ProgressAccumulator::New(); + progress->SetMiniPipelineFilter(this); + InputImageConstPointer inputImage; + inputImage = this->GetInput(); + + progress->RegisterInternalFilter(m_RectPara, 0.8f); + progress->RegisterInternalFilter(m_RectCast, 0.2f); + + m_RectPara->SetInput(inputImage); + m_RectCast->SetInput(m_RectPara->GetOutput()); + m_RectCast->SetVal(1); + m_RectCast->GraftOutput(this->GetOutput()); + m_RectCast->Update(); + this->GraftOutput(m_RectCast->GetOutput()); + } +} + +template +void +BinaryErodeParabolicImageFilter::PrintSelf(std::ostream & os, Indent indent) const +{ + Superclass::PrintSelf(os, indent); + if (this->m_CircPara->GetUseImageSpacing()) + { + os << "Radius in world units: " << this->GetRadius() << std::endl; + } + else + { + os << "Radius in voxels: " << this->GetRadius() << std::endl; + } +} + +template +void +BinaryErodeParabolicImageFilter::Modified() const +{ + Superclass::Modified(); + m_CircPara->Modified(); + m_CircCast->Modified(); + m_RectPara->Modified(); + m_RectCast->Modified(); +} +} // namespace itk +#endif diff --git a/Modules/Filtering/ParabolicMorphology/include/itkBinaryOpenParabolicImageFilter.h b/Modules/Filtering/ParabolicMorphology/include/itkBinaryOpenParabolicImageFilter.h new file mode 100644 index 00000000000..8583bad62bc --- /dev/null +++ b/Modules/Filtering/ParabolicMorphology/include/itkBinaryOpenParabolicImageFilter.h @@ -0,0 +1,186 @@ +/*========================================================================= + * + * Copyright NumFOCUS + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * https://www.apache.org/licenses/LICENSE-2.0.txt + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + *=========================================================================*/ +#ifndef itkBinaryOpenParabolicImageFilter_h +#define itkBinaryOpenParabolicImageFilter_h + +#include "itkParabolicErodeImageFilter.h" +#include "itkParabolicDilateImageFilter.h" +#include "itkGreaterEqualValImageFilter.h" +#include "itkBinaryThresholdImageFilter.h" + +namespace itk +{ +/** + * \class BinaryOpenParabolicImageFilter + * \brief Class for binary morphological opening operation. + * + * This class uses the parabolic morphology operations to do very + * efficient erosions by circles/spheres. The operations are efficient + * because the underlying parabolic operations are separable and + * the operations are implicitly short circuited in comparison to a + * full distance transform approach. + * + * The basic idea is that a binary erosion or dilation by a circle/sphere + * can be carried out by thresholding a distance transform. By using + * the parabolic filters we can avoid computing the entire distance + * transform and instead only compute the subset we are interested + * in. + * + * Note that the circles and spheres may not be quite what you + * expect, because this class doesn't explicitly use Bresenham circles + * as most of the others do. A voxel's centre needs to be less than or + * equal to the circle radius, rather than any part of the voxel + * inside the circle. + * + * This filter was developed as a result of discussions with + * M.Starring on the ITK mailing list. + * + * Also note that the inputs must be 0/1 not 0/max for pixel type. + * + * Core methods described in the InsightJournal article: + * "Morphology with parabolic structuring elements" + * + * https://hdl.handle.net/1926/1370 + * \sa itkParabolicErodeImageFilter + * + * + * \ingroup ParabolicMorphology + * + * \author Richard Beare, Department of Medicine, Monash University, + * Australia. + **/ + +template +class ITK_TEMPLATE_EXPORT BinaryOpenParabolicImageFilter : public ImageToImageFilter + +{ +public: + ITK_DISALLOW_COPY_AND_MOVE(BinaryOpenParabolicImageFilter); + + /** Standard class type alias. */ + using Self = BinaryOpenParabolicImageFilter; + using Superclass = ImageToImageFilter; + using Pointer = SmartPointer; + using ConstPointer = SmartPointer; + + /** Method for creation through the object factory. */ + itkNewMacro(Self); + + /** Runtime information support. */ + itkTypeMacro(BinaryOpenParabolicImageFilter, ImageToImageFilter); + + /** Pixel Type of the input image */ + using InputImageType = TInputImage; + using OutputImageType = TOutputImage; + using PixelType = typename TInputImage::PixelType; + using RealType = typename NumericTraits::RealType; + using ScalarRealType = typename NumericTraits::ScalarRealType; + using OutputPixelType = typename TOutputImage::PixelType; + + /** Smart pointer type alias support. */ + using InputImagePointer = typename TInputImage::Pointer; + using InputImageConstPointer = typename TInputImage::ConstPointer; + + using InternalRealType = typename NumericTraits::FloatType; + // perhaps a bit dodgy, change to int if you want to do enormous + // binary operations + using InternalIntType = short; + + /** Image dimension. */ + static constexpr unsigned int ImageDimension = TInputImage::ImageDimension; + + using RadiusType = typename itk::FixedArray; + + void + SetRadius(ScalarRealType radius); + + itkSetMacro(Radius, RadiusType); + itkGetConstReferenceMacro(Radius, RadiusType); + + void + SetUseImageSpacing(bool g) + { + m_RectErode->SetUseImageSpacing(g); + m_RectDilate->SetUseImageSpacing(g); + m_CircErode->SetUseImageSpacing(g); + m_CircDilate->SetUseImageSpacing(g); + } + + /** + * Set/Get whether the erosion is circular/rectangular - + * default is true (circular) + */ + itkSetMacro(Circular, bool); + itkGetConstReferenceMacro(Circular, bool); + itkBooleanMacro(Circular); + + /** A safe border is added to input image to avoid borders effects + * and remove it once the closing is done */ + itkSetMacro(SafeBorder, bool); + itkGetConstReferenceMacro(SafeBorder, bool); + itkBooleanMacro(SafeBorder); + + /** Image related type alias. */ + + /* add in the traits here */ +protected: + void + GenerateData() override; + + BinaryOpenParabolicImageFilter(); + ~BinaryOpenParabolicImageFilter() override = default; + void + PrintSelf(std::ostream & os, Indent indent) const override; + + using InternalRealImageType = typename itk::Image; + using InternalIntImageType = typename itk::Image; + using CircErodeType = typename itk::ParabolicErodeImageFilter; + using RectErodeType = typename itk::ParabolicErodeImageFilter; + using CircDilateType = typename itk::ParabolicDilateImageFilter; + using RectDilateType = typename itk::ParabolicDilateImageFilter; + + using CCastTypeA = typename itk::GreaterEqualValImageFilter; + using CCastTypeB = typename itk::BinaryThresholdImageFilter; + + using RCastTypeA = typename itk::GreaterEqualValImageFilter; + using RCastTypeB = typename itk::BinaryThresholdImageFilter; + +private: + RadiusType m_Radius; + bool m_Circular; + bool m_SafeBorder; + + typename CircErodeType::Pointer m_CircErode; + typename CircDilateType::Pointer m_CircDilate; + + typename CCastTypeA::Pointer m_CircCastA; + typename CCastTypeB::Pointer m_CircCastB; + + typename RectErodeType::Pointer m_RectErode; + typename RectDilateType::Pointer m_RectDilate; + + typename RCastTypeA::Pointer m_RectCastA; + typename RCastTypeB::Pointer m_RectCastB; +}; +} // end namespace itk + +#ifndef ITK_MANUAL_INSTANTIATION +# include "itkBinaryOpenParabolicImageFilter.hxx" +#endif + +#endif //__itkBinaryOpenParabolicImageFilter_h diff --git a/Modules/Filtering/ParabolicMorphology/include/itkBinaryOpenParabolicImageFilter.hxx b/Modules/Filtering/ParabolicMorphology/include/itkBinaryOpenParabolicImageFilter.hxx new file mode 100644 index 00000000000..634751a3565 --- /dev/null +++ b/Modules/Filtering/ParabolicMorphology/include/itkBinaryOpenParabolicImageFilter.hxx @@ -0,0 +1,255 @@ +/*========================================================================= + * + * Copyright NumFOCUS + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * https://www.apache.org/licenses/LICENSE-2.0.txt + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + *=========================================================================*/ +#ifndef itkBinaryOpenParabolicImageFilter_hxx +#define itkBinaryOpenParabolicImageFilter_hxx + +#include "itkProgressAccumulator.h" +#include "itkParabolicErodeImageFilter.h" +#include "itkCropImageFilter.h" +#include "itkConstantPadImageFilter.h" +#include "itkMath.h" + +namespace itk +{ +template +BinaryOpenParabolicImageFilter::BinaryOpenParabolicImageFilter() +{ + this->SetNumberOfRequiredOutputs(1); + this->SetNumberOfRequiredInputs(1); + this->m_CircErode = CircErodeType::New(); + this->m_CircDilate = CircDilateType::New(); + this->m_CircCastA = CCastTypeA::New(); + this->m_CircCastB = CCastTypeB::New(); + + this->m_RectErode = RectErodeType::New(); + this->m_RectDilate = RectDilateType::New(); + this->m_RectCastA = RCastTypeA::New(); + this->m_RectCastB = RCastTypeB::New(); + this->m_Circular = true; + // Need to call this after filters are created + this->SetUseImageSpacing(false); + this->SetSafeBorder(true); +} + +template +void +BinaryOpenParabolicImageFilter::SetRadius(ScalarRealType radius) +{ + RadiusType s; + + s.Fill(radius); + this->SetRadius(s); +} + +template +void +BinaryOpenParabolicImageFilter::GenerateData() +{ + // Allocate the output + this->AllocateOutputs(); + typename TInputImage::SizeType Pad; + + // numerical errors do seem to build up, so we need a margin on the + // thresholding steps. + // ScalarRealType margin = 0.0; + + // ScalarRealType mxRad = (ScalarRealType)(*std::max_element(m_Radius.Begin(), + // m_Radius.End())); + // // this needs to be examined more closely + // margin = 1.0/(pow(mxRad, TInputImage::ImageDimension) * 10); + // margin = std::min(margin, 0.00001); + // set up the scaling before we pass control over to superclass + if (this->m_RectErode->GetUseImageSpacing()) + { + // radius is in mm - need to do an adjustment to make sure that we + // end up with an odd number of voxels for the radius + RadiusType R; + for (unsigned P = 0; P < InputImageType::ImageDimension; P++) + { + typename TInputImage::SpacingValueType tsp = this->GetInput()->GetSpacing()[P]; + + // int thisvox=(int)round(m_Radius[P]/this->GetInput()->GetSpacing()[P]); + // if (thisvox % 2 == 0) ++thisvox; + // std::cout << thisvox << std::endl; + // float thisRad = thisvox * this->GetInput()->GetSpacing()[P]; + // R[P] = 0.5 * thisRad * thisRad + + // this->GetInput()->GetSpacing()[P]; + R[P] = 0.5 * (m_Radius[P] * m_Radius[P]) + tsp * tsp; + Pad[P] = (typename TInputImage::SizeType::SizeValueType)(itk::Math::rnd_halfinttoeven(m_Radius[P] / tsp) + 2); + } + m_RectErode->SetScale(R); + m_CircErode->SetScale(R); + m_RectDilate->SetScale(R); + m_CircDilate->SetScale(R); + } + else + { + // radius is in pixels + RadiusType R; + // this gives us a little bit of a margin + for (unsigned P = 0; P < InputImageType::ImageDimension; P++) + { + R[P] = (0.5 * m_Radius[P] * m_Radius[P] + 1); + Pad[P] = (typename TInputImage::SizeType::SizeValueType)m_Radius[P] + 1; + } + // std::cout << "no image spacing " << m_Radius << R << std::endl; + std::cout << Pad << R << std::endl; + m_RectErode->SetScale(R); + m_CircErode->SetScale(R); + m_RectDilate->SetScale(R); + m_CircDilate->SetScale(R); + } + + if (m_Circular) + { + ProgressAccumulator::Pointer progress = ProgressAccumulator::New(); + progress->SetMiniPipelineFilter(this); + InputImageConstPointer inputImage; + inputImage = this->GetInput(); + + progress->RegisterInternalFilter(m_CircErode, 0.4f); + progress->RegisterInternalFilter(m_CircCastA, 0.1f); + progress->RegisterInternalFilter(m_CircDilate, 0.4f); + progress->RegisterInternalFilter(m_CircCastB, 0.1f); + + m_CircCastA->SetInput(m_CircErode->GetOutput()); + // m_CircCastA->SetVal(1.0 - margin); + m_CircCastA->SetVal(1.0); + + m_CircDilate->SetInput(m_CircCastA->GetOutput()); + + m_CircCastB->SetInput(m_CircDilate->GetOutput()); + m_CircCastB->SetUpperThreshold(0.0); + m_CircCastB->SetOutsideValue(1); + m_CircCastB->SetInsideValue(0); + + if (m_SafeBorder) + { + using PadType = typename itk::ConstantPadImageFilter; + typename PadType::Pointer pad = PadType::New(); + pad->SetPadLowerBound(Pad); + pad->SetPadUpperBound(Pad); + pad->SetConstant(1); + pad->SetInput(inputImage); + m_CircErode->SetInput(pad->GetOutput()); + using CropType = typename itk::CropImageFilter; + typename CropType::Pointer crop = CropType::New(); + crop->SetInput(m_CircCastB->GetOutput()); + crop->SetUpperBoundaryCropSize(Pad); + crop->SetLowerBoundaryCropSize(Pad); + + crop->GraftOutput(this->GetOutput()); + crop->Update(); + this->GraftOutput(crop->GetOutput()); + } + else + { + m_CircErode->SetInput(inputImage); + + m_CircCastB->GraftOutput(this->GetOutput()); + m_CircCastB->Update(); + this->GraftOutput(m_CircCastB->GetOutput()); + } + } + else + { + ProgressAccumulator::Pointer progress = ProgressAccumulator::New(); + progress->SetMiniPipelineFilter(this); + InputImageConstPointer inputImage; + inputImage = this->GetInput(); + + progress->RegisterInternalFilter(m_RectErode, 0.4f); + progress->RegisterInternalFilter(m_RectCastA, 0.1f); + progress->RegisterInternalFilter(m_RectDilate, 0.4f); + progress->RegisterInternalFilter(m_RectCastB, 0.1f); + + m_RectCastA->SetInput(m_RectErode->GetOutput()); + m_RectCastA->SetVal(1); + + m_RectDilate->SetInput(m_RectCastA->GetOutput()); + + m_RectCastB->SetInput(m_RectDilate->GetOutput()); + m_RectCastB->SetUpperThreshold(0); + m_RectCastB->SetOutsideValue(1); + m_RectCastB->SetInsideValue(0); + + if (m_SafeBorder) + { + using PadType = typename itk::ConstantPadImageFilter; + typename PadType::Pointer pad = PadType::New(); + pad->SetPadLowerBound(Pad); + pad->SetPadUpperBound(Pad); + pad->SetConstant(1); + pad->SetInput(inputImage); + m_RectErode->SetInput(pad->GetOutput()); + + using CropType = typename itk::CropImageFilter; + typename CropType::Pointer crop = CropType::New(); + crop->SetInput(m_RectCastB->GetOutput()); + crop->SetUpperBoundaryCropSize(Pad); + crop->SetLowerBoundaryCropSize(Pad); + + crop->GraftOutput(this->GetOutput()); + crop->Update(); + this->GraftOutput(crop->GetOutput()); + } + else + { + m_RectErode->SetInput(inputImage); + m_RectCastB->GraftOutput(this->GetOutput()); + m_RectCastB->Update(); + + this->GraftOutput(m_RectCastB->GetOutput()); + } + } +} + +template +void +BinaryOpenParabolicImageFilter::PrintSelf(std::ostream & os, Indent indent) const +{ + Superclass::PrintSelf(os, indent); + if (this->m_Circular) + { + os << "Circular opening, "; + } + else + { + os << "Rectangular opening, "; + } + + if (this->m_SafeBorder) + { + os << "safe border" << std::endl; + } + else + { + os << "unsafe border" << std::endl; + } + + if (this->m_CircErode->GetUseImageSpacing()) + { + os << "Radius in world units: " << this->GetRadius() << std::endl; + } + else + { + os << "Radius in voxels: " << this->GetRadius() << std::endl; + } +} +} // namespace itk +#endif diff --git a/Modules/Filtering/ParabolicMorphology/include/itkGreaterEqualValImageFilter.h b/Modules/Filtering/ParabolicMorphology/include/itkGreaterEqualValImageFilter.h new file mode 100644 index 00000000000..b13426da017 --- /dev/null +++ b/Modules/Filtering/ParabolicMorphology/include/itkGreaterEqualValImageFilter.h @@ -0,0 +1,108 @@ +/*========================================================================= + * + * Copyright NumFOCUS + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * https://www.apache.org/licenses/LICENSE-2.0.txt + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + *=========================================================================*/ +#ifndef itkGreaterEqualValImageFilter_h +#define itkGreaterEqualValImageFilter_h +#include +namespace itk +{ +/** \class GreaterEqualValImageFilter + * \brief Computes the absolute difference between an image and a + * constant. Can be done with ShiftScale and AbsIamgeFilters. + * + * \ingroup ParabolicMorphology + * + */ + +namespace Functor +{ +template +class GEConst +{ +public: + GEConst() { m_Val = (TInput)0.0; } + void + SetVal(const TInput i) + { + m_Val = i; + } + + ~GEConst() = default; + bool + operator!=(const GEConst &) const + { + return false; + } + + bool + operator==(const GEConst & other) const + { + return !(*this != other); + } + + inline TOutput + operator()(const TInput & A) + { + return static_cast(A >= m_Val); + } + +private: + TInput m_Val; +}; +} // namespace Functor + +template +class ITK_TEMPLATE_EXPORT GreaterEqualValImageFilter + : public UnaryFunctorImageFilter> +{ +public: + ITK_DISALLOW_COPY_AND_MOVE(GreaterEqualValImageFilter); + + /** Standard class type alias. */ + using Self = GreaterEqualValImageFilter; + using Superclass = + UnaryFunctorImageFilter>; + using Pointer = SmartPointer; + using ConstPointer = SmartPointer; + + /** Method for creation through the object factory. */ + itkNewMacro(Self); + + void + SetVal(typename TInputImage::PixelType val) + { + this->GetFunctor().SetVal(val); + this->Modified(); + } + +#ifdef ITK_USE_CONCEPT_CHECKING + /** Begin concept checking */ + itkConceptMacro(InputConvertibleToDoubleCheck, (Concept::Convertible)); + itkConceptMacro(DoubleConvertibleToOutputCheck, (Concept::Convertible)); + /** End concept checking */ +#endif +protected: + GreaterEqualValImageFilter() = default; + ~GreaterEqualValImageFilter() override = default; +}; +} // end namespace itk + +#endif diff --git a/Modules/Filtering/ParabolicMorphology/include/itkMorphSDTHelperImageFilter.h b/Modules/Filtering/ParabolicMorphology/include/itkMorphSDTHelperImageFilter.h new file mode 100644 index 00000000000..372d6751e12 --- /dev/null +++ b/Modules/Filtering/ParabolicMorphology/include/itkMorphSDTHelperImageFilter.h @@ -0,0 +1,141 @@ +/*========================================================================= + * + * Copyright NumFOCUS + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * https://www.apache.org/licenses/LICENSE-2.0.txt + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + *=========================================================================*/ +#ifndef itkMorphSDTHelperImageFilter_h +#define itkMorphSDTHelperImageFilter_h + +#include "itkTernaryFunctorImageFilter.h" +#include "itkMath.h" + +namespace itk +{ +/** \class MorphSDTHelperImageFilter + * \brief Implements a pixel-wise operator to form a signed distance transform. + * + * Numeric conversions (castings) are done by the C++ defaults. + * + * \ingroup ParabolicMorphology + * + * \ingroup IntensityImageFilters Multithreaded + */ +namespace Function +{ +template +class MorphSDTHelper +{ +public: + MorphSDTHelper() {} + ~MorphSDTHelper() {} + void + SetVal(double i) + { + m_Val = i; + } + bool + operator!=(const MorphSDTHelper &) const + { + return false; + } + + bool + operator==(const MorphSDTHelper & other) const + { + return !(*this != other); + } + + inline TOutput + operator()(const TInput1 & A, const TInput2 & B, const TInput3 & C) + { + // A should be the output of the erosion, B the dilation, C the mask + if (C > 0) + { + // inside the mask + return static_cast(std::sqrt((double)A + m_Val)); + } + else + { + // outside the mask + return static_cast(-std::sqrt(m_Val - (double)B)); + } + } + +private: + double m_Val; +}; +} // namespace Function + +template +class ITK_TEMPLATE_EXPORT MorphSDTHelperImageFilter + : public TernaryFunctorImageFilter> +{ +public: + ITK_DISALLOW_COPY_AND_MOVE(MorphSDTHelperImageFilter); + + /** Standard class type alias. */ + using Self = MorphSDTHelperImageFilter; + using Superclass = TernaryFunctorImageFilter>; + + using Pointer = SmartPointer; + using ConstPointer = SmartPointer; + + /** Method for creation through the object factory. */ + itkNewMacro(Self); + + /** Runtime information support. */ + itkTypeMacro(MorphSDTHelperImageFilter, TernaryFunctorImageFilter); + + void + SetVal(double val) + { + this->GetFunctor().SetVal(val); + this->Modified(); + } + +#ifdef ITK_USE_CONCEPT_CHECKING + /** Begin concept checking */ + itkConceptMacro(Input1ConvertibleToOutputCheck, + (Concept::Convertible)); + itkConceptMacro(Input2ConvertibleToOutputCheck, + (Concept::Convertible)); + itkConceptMacro(GreaterThanComparable, + (Concept::GreaterThanComparable)); + /** End concept checking */ +#endif +protected: + MorphSDTHelperImageFilter() {} + virtual ~MorphSDTHelperImageFilter() {} +}; +} // end namespace itk + +#endif diff --git a/Modules/Filtering/ParabolicMorphology/include/itkMorphologicalDistanceTransformImageFilter.h b/Modules/Filtering/ParabolicMorphology/include/itkMorphologicalDistanceTransformImageFilter.h new file mode 100644 index 00000000000..2e0c66e4a3c --- /dev/null +++ b/Modules/Filtering/ParabolicMorphology/include/itkMorphologicalDistanceTransformImageFilter.h @@ -0,0 +1,159 @@ +/*========================================================================= + * + * Copyright NumFOCUS + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * https://www.apache.org/licenses/LICENSE-2.0.txt + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + *=========================================================================*/ +#ifndef itkMorphologicalDistanceTransformImageFilter_h +#define itkMorphologicalDistanceTransformImageFilter_h + +#include "itkImageToImageFilter.h" +#include "itkProgressReporter.h" + +#include "itkBinaryThresholdImageFilter.h" +#include "itkParabolicErodeImageFilter.h" +#include "itkSqrtImageFilter.h" + +namespace itk +{ +/** + * \class MorphologicalDistanceTransformImageFilter + * \brief Distance transform of a mask using parabolic morphological + * methods + * + * Morphological erosions using a parabolic structuring element can be + * used to compute a distance transform of a mask by setting the + * "Outside" value to 0 and the "inside" value to +infinity. The + * output of the parabolic erosion are the square of the distance to + * the nearest zero valued voxel. Thus we can compute the distance + * transform by taking the sqrt of the erosion. + * + * The output pixel type needs to support values as large as the + * square of the largest value of the distance - just use float to be + * safe. + * + * Core methods described in the InsightJournal article: + * "Morphology with parabolic structuring elements" + * + * https://hdl.handle.net/1926/1370 + * + * \ingroup ParabolicMorphology + * + * \author Richard Beare, Monash University, Department of Medicine, + * Melbourne, Australia. + * + **/ + +template +class ITK_TEMPLATE_EXPORT MorphologicalDistanceTransformImageFilter + : public ImageToImageFilter +{ +public: + ITK_DISALLOW_COPY_AND_MOVE(MorphologicalDistanceTransformImageFilter); + + /** Standard class type alias. */ + using Self = MorphologicalDistanceTransformImageFilter; + using Superclass = ImageToImageFilter; + using Pointer = SmartPointer; + using ConstPointer = SmartPointer; + + /** Method for creation through the object factory. */ + itkNewMacro(Self); + + /** Runtime information support. */ + itkTypeMacro(MorphologicalDistanceTransformImageFilter, ImageToImageFilter); + + /** Pixel Type of the input image */ + using InputImageType = TInputImage; + using OutputImageType = TOutputImage; + using InputPixelType = typename TInputImage::PixelType; + using RealType = typename NumericTraits::RealType; + using ScalarRealType = typename NumericTraits::ScalarRealType; + using OutputPixelType = typename TOutputImage::PixelType; + + /** Smart pointer type alias support. */ + using InputImagePointer = typename TInputImage::Pointer; + using InputImageConstPointer = typename TInputImage::ConstPointer; + + /** Image related type alias. */ + static constexpr unsigned int OutputImageDimension = TOutputImage::ImageDimension; + static constexpr unsigned int InputImageDimension = TInputImage::ImageDimension; + static constexpr unsigned int ImageDimension = TInputImage::ImageDimension; + + /** a type to represent the "kernel radius" */ + using RadiusType = typename itk::FixedArray; + void + Modified() const override; + + /** this describes the input mask - default value 0 - we compute the + distance from all voxels with value not equal to "OutsideValue" to + the nearest voxel with value "OutsideValue" */ + itkSetMacro(OutsideValue, InputPixelType); + itkGetConstReferenceMacro(OutsideValue, InputPixelType); + + /** Is the transform in world or voxel units - default is world */ + void + SetUseImageSpacing(bool uis) + { + m_Erode->SetUseImageSpacing(uis); + } + + const bool & + GetUseImageSpacing() + { + return m_Erode->GetUseImageSpacing(); + } + + itkSetMacro(SqrDist, bool); + itkGetConstReferenceMacro(SqrDist, bool); + itkBooleanMacro(SqrDist); + +#ifdef ITK_USE_CONCEPT_CHECKING + /** Begin concept checking */ + itkConceptMacro(SameDimension, + (Concept::SameDimension)); + + itkConceptMacro(Comparable, (Concept::Comparable)); + + /** End concept checking */ +#endif +protected: + MorphologicalDistanceTransformImageFilter(); + ~MorphologicalDistanceTransformImageFilter() override = default; + void + PrintSelf(std::ostream & os, Indent indent) const override; + + /** Generate Data */ + void + GenerateData() override; + + // do everything in the output image type, which should have high precision + using ThreshType = typename itk::BinaryThresholdImageFilter; + using ErodeType = typename itk::ParabolicErodeImageFilter; + using SqrtType = typename itk::SqrtImageFilter; + +private: + InputPixelType m_OutsideValue; + typename ErodeType::Pointer m_Erode; + typename ThreshType::Pointer m_Thresh; + typename SqrtType::Pointer m_Sqrt; + bool m_SqrDist; +}; +} // namespace itk +#ifndef ITK_MANUAL_INSTANTIATION +# include "itkMorphologicalDistanceTransformImageFilter.hxx" +#endif + +#endif diff --git a/Modules/Filtering/ParabolicMorphology/include/itkMorphologicalDistanceTransformImageFilter.hxx b/Modules/Filtering/ParabolicMorphology/include/itkMorphologicalDistanceTransformImageFilter.hxx new file mode 100644 index 00000000000..77691e196ed --- /dev/null +++ b/Modules/Filtering/ParabolicMorphology/include/itkMorphologicalDistanceTransformImageFilter.hxx @@ -0,0 +1,136 @@ +/*========================================================================= + * + * Copyright NumFOCUS + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * https://www.apache.org/licenses/LICENSE-2.0.txt + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + *=========================================================================*/ +#ifndef itkMorphologicalDistanceTransformImageFilter_hxx +#define itkMorphologicalDistanceTransformImageFilter_hxx + +#include "itkProgressAccumulator.h" + +namespace itk +{ +template +MorphologicalDistanceTransformImageFilter::MorphologicalDistanceTransformImageFilter() +{ + this->SetNumberOfRequiredOutputs(1); + this->SetNumberOfRequiredInputs(1); + + m_Erode = ErodeType::New(); + m_Thresh = ThreshType::New(); + m_Sqrt = SqrtType::New(); + m_OutsideValue = 0; + m_Erode->SetScale(0.5); + this->SetUseImageSpacing(true); + m_SqrDist = false; +} + +template +void +MorphologicalDistanceTransformImageFilter::Modified() const +{ + Superclass::Modified(); + m_Erode->Modified(); + m_Thresh->Modified(); + m_Sqrt->Modified(); +} + +template +void +MorphologicalDistanceTransformImageFilter::GenerateData() +{ + ProgressAccumulator::Pointer progress = ProgressAccumulator::New(); + + progress->SetMiniPipelineFilter(this); + // these values are guesses at present - need to profile to get a + // real idea + progress->RegisterInternalFilter(m_Thresh, 0.1f); + progress->RegisterInternalFilter(m_Erode, 0.8f); + progress->RegisterInternalFilter(m_Sqrt, 0.1f); + + // std::cout << "DT" << std::endl; + + double MaxDist = 0.0; + typename TOutputImage::SpacingType sp = this->GetOutput()->GetSpacing(); + typename TOutputImage::SizeType sz = this->GetOutput()->GetLargestPossibleRegion().GetSize(); + if (this->GetUseImageSpacing()) + { + for (unsigned k = 0; k < TOutputImage::ImageDimension; k++) + { + double thisdim = (sz[k] * sp[k]); + MaxDist += thisdim * thisdim; + } + } + else + { + for (unsigned k = 0; k < TOutputImage::ImageDimension; k++) + { + double thisdim = sz[k]; + MaxDist += thisdim * thisdim; + } + } + + // double Wt = 0.0; + // if (this->GetUseImageSpacing()) + // { + // for (unsigned k = 0; k < TOutputImage::ImageDimension; k++) + // { + // Wt += sp[k] * sp[k]; + // } + // } + // else + // { + // for (unsigned k = 0; k < TOutputImage::ImageDimension; k++) + // { + // Wt += 1.0; + // } + // } + // Wt = sqrt(Wt); + this->AllocateOutputs(); + + m_Thresh->SetLowerThreshold(m_OutsideValue); + m_Thresh->SetUpperThreshold(m_OutsideValue); + m_Thresh->SetOutsideValue(MaxDist); + m_Thresh->SetInsideValue(0); + + m_Thresh->SetInput(this->GetInput()); + m_Erode->SetInput(m_Thresh->GetOutput()); + + if (m_SqrDist) + { + m_Erode->GraftOutput(this->GetOutput()); + m_Erode->Update(); + this->GraftOutput(m_Erode->GetOutput()); + } + else + { + m_Sqrt->SetInput(m_Erode->GetOutput()); + m_Sqrt->GraftOutput(this->GetOutput()); + m_Sqrt->Update(); + this->GraftOutput(m_Sqrt->GetOutput()); + } +} + +template +void +MorphologicalDistanceTransformImageFilter::PrintSelf(std::ostream & os, Indent indent) const +{ + Superclass::PrintSelf(os, indent); + os << "Outside Value = " << (OutputPixelType)m_OutsideValue << std::endl; + os << "ImageScale = " << m_Erode->GetUseImageSpacing() << std::endl; +} +} // namespace itk + +#endif diff --git a/Modules/Filtering/ParabolicMorphology/include/itkMorphologicalSharpeningImageFilter.h b/Modules/Filtering/ParabolicMorphology/include/itkMorphologicalSharpeningImageFilter.h new file mode 100644 index 00000000000..bda9f1e7898 --- /dev/null +++ b/Modules/Filtering/ParabolicMorphology/include/itkMorphologicalSharpeningImageFilter.h @@ -0,0 +1,182 @@ +/*========================================================================= + * + * Copyright NumFOCUS + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * https://www.apache.org/licenses/LICENSE-2.0.txt + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + *=========================================================================*/ +#ifndef itkMorphologicalSharpeningImageFilter_h +#define itkMorphologicalSharpeningImageFilter_h + +#include "itkImageToImageFilter.h" +// #include "itkProgressReporter.h" +#include "itkCastImageFilter.h" +#include "itkParabolicErodeImageFilter.h" +#include "itkParabolicDilateImageFilter.h" +#include "itkSharpenOpImageFilter.h" + +namespace itk +{ +/** + * \class MorphologicalSharpeningImageFilter + * \brief Image sharpening using methods based on parabolic + * structuring elements. + * + * This is an implemtentation of the method of Schavemaker for testing + * the parabolic morphology routines. No particular efforts have been + * made to minimize memory consumption. + * + * + * \@article{Schavemaker2000, + * author = {Schavemaker, J. and Reinders, M. and Gerbrands, J. and Backer, E. +}, + * title = {Image sharpening by morphological filtering}, + * journal = {Pattern Recognition}, + * volume = {33}, + * number = {6}, + * year = {2000}, + * pages = {997-1012}, + * ee = {https://dx.doi.org/10.1016/S0031-3203(99)00160-0}, + * bibsource = {DBLP, https://dblp.uni-trier.de} +} + + * Core methods described in the InsightJournal article: + * "Morphology with parabolic structuring elements" + * + * https://hdl.handle.net/1926/1370 + * + * + * \ingroup ParabolicMorphology + * + * \author Richard Beare, Monash University, Department of Medicine, + * Melbourne, Australia. + * +**/ + +template +class ITK_TEMPLATE_EXPORT MorphologicalSharpeningImageFilter : public ImageToImageFilter +{ +public: + ITK_DISALLOW_COPY_AND_MOVE(MorphologicalSharpeningImageFilter); + + /** Standard class type alias. */ + using Self = MorphologicalSharpeningImageFilter; + using Superclass = ImageToImageFilter; + using Pointer = SmartPointer; + using ConstPointer = SmartPointer; + + /** Method for creation through the object factory. */ + itkNewMacro(Self); + + /** Runtime information support. */ + itkTypeMacro(MorphologicalSharpeningImageFilter, ImageToImageFilter); + + /** Pixel Type of the input image */ + using InputImageType = TInputImage; + using OutputImageType = TOutputImage; + using InputPixelType = typename TInputImage::PixelType; + using RealType = typename NumericTraits::RealType; + using ScalarRealType = typename NumericTraits::ScalarRealType; + using OutputPixelType = typename TOutputImage::PixelType; + + /** Smart pointer type alias support. */ + using InputImagePointer = typename TInputImage::Pointer; + using InputImageConstPointer = typename TInputImage::ConstPointer; + + /** Image related type alias. */ + static constexpr unsigned int OutputImageDimension = TOutputImage::ImageDimension; + static constexpr unsigned int InputImageDimension = TInputImage::ImageDimension; + static constexpr unsigned int ImageDimension = TInputImage::ImageDimension; + + /** a type to represent the "kernel radius" */ + using RadiusType = typename itk::FixedArray; + + itkSetMacro(Iterations, int); + itkGetConstReferenceMacro(Iterations, int); + + void + SetScale(ScalarRealType scale) + { + m_Erode->SetScale(scale); + m_Dilate->SetScale(scale); + } + + void + SetScale(RadiusType scale) + { + m_Erode->SetScale(scale); + m_Dilate->SetScale(scale); + } + + void + SetUseImageSpacing(bool uis) + { + m_Erode->SetUseImageSpacing(uis); + m_Dilate->SetUseImageSpacing(uis); + } + + // need to include the Get methods + const RadiusType & + GetScale() + { + return m_Erode->GetScale(); + } + + const bool & + GetUseImageSpacing() + { + return m_Erode->GetUseImageSpacing(); + } + +#ifdef ITK_USE_CONCEPT_CHECKING + /** Begin concept checking */ + itkConceptMacro(SameDimension, + (Concept::SameDimension)); + + itkConceptMacro(Comparable, (Concept::Comparable)); + + /** End concept checking */ +#endif +protected: + MorphologicalSharpeningImageFilter(); + ~MorphologicalSharpeningImageFilter() override = default; + void + PrintSelf(std::ostream & os, Indent indent) const override; + + /** Generate Data */ + void + GenerateData() override; + + // do everything in the output image type, which should have high precision + using ErodeType = typename itk::ParabolicErodeImageFilter; + using DilateType = typename itk::ParabolicDilateImageFilter; + using CastType = typename itk::CastImageFilter; + + using SharpenOpType = + typename itk::SharpenOpImageFilter; + +private: + int m_Iterations; + + typename ErodeType::Pointer m_Erode; + typename DilateType::Pointer m_Dilate; + typename CastType::Pointer m_Cast; + typename SharpenOpType::Pointer m_SharpenOp; +}; +} // namespace itk +#ifndef ITK_MANUAL_INSTANTIATION +# include "itkMorphologicalSharpeningImageFilter.hxx" +#endif + +#endif diff --git a/Modules/Filtering/ParabolicMorphology/include/itkMorphologicalSharpeningImageFilter.hxx b/Modules/Filtering/ParabolicMorphology/include/itkMorphologicalSharpeningImageFilter.hxx new file mode 100644 index 00000000000..2129dc11eea --- /dev/null +++ b/Modules/Filtering/ParabolicMorphology/include/itkMorphologicalSharpeningImageFilter.hxx @@ -0,0 +1,102 @@ +/*========================================================================= + * + * Copyright NumFOCUS + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * https://www.apache.org/licenses/LICENSE-2.0.txt + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + *=========================================================================*/ +#ifndef itkMorphologicalSharpeningImageFilter_hxx +#define itkMorphologicalSharpeningImageFilter_hxx + +#include "itkProgressAccumulator.h" + +namespace itk + +{ +template +MorphologicalSharpeningImageFilter::MorphologicalSharpeningImageFilter() +{ + this->SetNumberOfRequiredOutputs(1); + this->SetNumberOfRequiredInputs(1); + + m_Erode = ErodeType::New(); + m_Dilate = DilateType::New(); + m_Cast = CastType::New(); + m_SharpenOp = SharpenOpType::New(); + m_Iterations = 1; + this->SetScale(1); + this->SetUseImageSpacing(false); +} + +template +void +MorphologicalSharpeningImageFilter::GenerateData() +{ + ProgressAccumulator::Pointer progress = ProgressAccumulator::New(); + + progress->SetMiniPipelineFilter(this); + + // Allocate the output + this->AllocateOutputs(); + + InputImageConstPointer inputImage = this->GetInput(); + m_Cast->SetInput(inputImage); + + // set the input to the morph operations + m_Erode->SetInput(m_Cast->GetOutput()); + m_Dilate->SetInput(m_Cast->GetOutput()); + m_SharpenOp->SetInput(m_Dilate->GetOutput()); + m_SharpenOp->SetInput2(m_Cast->GetOutput()); + m_SharpenOp->SetInput3(m_Erode->GetOutput()); + + progress->RegisterInternalFilter(m_Erode, 1.0f); + progress->RegisterInternalFilter(m_Dilate, 1.0f); + progress->RegisterInternalFilter(m_SharpenOp, 1.0f); + + // set up the progrss monitor + // WatershedMiniPipelineProgressCommand::Pointer c = + // WatershedMiniPipelineProgressCommand::New(); + // c->SetFilter(this); + // c->SetNumberOfFilters(3 * m_Iterations); + + // m_Erode->AddObserver(ProgressEvent(), c); + // m_Dilate->AddObserver(ProgressEvent(), c); + // m_SharpenOp->AddObserver(ProgressEvent(), c); + + for (int i = 0; i < m_Iterations; i++) + { + if (i != 0) + { + m_Erode->SetInput(this->GetOutput()); + m_Dilate->SetInput(this->GetOutput()); + m_SharpenOp->SetInput2(this->GetOutput()); + m_Erode->Modified(); + m_Dilate->Modified(); + } + + m_SharpenOp->GraftOutput(this->GetOutput()); + m_SharpenOp->Update(); + this->GraftOutput(m_SharpenOp->GetOutput()); + } +} + +template +void +MorphologicalSharpeningImageFilter::PrintSelf(std::ostream & os, Indent indent) const +{ + Superclass::PrintSelf(os, indent); + os << "Iterations = " << m_Iterations << std::endl; +} +} // end namespace itk + +#endif diff --git a/Modules/Filtering/ParabolicMorphology/include/itkMorphologicalSignedDistanceTransformImageFilter.h b/Modules/Filtering/ParabolicMorphology/include/itkMorphologicalSignedDistanceTransformImageFilter.h new file mode 100644 index 00000000000..5be90086c09 --- /dev/null +++ b/Modules/Filtering/ParabolicMorphology/include/itkMorphologicalSignedDistanceTransformImageFilter.h @@ -0,0 +1,206 @@ +/*========================================================================= + * + * Copyright NumFOCUS + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * https://www.apache.org/licenses/LICENSE-2.0.txt + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + *=========================================================================*/ +#ifndef itkMorphologicalSignedDistanceTransformImageFilter_h +#define itkMorphologicalSignedDistanceTransformImageFilter_h + +#include "itkImageToImageFilter.h" +#include "itkProgressReporter.h" + +#include "itkBinaryThresholdImageFilter.h" +#include "itkParabolicErodeImageFilter.h" +#include "itkParabolicDilateImageFilter.h" +#include "itkMorphSDTHelperImageFilter.h" + +namespace itk +{ +/** + * \class MorphologicalSignedDistanceTransformImageFilter + * \brief Signed distance transform of a mask using parabolic morphological + * methods + * + * Morphological erosions using a parabolic structuring element can be + * used to compute a distance transform of a mask by setting the + * "Outside" value to 0 and the "inside" value to +infinity (or beyond + * the maximum possible value). The + * output of the parabolic erosion are the square of the distance to + * the nearest zero valued voxel. Thus we can compute the distance + * transform by taking the sqrt of the erosion. + * + * The output pixel type needs to support values as large as the + * square of the largest value of the distance - just use float to be + * safe. + * + * The inside is considered to have negative distances. Use + * InsideIsPositive(bool) to change. + * + * There are also OutsideValue methods which can be used in similar + * ways. + * + * Otherwise it is meant to have an interface to the other + * DistanceTransforms filters. + * + * Core methods described in the InsightJournal article: + * "Morphology with parabolic structuring elements" + * + * https://hdl.handle.net/1926/1370 + * + * + * \ingroup ParabolicMorphology + * + * \author Richard Beare, Monash University, Department of Medicine, + * Melbourne, Australia. + * + **/ + +template +class ITK_TEMPLATE_EXPORT MorphologicalSignedDistanceTransformImageFilter + : public ImageToImageFilter +{ +public: + ITK_DISALLOW_COPY_AND_MOVE(MorphologicalSignedDistanceTransformImageFilter); + + /** Standard class type alias. */ + using Self = MorphologicalSignedDistanceTransformImageFilter; + using Superclass = ImageToImageFilter; + using Pointer = SmartPointer; + using ConstPointer = SmartPointer; + + /** Method for creation through the object factory. */ + itkNewMacro(Self); + + /** Runtime information support. */ + itkTypeMacro(MorphologicalSignedDistanceTransformImageFilter, ImageToImageFilter); + + /** Pixel Type of the input image */ + using InputImageType = TInputImage; + using OutputImageType = TOutputImage; + using InputPixelType = typename TInputImage::PixelType; + using RealType = typename NumericTraits::RealType; + using ScalarRealType = typename NumericTraits::ScalarRealType; + using OutputPixelType = typename TOutputImage::PixelType; + + /** Smart pointer type alias support. */ + using InputImagePointer = typename TInputImage::Pointer; + using InputImageConstPointer = typename TInputImage::ConstPointer; + + /** Image related type alias. */ + static constexpr unsigned int OutputImageDimension = TOutputImage::ImageDimension; + static constexpr unsigned int InputImageDimension = TInputImage::ImageDimension; + static constexpr unsigned int ImageDimension = TInputImage::ImageDimension; + + /** a type to represent the "kernel radius" */ + using RadiusType = typename itk::FixedArray; + + virtual void + Modified() const; + + /** this describes the input mask - default value 0 - we compute the + distance from all voxels with value not equal to "OutsideValue" to + the nearest voxel with value "OutsideValue" */ + itkSetMacro(OutsideValue, InputPixelType); + itkGetConstReferenceMacro(OutsideValue, InputPixelType); + + /** Set On/Off whether spacing is used. */ + itkBooleanMacro(UseImageSpacing); + + /** Set if the inside represents positive values in the signed distance + * map. By convention ON pixels are treated as inside pixels. */ + itkSetMacro(InsideIsPositive, bool); + + /** Get if the inside represents positive values in the signed distance map. + * See GetInsideIsPositive() */ + itkGetConstReferenceMacro(InsideIsPositive, bool); + + /** Set if the inside represents positive values in the signed distance + * map. By convention ON pixels are treated as inside pixels. Default is + * true. */ + itkBooleanMacro(InsideIsPositive); + /** Is the transform in world or voxel units - default is world */ + void + SetUseImageSpacing(bool uis) + { + m_Erode->SetUseImageSpacing(uis); + m_Dilate->SetUseImageSpacing(uis); + this->Modified(); + } + + enum ParabolicAlgorithm + { + NOCHOICE = 0, // decices based on scale - experimental + CONTACTPOINT = 1, // sometimes faster at low scale + INTERSECTION = 2 // default + }; + + /** + * Set/Get the method used. Choices are contact point or + * intersection. Intersection is the default. Contact point can be + * faster at small scales. This is very unlikely to be the case for + * a distance transform. + */ + + itkSetMacro(ParabolicAlgorithm, int); + itkGetConstReferenceMacro(ParabolicAlgorithm, int); + + const bool & + GetUseImageSpacing() + { + return m_Erode->GetUseImageSpacing(); + } + +#ifdef ITK_USE_CONCEPT_CHECKING + /** Begin concept checking */ + itkConceptMacro(SameDimension, + (Concept::SameDimension)); + + itkConceptMacro(Comparable, (Concept::Comparable)); + + /** End concept checking */ +#endif +protected: + MorphologicalSignedDistanceTransformImageFilter(); + virtual ~MorphologicalSignedDistanceTransformImageFilter() {} + void + PrintSelf(std::ostream & os, Indent indent) const; + + /** Generate Data */ + void + GenerateData(void); + + int m_ParabolicAlgorithm; + + // do everything in the output image type, which should have high precision + using ThreshType = typename itk::BinaryThresholdImageFilter; + using ErodeType = typename itk::ParabolicErodeImageFilter; + using DilateType = typename itk::ParabolicDilateImageFilter; + using HelperType = typename itk::MorphSDTHelperImageFilter; + +private: + InputPixelType m_OutsideValue; + bool m_InsideIsPositive; + typename ErodeType::Pointer m_Erode; + typename DilateType::Pointer m_Dilate; + typename ThreshType::Pointer m_Thresh; + typename HelperType::Pointer m_Helper; +}; +} // namespace itk +#ifndef ITK_MANUAL_INSTANTIATION +# include "itkMorphologicalSignedDistanceTransformImageFilter.hxx" +#endif + +#endif diff --git a/Modules/Filtering/ParabolicMorphology/include/itkMorphologicalSignedDistanceTransformImageFilter.hxx b/Modules/Filtering/ParabolicMorphology/include/itkMorphologicalSignedDistanceTransformImageFilter.hxx new file mode 100644 index 00000000000..f3ce8ee9b24 --- /dev/null +++ b/Modules/Filtering/ParabolicMorphology/include/itkMorphologicalSignedDistanceTransformImageFilter.hxx @@ -0,0 +1,138 @@ +/*========================================================================= + * + * Copyright NumFOCUS + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * https://www.apache.org/licenses/LICENSE-2.0.txt + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + *=========================================================================*/ +#ifndef itkMorphologicalSignedDistanceTransformImageFilter_hxx +#define itkMorphologicalSignedDistanceTransformImageFilter_hxx + +#include "itkProgressAccumulator.h" + +namespace itk +{ +template +MorphologicalSignedDistanceTransformImageFilter::MorphologicalSignedDistanceTransformImageFilter() +{ + this->SetNumberOfRequiredOutputs(1); + this->SetNumberOfRequiredInputs(1); + + m_Erode = ErodeType::New(); + m_Dilate = DilateType::New(); + m_Thresh = ThreshType::New(); + m_Helper = HelperType::New(); + m_Erode->SetScale(0.5); + m_Dilate->SetScale(0.5); + this->SetUseImageSpacing(true); + this->SetInsideIsPositive(false); + m_OutsideValue = 0; + m_ParabolicAlgorithm = INTERSECTION; +} + +template +void +MorphologicalSignedDistanceTransformImageFilter::Modified() const +{ + Superclass::Modified(); + m_Erode->Modified(); + m_Dilate->Modified(); + m_Thresh->Modified(); + m_Helper->Modified(); +} + +template +void +MorphologicalSignedDistanceTransformImageFilter::GenerateData(void) +{ + ProgressAccumulator::Pointer progress = ProgressAccumulator::New(); + + progress->SetMiniPipelineFilter(this); + // these values are guesses at present - need to profile to get a + // real idea + progress->RegisterInternalFilter(m_Thresh, 0.1f); + progress->RegisterInternalFilter(m_Erode, 0.4f); + progress->RegisterInternalFilter(m_Dilate, 0.4f); + progress->RegisterInternalFilter(m_Helper, 0.1f); + + m_Erode->SetParabolicAlgorithm(m_ParabolicAlgorithm); + m_Dilate->SetParabolicAlgorithm(m_ParabolicAlgorithm); + + this->AllocateOutputs(); + // figure out the maximum value of distance transform using the + // image dimensions + typename TOutputImage::SizeType sz = this->GetOutput()->GetRequestedRegion().GetSize(); + typename TOutputImage::SpacingType sp = this->GetOutput()->GetSpacing(); + + double MaxDist = 0.0; + if (this->GetUseImageSpacing()) + { + for (unsigned k = 0; k < TOutputImage::ImageDimension; k++) + { + double thisdim = (sz[k] * sp[k]); + MaxDist += thisdim * thisdim; + } + } + else + { + for (unsigned k = 0; k < TOutputImage::ImageDimension; k++) + { + double thisdim = sz[k]; + MaxDist += thisdim * thisdim; + } + } + + m_Thresh->SetLowerThreshold(m_OutsideValue); + m_Thresh->SetUpperThreshold(m_OutsideValue); + if (this->GetInsideIsPositive()) + { + m_Thresh->SetOutsideValue(MaxDist); + m_Thresh->SetInsideValue(-MaxDist); + } + else + { + m_Thresh->SetOutsideValue(-MaxDist); + m_Thresh->SetInsideValue(MaxDist); + } + + m_Thresh->SetInput(this->GetInput()); + m_Erode->SetInput(m_Thresh->GetOutput()); + m_Dilate->SetInput(m_Thresh->GetOutput()); +#if 1 + m_Helper->SetInput(m_Erode->GetOutput()); + m_Helper->SetInput2(m_Dilate->GetOutput()); + m_Helper->SetInput3(m_Thresh->GetOutput()); + m_Helper->SetVal(MaxDist); + m_Helper->GraftOutput(this->GetOutput()); + m_Helper->Update(); + this->GraftOutput(m_Helper->GetOutput()); +#else + m_Dilate->GraftOutput(this->GetOutput()); + m_Dilate->Update(); + this->GraftOutput(m_Dilate->GetOutput()); +#endif +} + +template +void +MorphologicalSignedDistanceTransformImageFilter::PrintSelf(std::ostream & os, + Indent indent) const +{ + Superclass::PrintSelf(os, indent); + os << "Outside Value = " << (OutputPixelType)m_OutsideValue << std::endl; + os << "ImageScale = " << m_Erode->GetUseImageSpacing() << std::endl; +} +} // namespace itk + +#endif diff --git a/Modules/Filtering/ParabolicMorphology/include/itkParabolicCloseImageFilter.h b/Modules/Filtering/ParabolicMorphology/include/itkParabolicCloseImageFilter.h new file mode 100644 index 00000000000..c51bc719c3b --- /dev/null +++ b/Modules/Filtering/ParabolicMorphology/include/itkParabolicCloseImageFilter.h @@ -0,0 +1,95 @@ +/*========================================================================= + * + * Copyright NumFOCUS + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * https://www.apache.org/licenses/LICENSE-2.0.txt + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + *=========================================================================*/ +#ifndef itkParabolicCloseImageFilter_h +#define itkParabolicCloseImageFilter_h + +#include "itkParabolicOpenCloseSafeBorderImageFilter.h" +#include "itkNumericTraits.h" + +namespace itk +{ +/** + * \class ParabolicCloseImageFilter + * \brief Class for morphological closing + * operations with parabolic structuring elements. + * + * This filter provides options for padded borders + * + * This filter is threaded. + * + * Core methods described in the InsightJournal article: + * "Morphology with parabolic structuring elements" + * + * https://hdl.handle.net/1926/1370 + * + * \sa itkParabolicOpenCloseImageFilter + * + * + * \ingroup ParabolicMorphology + * + * \author Richard Beare, Department of Medicine, Monash University, + * Australia. + **/ + +template +class ITK_TEMPLATE_EXPORT ParabolicCloseImageFilter + : public ParabolicOpenCloseSafeBorderImageFilter +{ +public: + ITK_DISALLOW_COPY_AND_MOVE(ParabolicCloseImageFilter); + + /** Standard class type alias. */ + using Self = ParabolicCloseImageFilter; + using Superclass = ParabolicOpenCloseSafeBorderImageFilter; + using Pointer = SmartPointer; + using ConstPointer = SmartPointer; + + /** Method for creation through the object factory. */ + itkNewMacro(Self); + + /** Runtime information support. */ + itkTypeMacro(ParabolicCloseImageFilter, ParabolicOpenCloseSafeBorderImageFilter); + + /** Pixel Type of the input image */ + using InputImageType = TInputImage; + using OutputImageType = TOutputImage; + using PixelType = typename TInputImage::PixelType; + using RealType = typename NumericTraits::RealType; + using ScalarRealType = typename NumericTraits::ScalarRealType; + using OutputPixelType = typename TOutputImage::PixelType; + + /** Smart pointer type alias support. */ + using InputImagePointer = typename TInputImage::Pointer; + using InputImageConstPointer = typename TInputImage::ConstPointer; + + /** a type to represent the "kernel radius" */ + using RadiusType = typename itk::FixedArray; + + /** Image dimension. */ + static constexpr unsigned int ImageDimension = TInputImage::ImageDimension; + /** Define the image type for internal computations + RealType is usually 'double' in NumericTraits. + Here we prefer float in order to save memory. */ +protected: + ParabolicCloseImageFilter() {} + virtual ~ParabolicCloseImageFilter() {} + // void PrintSelf(std::ostream& os, Indent indent) const; +}; +} // end namespace itk + +#endif diff --git a/Modules/Filtering/ParabolicMorphology/include/itkParabolicDilateImageFilter.h b/Modules/Filtering/ParabolicMorphology/include/itkParabolicDilateImageFilter.h new file mode 100644 index 00000000000..8d3022025f7 --- /dev/null +++ b/Modules/Filtering/ParabolicMorphology/include/itkParabolicDilateImageFilter.h @@ -0,0 +1,92 @@ +/*========================================================================= + * + * Copyright NumFOCUS + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * https://www.apache.org/licenses/LICENSE-2.0.txt + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + *=========================================================================*/ +#ifndef itkParabolicDilateImageFilter_h +#define itkParabolicDilateImageFilter_h + +#include "itkParabolicErodeDilateImageFilter.h" +#include "itkNumericTraits.h" + +namespace itk +{ +/** + * \class ParabolicDilateImageFilter + * \brief Class for morphological dilation + * operations with parabolic structuring elements. + * + * This filter is threaded. + * + * Core methods described in the InsightJournal article: + * "Morphology with parabolic structuring elements" + * + * https://hdl.handle.net/1926/1370 + * + * \sa itkParabolicOpenCloseImageFilter + * + * + * \ingroup ParabolicMorphology + * + * \author Richard Beare, Department of Medicine, Monash University, + * Australia. + * + **/ + +template +class ITK_TEMPLATE_EXPORT ParabolicDilateImageFilter + : public ParabolicErodeDilateImageFilter +{ +public: + ITK_DISALLOW_COPY_AND_MOVE(ParabolicDilateImageFilter); + + /** Standard class type alias. */ + using Self = ParabolicDilateImageFilter; + using Superclass = ParabolicErodeDilateImageFilter; + using Pointer = SmartPointer; + using ConstPointer = SmartPointer; + + /** Method for creation through the object factory. */ + itkNewMacro(Self); + + /** Runtime information support. */ + itkTypeMacro(ParabolicDilateImageFilter, ParabolicErodeDilateImageFilter); + + /** Pixel Type of the input image */ + using InputImageType = TInputImage; + using OutputImageType = TOutputImage; + using PixelType = typename TInputImage::PixelType; + using RealType = typename NumericTraits::RealType; + using ScalarRealType = typename NumericTraits::ScalarRealType; + using OutputPixelType = typename TOutputImage::PixelType; + + /** Smart pointer type alias support. */ + using InputImagePointer = typename TInputImage::Pointer; + using InputImageConstPointer = typename TInputImage::ConstPointer; + + /** a type to represent the "kernel radius" */ + using RadiusType = typename itk::FixedArray; + + /** Image dimension. */ + static constexpr unsigned int ImageDimension = TInputImage::ImageDimension; + +protected: + ParabolicDilateImageFilter() = default; + ~ParabolicDilateImageFilter() override = default; + // void PrintSelf(std::ostream& os, Indent indent) const; +}; +} // end namespace itk + +#endif diff --git a/Modules/Filtering/ParabolicMorphology/include/itkParabolicErodeDilateImageFilter.h b/Modules/Filtering/ParabolicMorphology/include/itkParabolicErodeDilateImageFilter.h new file mode 100644 index 00000000000..c90d88a0efe --- /dev/null +++ b/Modules/Filtering/ParabolicMorphology/include/itkParabolicErodeDilateImageFilter.h @@ -0,0 +1,211 @@ +/*========================================================================= + * + * Copyright NumFOCUS + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * https://www.apache.org/licenses/LICENSE-2.0.txt + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + *=========================================================================*/ +#ifndef itkParabolicErodeDilateImageFilter_h +#define itkParabolicErodeDilateImageFilter_h + +#include "itkImageToImageFilter.h" +#include "itkNumericTraits.h" +#include "itkProgressReporter.h" + +namespace itk +{ +/** + * \class ParabolicErodeDilateImageFilter + * \brief Parent class for morphological operations with parabolic + * structuring elements. + * + * Parabolic structuring elements are the morphological counterpart of + * the gaussian function in linear processing. Parabolic structuring + * elements are dimensionally decomposable and fast algorithms are + * available for computing erosions and dilations along lines. + * This class implements the "point of contact" algorithm and the + * "intersection" algorithm. The contact point algorithm is faster at + * small scales. The intersection algorithm was rediscovered by + * Felzenszwalb & Huttenlocher, but was actually described and tested + * by van den Boomgaard many years earlier. The intersection algorithm + * is faster for scales > 1, and independent of scale. It also seems + * to perform significantly better in the distance transform + * application. The intersection algorithm is also used in the IJ + * article on generalised distance transforms. + * + * Parabolic structuring functions can be used as a fast alternative + * to the "rolling ball" structuring element classically used in + * background estimation, for example in ImageJ, have applications in + * image sharpening and distance transform computation. + * + * This class uses an internal buffer of RealType pixels for each + * line. This line is cast to the output pixel type when written back + * to the output image. Since the filter uses dimensional + * decomposition this approach could result in inaccuracy as pixels + * are cast back and forth between low and high precision types. Use a + * high precision output type and cast manually if this is a problem. + * + * Boomgaard, R. van den and Dorst, L. and Makram-Ebeid, L.S. and + * Schavemaker, J. Quadratic structuring functions in mathematical + * morphology. Mathematical Morphology and its Applications to Image + * and Signal Processing. + * + * Felzenszwalb, P.F. & Huttenlocher, D.P. Distance Transforms of Sampled Functions. + * Techreport: Cornell Computing and Information Science, 2004. + * + * This filter is threaded. Threading mechanism derived from + * SignedMaurerDistanceMap extensions by Gaetan Lehman + * + * Core methods described in the InsightJournal article: + * "Morphology with parabolic structuring elements" + * + * https://hdl.handle.net/1926/1370 + * + * + * \ingroup ParabolicMorphology + * + * \author Richard Beare, Department of Medicine, Monash University, + * Australia. + * + **/ + +template +class ITK_TEMPLATE_EXPORT ParabolicErodeDilateImageFilter : public ImageToImageFilter +{ +public: + ITK_DISALLOW_COPY_AND_MOVE(ParabolicErodeDilateImageFilter); + + /** Standard class type alias. */ + using Self = ParabolicErodeDilateImageFilter; + using Superclass = ImageToImageFilter; + using Pointer = SmartPointer; + using ConstPointer = SmartPointer; + + /** Method for creation through the object factory. */ + itkNewMacro(Self); + + /** Runtime information support. */ + itkTypeMacro(ParabolicErodeDilateImageFilter, ImageToImageFilter); + + /** Pixel Type of the input image */ + using InputImageType = TInputImage; + using OutputImageType = TOutputImage; + using PixelType = typename TInputImage::PixelType; + using RealType = typename NumericTraits::RealType; + using ScalarRealType = typename NumericTraits::ScalarRealType; + using OutputPixelType = typename TOutputImage::PixelType; + + /** Smart pointer type alias support. */ + using InputImagePointer = typename TInputImage::Pointer; + using InputImageConstPointer = typename TInputImage::ConstPointer; + using InputSizeType = typename TInputImage::SizeType; + using OutputSizeType = typename TOutputImage::SizeType; + + using OutputIndexType = typename OutputImageType::IndexType; + + /** a type to represent the "kernel radius" */ + using RadiusType = typename itk::FixedArray; + + /** Image dimension. */ + static constexpr unsigned int ImageDimension = TInputImage::ImageDimension; + static constexpr unsigned int OutputImageDimension = TOutputImage::ImageDimension; + static constexpr unsigned int InputImageDimension = TInputImage::ImageDimension; + + using OutputImageRegionType = typename OutputImageType::RegionType; + /** Define the image type for internal computations + RealType is usually 'double' in NumericTraits. + Here we prefer float in order to save memory. */ + + using InternalRealType = typename NumericTraits::FloatType; + // using RealImageType = typename Image; + + // set all of the scales the same + void + SetScale(ScalarRealType scale); + + itkSetMacro(Scale, RadiusType); + itkGetConstReferenceMacro(Scale, RadiusType); + + enum ParabolicAlgorithm + { + NOCHOICE = 0, // decices based on scale - experimental + CONTACTPOINT = 1, // sometimes faster at low scale + INTERSECTION = 2 // default + }; + /** + * Set/Get the method used. Choices are contact point or + * intersection. Intersection is the default. Contact point can be + * faster at small scales. + */ + + itkSetMacro(ParabolicAlgorithm, int); + itkGetConstReferenceMacro(ParabolicAlgorithm, int); + + /** + * Set/Get whether the scale refers to pixels or world units - + * default is false + */ + itkSetMacro(UseImageSpacing, bool); + itkGetConstReferenceMacro(UseImageSpacing, bool); + itkBooleanMacro(UseImageSpacing); + /** Image related type alias. */ + +#ifdef ITK_USE_CONCEPT_CHECKING + /** Begin concept checking */ + itkConceptMacro(SameDimension, + (Concept::SameDimension)); + + itkConceptMacro(Comparable, (Concept::Comparable)); + + /** End concept checking */ +#endif +protected: + ParabolicErodeDilateImageFilter(); + ~ParabolicErodeDilateImageFilter() override = default; + void + PrintSelf(std::ostream & os, Indent indent) const override; + + /** Generate Data */ + void + GenerateData() override; + + unsigned int + SplitRequestedRegion(unsigned int i, unsigned int num, OutputImageRegionType & splitRegion) override; + + void + ThreadedGenerateData(const OutputImageRegionType & outputRegionForThread, ThreadIdType threadId) override; + + void + GenerateInputRequestedRegion() override; + + // Override since the filter produces the entire dataset. + void + EnlargeOutputRequestedRegion(DataObject * output) override; + + bool m_UseImageSpacing; + int m_ParabolicAlgorithm; + +private: + RadiusType m_Scale; + + int m_CurrentDimension; +}; +} // end namespace itk + +#ifndef ITK_MANUAL_INSTANTIATION +# include "itkParabolicErodeDilateImageFilter.hxx" +#endif + +#endif diff --git a/Modules/Filtering/ParabolicMorphology/include/itkParabolicErodeDilateImageFilter.hxx b/Modules/Filtering/ParabolicMorphology/include/itkParabolicErodeDilateImageFilter.hxx new file mode 100644 index 00000000000..77b49172ae2 --- /dev/null +++ b/Modules/Filtering/ParabolicMorphology/include/itkParabolicErodeDilateImageFilter.hxx @@ -0,0 +1,298 @@ +/*========================================================================= + * + * Copyright NumFOCUS + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * https://www.apache.org/licenses/LICENSE-2.0.txt + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + *=========================================================================*/ +#ifndef itkParabolicErodeDilateImageFilter_hxx +#define itkParabolicErodeDilateImageFilter_hxx + +#include "itkImageRegionConstIterator.h" +#include "itkImageRegionIterator.h" + +// #define NOINDEX +#ifndef NOINDEX +# include "itkImageLinearIteratorWithIndex.h" +# include "itkImageLinearConstIteratorWithIndex.h" +#else +# include "itkImageLinearIterator.h" +# include "itkImageLinearConstIterator.h" +#endif +#include "itkParabolicMorphUtils.h" + +namespace itk +{ +template +ParabolicErodeDilateImageFilter::ParabolicErodeDilateImageFilter() +{ + this->SetNumberOfRequiredOutputs(1); + this->SetNumberOfRequiredInputs(1); + + m_UseImageSpacing = false; + m_ParabolicAlgorithm = INTERSECTION; + + this->DynamicMultiThreadingOff(); +} + +template +unsigned int +ParabolicErodeDilateImageFilter::SplitRequestedRegion( + unsigned int i, + unsigned int num, + OutputImageRegionType & splitRegion) +{ + // Get the output pointer + OutputImageType * outputPtr = this->GetOutput(); + + // Initialize the splitRegion to the output requested region + splitRegion = outputPtr->GetRequestedRegion(); + + const OutputSizeType & requestedRegionSize = splitRegion.GetSize(); + + OutputIndexType splitIndex = splitRegion.GetIndex(); + OutputSizeType splitSize = splitRegion.GetSize(); + + // split on the outermost dimension available + // and avoid the current dimension + int splitAxis = static_cast(outputPtr->GetImageDimension()) - 1; + while ((requestedRegionSize[splitAxis] == 1) || (splitAxis == static_cast(m_CurrentDimension))) + { + --splitAxis; + if (splitAxis < 0) + { // cannot split + itkDebugMacro("Cannot Split"); + return 1; + } + } + + // determine the actual number of pieces that will be generated + auto range = static_cast(requestedRegionSize[splitAxis]); + + auto valuesPerThread = static_cast(std::ceil(range / static_cast(num))); + unsigned int maxThreadIdUsed = static_cast(std::ceil(range / static_cast(valuesPerThread))) - 1; + + // Split the region + if (i < maxThreadIdUsed) + { + splitIndex[splitAxis] += i * valuesPerThread; + splitSize[splitAxis] = valuesPerThread; + } + if (i == maxThreadIdUsed) + { + splitIndex[splitAxis] += i * valuesPerThread; + // last thread needs to process the "rest" dimension being split + splitSize[splitAxis] = splitSize[splitAxis] - i * valuesPerThread; + } + + // set the split region ivars + splitRegion.SetIndex(splitIndex); + splitRegion.SetSize(splitSize); + + itkDebugMacro("Split Piece: " << splitRegion); + + return maxThreadIdUsed + 1; +} + +template +void +ParabolicErodeDilateImageFilter::SetScale(ScalarRealType scale) +{ + RadiusType s; + + s.Fill(scale); + this->SetScale(s); +} + +#if 1 +template +void +ParabolicErodeDilateImageFilter::GenerateInputRequestedRegion() +{ + // call the superclass' implementation of this method. this should + // copy the output requested region to the input requested region + Superclass::GenerateInputRequestedRegion(); + + // This filter needs all of the input + InputImagePointer image = const_cast(this->GetInput()); + if (image) + { + image->SetRequestedRegion(this->GetInput()->GetLargestPossibleRegion()); + } +} + +#endif +#if 1 +template +void +ParabolicErodeDilateImageFilter::EnlargeOutputRequestedRegion(DataObject * output) +{ + auto * out = dynamic_cast(output); + + if (out) + { + out->SetRequestedRegion(out->GetLargestPossibleRegion()); + } +} + +#endif + +template +void +ParabolicErodeDilateImageFilter::GenerateData() +{ + ThreadIdType nbthreads = this->GetNumberOfWorkUnits(); + + typename TInputImage::ConstPointer inputImage(this->GetInput()); + typename TOutputImage::Pointer outputImage(this->GetOutput()); + + // const unsigned int imageDimension = inputImage->GetImageDimension(); + outputImage->SetBufferedRegion(outputImage->GetRequestedRegion()); + outputImage->Allocate(); + + // Set up the multithreaded processing + typename ImageSource::ThreadStruct str; + str.Filter = this; + + itk::MultiThreaderBase * multithreader = this->GetMultiThreader(); + multithreader->SetNumberOfWorkUnits(nbthreads); + multithreader->SetSingleMethod(this->ThreaderCallback, &str); + + // multithread the execution + for (unsigned int d = 0; d < ImageDimension; d++) + { + m_CurrentDimension = d; + multithreader->SingleMethodExecute(); + } +} + +template +void +ParabolicErodeDilateImageFilter::ThreadedGenerateData( + const OutputImageRegionType & outputRegionForThread, + ThreadIdType threadId) +{ + // compute the number of rows first, so we can setup a progress reporter + unsigned int numberOfRows = 1; + InputSizeType size = outputRegionForThread.GetSize(); + + for (unsigned int d = 0; d < InputImageDimension; ++d) + { + if (d != m_CurrentDimension) + { + numberOfRows *= size[d]; + } + } + + float progressPerDimension = 1.0 / ImageDimension; + + ProgressReporter progress( + this, threadId, numberOfRows, 30, m_CurrentDimension * progressPerDimension, progressPerDimension); + + using InputConstIteratorType = ImageLinearConstIteratorWithIndex; + using OutputIteratorType = ImageLinearIteratorWithIndex; + + // for stages after the first + using OutputConstIteratorType = ImageLinearConstIteratorWithIndex; + + using RegionType = ImageRegion; + + typename TInputImage::ConstPointer inputImage(this->GetInput()); + typename TOutputImage::Pointer outputImage(this->GetOutput()); + + // outputImage->SetBufferedRegion( outputImage->GetRequestedRegion() ); + // outputImage->Allocate(); + RegionType region = outputRegionForThread; + + InputConstIteratorType inputIterator(inputImage, region); + OutputIteratorType outputIterator(outputImage, region); + OutputConstIteratorType inputIteratorStage2(outputImage, region); + + // deal with the first dimension - this should be copied to the + // output if the scale is 0 + if (m_CurrentDimension == 0) + { + if (m_Scale[0] > 0) + { + // Perform as normal + // RealType magnitude = 1.0/(2.0 * m_Scale[0]); + unsigned long LineLength = region.GetSize()[0]; + RealType image_scale = this->GetInput()->GetSpacing()[0]; + + doOneDimension( + inputIterator, + outputIterator, + progress, + LineLength, + 0, + this->m_UseImageSpacing, + image_scale, + this->m_Scale[0], + m_ParabolicAlgorithm); + } + else + { + // copy to output + using InItType = ImageRegionConstIterator; + using OutItType = ImageRegionIterator; + + InItType InIt(inputImage, region); + OutItType OutIt(outputImage, region); + while (!InIt.IsAtEnd()) + { + OutIt.Set(static_cast(InIt.Get())); + ++InIt; + ++OutIt; + } + } + } + else + { + // other dimensions + if (m_Scale[m_CurrentDimension] > 0) + { + // create a vector to buffer lines + unsigned long LineLength = region.GetSize()[m_CurrentDimension]; + // RealType magnitude = 1.0/(2.0 * m_Scale[dd]); + RealType image_scale = this->GetInput()->GetSpacing()[m_CurrentDimension]; + + doOneDimension( + inputIteratorStage2, + outputIterator, + progress, + LineLength, + m_CurrentDimension, + this->m_UseImageSpacing, + image_scale, + this->m_Scale[m_CurrentDimension], + m_ParabolicAlgorithm); + } + } +} + +template +void +ParabolicErodeDilateImageFilter::PrintSelf(std::ostream & os, Indent indent) const +{ + Superclass::PrintSelf(os, indent); + if (m_UseImageSpacing) + { + os << "Scale in world units: " << m_Scale << std::endl; + } + else + { + os << "Scale in voxels: " << m_Scale << std::endl; + } +} +} // namespace itk +#endif diff --git a/Modules/Filtering/ParabolicMorphology/include/itkParabolicErodeImageFilter.h b/Modules/Filtering/ParabolicMorphology/include/itkParabolicErodeImageFilter.h new file mode 100644 index 00000000000..a1c872ea146 --- /dev/null +++ b/Modules/Filtering/ParabolicMorphology/include/itkParabolicErodeImageFilter.h @@ -0,0 +1,94 @@ +/*========================================================================= + * + * Copyright NumFOCUS + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * https://www.apache.org/licenses/LICENSE-2.0.txt + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + *=========================================================================*/ +#ifndef itkParabolicErodeImageFilter_h +#define itkParabolicErodeImageFilter_h + +#include "itkParabolicErodeDilateImageFilter.h" +#include "itkNumericTraits.h" + +namespace itk +{ +/** + * \class ParabolicErodeImageFilter + * \brief Class for morphological erosion + * operations with parabolic structuring elements. + * + * This filter doesn't use the erode/dilate classes directly so + * that multiple image copies aren't necessary. + * + * This filter is threaded. + * Core methods described in the InsightJournal article: + * "Morphology with parabolic structuring elements" + * + * https://hdl.handle.net/1926/1370 + * + * \sa itkParabolicErodeDilateImageFilter + * + * + * \ingroup ParabolicMorphology + * + * \author Richard Beare, Department of Medicine, Monash University, + * Australia. + **/ +template +class ITK_TEMPLATE_EXPORT ParabolicErodeImageFilter + : public ParabolicErodeDilateImageFilter +{ +public: + ITK_DISALLOW_COPY_AND_MOVE(ParabolicErodeImageFilter); + + /** Standard class type alias. */ + using Self = ParabolicErodeImageFilter; + using Superclass = ParabolicErodeDilateImageFilter; + using Pointer = SmartPointer; + using ConstPointer = SmartPointer; + + /** Method for creation through the object factory. */ + itkNewMacro(Self); + + /** Runtime information support. */ + itkTypeMacro(ParabolicErodeImageFilter, ParabolicErodeDilateImageFilter); + + /** Pixel Type of the input image */ + using InputImageType = TInputImage; + using OutputImageType = TOutputImage; + using PixelType = typename TInputImage::PixelType; + using RealType = typename NumericTraits::RealType; + using ScalarRealType = typename NumericTraits::ScalarRealType; + using OutputPixelType = typename TOutputImage::PixelType; + + /** Smart pointer type alias support. */ + using InputImagePointer = typename TInputImage::Pointer; + using InputImageConstPointer = typename TInputImage::ConstPointer; + + /** a type to represent the "kernel radius" */ + using RadiusType = typename itk::FixedArray; + + /** Image dimension. */ + static constexpr unsigned int ImageDimension = TInputImage::ImageDimension; + /** Define the image type for internal computations + RealType is usually 'double' in NumericTraits. + Here we prefer float in order to save memory. */ +protected: + ParabolicErodeImageFilter() = default; + ~ParabolicErodeImageFilter() override = default; + // void PrintSelf(std::ostream& os, Indent indent) const; +}; +} // end namespace itk + +#endif diff --git a/Modules/Filtering/ParabolicMorphology/include/itkParabolicMorphUtils.h b/Modules/Filtering/ParabolicMorphology/include/itkParabolicMorphUtils.h new file mode 100644 index 00000000000..157a2248d50 --- /dev/null +++ b/Modules/Filtering/ParabolicMorphology/include/itkParabolicMorphUtils.h @@ -0,0 +1,320 @@ +/*========================================================================= + * + * Copyright NumFOCUS + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * https://www.apache.org/licenses/LICENSE-2.0.txt + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + *=========================================================================*/ +#ifndef itkParabolicMorphUtils_h +#define itkParabolicMorphUtils_h + +#include + +#include "itkProgressReporter.h" + +namespace itk +{ +// contact point algorithm +template +void +DoLineCP(LineBufferType & LineBuf, LineBufferType & tmpLineBuf, const RealType magnitude) +{ + static constexpr RealType extreme = + doDilate ? NumericTraits::NonpositiveMin() : NumericTraits::max(); + + // contact point algorithm + long koffset = 0, newcontact = 0; // how far away the search starts. + + const long LineLength = LineBuf.size(); + + // negative half of the parabola + for (long pos = 0; pos < LineLength; pos++) + { + auto BaseVal = extreme; // the base value for comparison + for (long krange = koffset; krange <= 0; krange++) + { + // difference needs to be paramaterised + RealType T = LineBuf[pos + krange] - magnitude * krange * krange; + // switch on template parameter - hopefully gets optimized away. + if (doDilate ? (T >= BaseVal) : (T <= BaseVal)) + { + BaseVal = T; + newcontact = krange; + } + } + tmpLineBuf[pos] = BaseVal; + koffset = newcontact - 1; + } + // positive half of parabola + koffset = newcontact = 0; + for (long pos = LineLength - 1; pos >= 0; pos--) + { + auto BaseVal = extreme; // the base value for comparison + for (long krange = koffset; krange >= 0; krange--) + { + RealType T = tmpLineBuf[pos + krange] - magnitude * krange * krange; + if (doDilate ? (T >= BaseVal) : (T <= BaseVal)) + { + BaseVal = T; + newcontact = krange; + } + } + LineBuf[pos] = BaseVal; + koffset = newcontact + 1; + } +} + +// intersection algorithm +// This algorithm has been described a couple of times. First by van +// den Boomgaard and more recently by Felzenszwalb and Huttenlocher, +// in the context of generalized distance transform +template +void +DoLineIntAlg(LineBufferType & LineBuf, + EnvBufferType & F, + IndexBufferType & v, + EnvBufferType & z, + const RealType magnitude) +{ + int k; /* Index of rightmost parabola in lower envelope */ + /* Locations of parabolas in lower envelope */ + /* these need to be int, rather than unsigned, to make boundary + conditions easy to test */ + + // v stores locations of parabolas in the lower envelope + // z stores thr location of boundaries between parabolas + + // I've gone nuts with the static casts etc, because I seemed to + // have strange behaviour when I didn't do this. Also managed to get + // rid of all the warnings by sticking to size_t and equivalents. + RealType s; + + /* holds precomputed scale*f(q) + q^2 for speedup */ + // LineBufferType F(LineBuf.size()); + + // initialize + k = 0; + v[0] = 0; + z[0] = NumericTraits::NonpositiveMin(); + z[1] = NumericTraits::max(); + F[0] = LineBuf[0] / magnitude; + const size_t N(LineBuf.size()); + + for (size_t q = 1; q < N; q++) /* main loop */ + { + if (doDilate) + { + /* precompute f(q) + q^2 for speedup */ + F[q] = (LineBuf[q] / magnitude) - (static_cast(q) * static_cast(q)); + k++; + do + { + /* remove last parabola from surface */ + k--; + /* compute intersection */ + s = (F[q] - F[v[k]]) / (2.0 * (v[k] - static_cast(q))); + } while (s <= z[k]); + /* bump k to add new parabola */ + k++; + } + else + { + /* precompute f(q) + q^2 for speedup */ + F[q] = (LineBuf[q] / magnitude) + (static_cast(q) * static_cast(q)); + k++; + do + { + /* remove last parabola from surface */ + k--; + /* compute intersection */ + s = (F[q] - F[v[k]]) / (2.0 * (static_cast(q) - v[k])); + } while (s <= z[k]); + /* bump k to add new parabola */ + k++; + } + v[k] = q; + z[k] = s; + itkAssertInDebugAndIgnoreInReleaseMacro((size_t)(k + 1) <= N); + z[k + 1] = NumericTraits::max(); + } /* for q */ + /* now reconstruct output */ + if (doDilate) + { + k = 0; + for (size_t q = 0; q < N; q++) + { + while (z[k + 1] < static_cast(q)) + { + k++; + } + itkAssertInDebugAndIgnoreInReleaseMacro(static_cast(v[k]) < N); + itkAssertInDebugAndIgnoreInReleaseMacro(static_cast(v[k]) >= 0); + LineBuf[q] = static_cast( + (F[v[k]] - (static_cast(q) * (static_cast(q) - 2 * v[k]))) * magnitude); + } + } + else + { + k = 0; + for (size_t q = 0; q < N; q++) + { + while (z[k + 1] < static_cast(q)) + { + k++; + } + itkAssertInDebugAndIgnoreInReleaseMacro(static_cast(v[k]) < N); + itkAssertInDebugAndIgnoreInReleaseMacro(static_cast(v[k]) >= 0); + LineBuf[q] = ((static_cast(q) * (static_cast(q) - 2 * v[k]) + F[v[k]]) * magnitude); + } + } +} + +template +void +doOneDimension(TInIter & inputIterator, + TOutIter & outputIterator, + ProgressReporter & progress, + const long LineLength, + const unsigned direction, + const bool m_UseImageSpacing, + const RealType image_scale, + const RealType Sigma, + int ParabolicAlgorithmChoice) +{ + enum ParabolicAlgorithm + { + NOCHOICE = 0, // decices based on scale - experimental + CONTACTPOINT = 1, // sometimes faster at low scale + INTERSECTION = 2 // default + }; + + // using LineBufferType = typename std::vector; + + // message from M.Starring suggested performance gain using Array + // instead of std::vector. + using LineBufferType = typename itk::Array; + RealType iscale = 1.0; + if (m_UseImageSpacing) + { + iscale = image_scale; + } + if (ParabolicAlgorithmChoice == NOCHOICE) + { + // both set to true or false - use scale to figure it out + if ((2.0 * Sigma) < 0.2) + { + ParabolicAlgorithmChoice = CONTACTPOINT; + } + else + { + ParabolicAlgorithmChoice = INTERSECTION; + } + } + + if (ParabolicAlgorithmChoice == CONTACTPOINT) + { + // using the contact point algorithm + + // const RealType magnitude = magnitudeSign * 1.0/(2.0 * + // Sigma/(iscale*iscale)); + // restructure equation to reduce numerical error + constexpr int magnitudeSign = doDilate ? 1 : -1; + const RealType magnitudeCP = (magnitudeSign * iscale * iscale) / (2.0 * Sigma); + + LineBufferType LineBuf(LineLength); + LineBufferType tmpLineBuf(LineLength); + inputIterator.SetDirection(direction); + outputIterator.SetDirection(direction); + inputIterator.GoToBegin(); + outputIterator.GoToBegin(); + + while (!inputIterator.IsAtEnd() && !outputIterator.IsAtEnd()) + { + // process this direction + // fetch the line into the buffer - this methodology is like + // the gaussian filters + + unsigned int i = 0; + while (!inputIterator.IsAtEndOfLine()) + { + LineBuf[i++] = static_cast(inputIterator.Get()); + ++inputIterator; + } + + DoLineCP(LineBuf, tmpLineBuf, magnitudeCP); + // copy the line back + unsigned int j = 0; + while (!outputIterator.IsAtEndOfLine()) + { + outputIterator.Set(static_cast(LineBuf[j++])); + ++outputIterator; + } + + // now onto the next line + inputIterator.NextLine(); + outputIterator.NextLine(); + progress.CompletedPixel(); + } + } + else + { + // using the Intersection algorithm + using IndexBufferType = typename itk::Array; + + const RealType magnitudeInt = (iscale * iscale) / (2.0 * Sigma); + LineBufferType LineBuf(LineLength); + LineBufferType Fbuf(LineLength); + IndexBufferType Vbuf(LineLength); + LineBufferType Zbuf(LineLength + 1); + + inputIterator.SetDirection(direction); + outputIterator.SetDirection(direction); + inputIterator.GoToBegin(); + outputIterator.GoToBegin(); + + while (!inputIterator.IsAtEnd() && !outputIterator.IsAtEnd()) + { + // process this direction + // fetch the line into the buffer - this methodology is like + // the gaussian filters + + unsigned int i = 0; + while (!inputIterator.IsAtEndOfLine()) + { + LineBuf[i++] = static_cast(inputIterator.Get()); + ++inputIterator; + } + DoLineIntAlg( + LineBuf, Fbuf, Vbuf, Zbuf, magnitudeInt); + // copy the line back + unsigned int j = 0; + while (!outputIterator.IsAtEndOfLine()) + { + outputIterator.Set(static_cast(LineBuf[j++])); + ++outputIterator; + } + + // now onto the next line + inputIterator.NextLine(); + outputIterator.NextLine(); + progress.CompletedPixel(); + } + } +} +} // namespace itk +#endif diff --git a/Modules/Filtering/ParabolicMorphology/include/itkParabolicOpenCloseImageFilter.h b/Modules/Filtering/ParabolicMorphology/include/itkParabolicOpenCloseImageFilter.h new file mode 100644 index 00000000000..944cde9a56f --- /dev/null +++ b/Modules/Filtering/ParabolicMorphology/include/itkParabolicOpenCloseImageFilter.h @@ -0,0 +1,180 @@ +/*========================================================================= + * + * Copyright NumFOCUS + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * https://www.apache.org/licenses/LICENSE-2.0.txt + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + *=========================================================================*/ +#ifndef itkParabolicOpenCloseImageFilter_h +#define itkParabolicOpenCloseImageFilter_h + +#include "itkImageToImageFilter.h" +#include "itkNumericTraits.h" +#include "itkProgressReporter.h" + +namespace itk +{ +/** + * \class ParabolicOpenCloseImageFilter + * \brief Parent class for morphological opening and closing + * operations with parabolic structuring elements. + * + * This filter doesn't use the erode/dilate classes directly so + * that multiple image copies aren't necessary. + * + * This filter is threaded. Threading mechanism derived from + * SignedMaurerDistanceMap extensions by Gaetan Lehman + * + * Core methods described in the InsightJournal article: + * "Morphology with parabolic structuring elements" + * + * https://hdl.handle.net/1926/1370 + * + * \sa itkParabolicErodeDilateImageFilter + * + * + * \ingroup ParabolicMorphology + * + * \author Richard Beare, Department of Medicine, Monash University, + * Australia. + * + **/ +template +class ITK_TEMPLATE_EXPORT ParabolicOpenCloseImageFilter : public ImageToImageFilter +{ +public: + ITK_DISALLOW_COPY_AND_MOVE(ParabolicOpenCloseImageFilter); + + /** Standard class type alias. */ + using Self = ParabolicOpenCloseImageFilter; + using Superclass = ImageToImageFilter; + using Pointer = SmartPointer; + using ConstPointer = SmartPointer; + + /** Method for creation through the object factory. */ + itkNewMacro(Self); + + /** Runtime information support. */ + itkTypeMacro(ParabolicOpenCloseImageFilter, ImageToImageFilter); + + /** Pixel Type of the input image */ + using InputImageType = TInputImage; + using OutputImageType = TOutputImage; + using PixelType = typename TInputImage::PixelType; + using RealType = typename NumericTraits::RealType; + using ScalarRealType = typename NumericTraits::ScalarRealType; + using OutputPixelType = typename TOutputImage::PixelType; + using OutputImageRegionType = typename OutputImageType::RegionType; + using InputSizeType = typename TInputImage::SizeType; + using OutputSizeType = typename TOutputImage::SizeType; + + /** Smart pointer type alias support. */ + using InputImagePointer = typename TInputImage::Pointer; + using InputImageConstPointer = typename TInputImage::ConstPointer; + + using OutputIndexType = typename OutputImageType::IndexType; + + /** Image related type alias. */ + static constexpr unsigned int OutputImageDimension = TOutputImage::ImageDimension; + static constexpr unsigned int InputImageDimension = TInputImage::ImageDimension; + static constexpr unsigned int ImageDimension = TInputImage::ImageDimension; + + /** a type to represent the "kernel radius" */ + using RadiusType = typename itk::FixedArray; + + /** Define the image type for internal computations + RealType is usually 'double' in NumericTraits. + Here we prefer float in order to save memory. */ + + using InternalRealType = typename NumericTraits::FloatType; + // using RealImageType = typename Image; + + // set all of the scales the same + void + SetScale(ScalarRealType scale); + + itkSetMacro(Scale, RadiusType); + itkGetConstReferenceMacro(Scale, RadiusType); + /** + * Set/Get whether the scale refers to pixels or world units - + * default is false + */ + itkSetMacro(UseImageSpacing, bool); + itkGetConstReferenceMacro(UseImageSpacing, bool); + itkBooleanMacro(UseImageSpacing); + + enum ParabolicAlgorithm + { + NOCHOICE = 0, // decices based on scale - experimental + CONTACTPOINT = 1, // sometimes faster at low scale + INTERSECTION = 2 // default + }; + /** + * Set/Get the method used. Choices are contact point or + * intersection. Intersection is the default. Contact point can be + * faster at small scales. + */ + + itkSetMacro(ParabolicAlgorithm, int); + itkGetConstReferenceMacro(ParabolicAlgorithm, int); + +#ifdef ITK_USE_CONCEPT_CHECKING + /** Begin concept checking */ + itkConceptMacro(SameDimension, + (Concept::SameDimension)); + + itkConceptMacro(Comparable, (Concept::Comparable)); + + /** End concept checking */ +#endif +protected: + ParabolicOpenCloseImageFilter(); + ~ParabolicOpenCloseImageFilter() override = default; + void + PrintSelf(std::ostream & os, Indent indent) const override; + + /** Generate Data */ + void + GenerateData() override; + + unsigned int + SplitRequestedRegion(unsigned int i, unsigned int num, OutputImageRegionType & splitRegion) override; + + void + ThreadedGenerateData(const OutputImageRegionType & outputRegionForThread, ThreadIdType threadId) override; + + void + GenerateInputRequestedRegion() override; + + // Override since the filter produces the entire dataset. + void + EnlargeOutputRequestedRegion(DataObject * output) override; + + int m_ParabolicAlgorithm; + +private: + RadiusType m_Scale; + + int m_CurrentDimension; + int m_Stage; + bool m_UseImageSpacing; +}; +} // end namespace itk + +#ifndef ITK_MANUAL_INSTANTIATION +# include "itkParabolicOpenCloseImageFilter.hxx" +#endif + +#endif diff --git a/Modules/Filtering/ParabolicMorphology/include/itkParabolicOpenCloseImageFilter.hxx b/Modules/Filtering/ParabolicMorphology/include/itkParabolicOpenCloseImageFilter.hxx new file mode 100644 index 00000000000..10d2598a7a2 --- /dev/null +++ b/Modules/Filtering/ParabolicMorphology/include/itkParabolicOpenCloseImageFilter.hxx @@ -0,0 +1,381 @@ +/*========================================================================= + * + * Copyright NumFOCUS + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * https://www.apache.org/licenses/LICENSE-2.0.txt + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + *=========================================================================*/ +#ifndef itkParabolicOpenCloseImageFilter_hxx +#define itkParabolicOpenCloseImageFilter_hxx + +#include "itkImageRegionConstIterator.h" +#include "itkImageRegionIterator.h" + +// #define NOINDEX +#ifndef NOINDEX +# include "itkImageLinearIteratorWithIndex.h" +# include "itkImageLinearConstIteratorWithIndex.h" +#else +# include "itkImageLinearIterator.h" +# include "itkImageLinearConstIterator.h" +#endif +#include "itkStatisticsImageFilter.h" +#include "itkParabolicMorphUtils.h" + +namespace itk +{ +template +ParabolicOpenCloseImageFilter::ParabolicOpenCloseImageFilter() +{ + this->SetNumberOfRequiredOutputs(1); + this->SetNumberOfRequiredInputs(1); + // needs to be selected according to erosion/dilation + m_UseImageSpacing = false; + m_ParabolicAlgorithm = INTERSECTION; + m_Stage = 1; // indicate whether we are on the first pass or the + // second + + this->DynamicMultiThreadingOff(); +} + +template +unsigned int +ParabolicOpenCloseImageFilter::SplitRequestedRegion( + unsigned int i, + unsigned int num, + OutputImageRegionType & splitRegion) +{ + // Get the output pointer + OutputImageType * outputPtr = this->GetOutput(); + + // Initialize the splitRegion to the output requested region + splitRegion = outputPtr->GetRequestedRegion(); + + const OutputSizeType & requestedRegionSize = splitRegion.GetSize(); + + OutputIndexType splitIndex = splitRegion.GetIndex(); + OutputSizeType splitSize = splitRegion.GetSize(); + + // split on the outermost dimension available + // and avoid the current dimension + int splitAxis = static_cast(outputPtr->GetImageDimension()) - 1; + while ((requestedRegionSize[splitAxis] == 1) || (splitAxis == static_cast(m_CurrentDimension))) + { + --splitAxis; + if (splitAxis < 0) + { // cannot split + itkDebugMacro("Cannot Split"); + return 1; + } + } + + // determine the actual number of pieces that will be generated + auto range = static_cast(requestedRegionSize[splitAxis]); + + auto valuesPerThread = static_cast(std::ceil(range / static_cast(num))); + unsigned int maxThreadIdUsed = static_cast(std::ceil(range / static_cast(valuesPerThread))) - 1; + + // Split the region + if (i < maxThreadIdUsed) + { + splitIndex[splitAxis] += i * valuesPerThread; + splitSize[splitAxis] = valuesPerThread; + } + if (i == maxThreadIdUsed) + { + splitIndex[splitAxis] += i * valuesPerThread; + // last thread needs to process the "rest" dimension being split + splitSize[splitAxis] = splitSize[splitAxis] - i * valuesPerThread; + } + + // set the split region ivars + splitRegion.SetIndex(splitIndex); + splitRegion.SetSize(splitSize); + + itkDebugMacro("Split Piece: " << splitRegion); + + return maxThreadIdUsed + 1; +} + +template +void +ParabolicOpenCloseImageFilter::SetScale(ScalarRealType scale) +{ + RadiusType s; + + s.Fill(scale); + this->SetScale(s); +} + +#if 1 +template +void +ParabolicOpenCloseImageFilter::GenerateInputRequestedRegion() +{ + // call the superclass' implementation of this method. this should + // copy the output requested region to the input requested region + Superclass::GenerateInputRequestedRegion(); + + // This filter needs all of the input + InputImagePointer image = const_cast(this->GetInput()); + if (image) + { + image->SetRequestedRegion(this->GetInput()->GetLargestPossibleRegion()); + } +} + +#endif + +#if 1 +template +void +ParabolicOpenCloseImageFilter::EnlargeOutputRequestedRegion(DataObject * output) +{ + auto * out = dynamic_cast(output); + + if (out) + { + out->SetRequestedRegion(out->GetLargestPossibleRegion()); + } +} + +#endif + +template +void +ParabolicOpenCloseImageFilter::GenerateData() +{ + ThreadIdType nbthreads = this->GetNumberOfWorkUnits(); + + // using InputConstIteratorType = ImageLinearConstIteratorWithIndex< TInputImage > ; + // using OutputIteratorType = ImageLinearIteratorWithIndex< TOutputImage >; + + // for stages after the first + // using OutputConstIteratorType = ImageLinearConstIteratorWithIndex< TOutputImage > ; + + // using RegionType = ImageRegion< TInputImage::ImageDimension >; + + typename TInputImage::ConstPointer inputImage(this->GetInput()); + typename TOutputImage::Pointer outputImage(this->GetOutput()); + + // const unsigned int imageDimension = inputImage->GetImageDimension(); + + outputImage->SetBufferedRegion(outputImage->GetRequestedRegion()); + outputImage->Allocate(); + + typename ImageSource::ThreadStruct str; + str.Filter = this; + + itk::MultiThreaderBase * multithreader = this->GetMultiThreader(); + multithreader->SetNumberOfWorkUnits(nbthreads); + multithreader->SetSingleMethod(this->ThreaderCallback, &str); + + // multithread the execution + + // multithread the execution - stage 1 + m_Stage = 1; + + for (unsigned int d = 0; d < ImageDimension; d++) + { + m_CurrentDimension = d; + multithreader->SingleMethodExecute(); + } + + // multithread the execution - stage 2 + m_Stage = 2; + for (unsigned int d = 0; d < ImageDimension; d++) + { + m_CurrentDimension = d; + multithreader->SingleMethodExecute(); + } + + m_Stage = 1; + +#if 0 + // Set up the multithreaded processing + typename ImageSource< TOutputImage >::ThreadStruct str; + str.Filter = this; + this->GetMultiThreader()->SetNumberOfWorkUnits( this->GetNumberOfWorkUnits() ); + this->GetMultiThreader()->SetSingleMethod(this->ThreaderCallback, &str); + + // multithread the execution - stage 1 + m_Stage = 1; + for ( unsigned int d = 0; d < ImageDimension; d++ ) + { + m_CurrentDimension = d; + this->GetMultiThreader()->SingleMethodExecute(); + } + // swap over the parameters controlling erosion/dilation + m_Extreme = m_Extreme2; + m_MagnitudeSign = m_MagnitudeSign2; + + // multithread the execution - stage 2 + m_Stage = 2; + for ( unsigned int d = 0; d < ImageDimension; d++ ) + { + m_CurrentDimension = d; + this->GetMultiThreader()->SingleMethodExecute(); + } + // swap them back + m_Extreme = m_Extreme1; + m_MagnitudeSign = m_MagnitudeSign1; + m_Stage = 1; +#endif +} + +//////////////////////////////////////////////////////////// + +template +void +ParabolicOpenCloseImageFilter::ThreadedGenerateData( + const OutputImageRegionType & outputRegionForThread, + ThreadIdType threadId) +{ + // compute the number of rows first, so we can setup a progress reporter + unsigned int numberOfRows = 1; + InputSizeType size = outputRegionForThread.GetSize(); + + for (unsigned int d = 0; d < InputImageDimension; ++d) + { + if (d != m_CurrentDimension) + { + numberOfRows *= size[d]; + } + } + + float progressPerDimension = 1.0 / ImageDimension; + + ProgressReporter progress( + this, threadId, numberOfRows, 30, m_CurrentDimension * progressPerDimension, progressPerDimension); + + using InputConstIteratorType = ImageLinearConstIteratorWithIndex; + using OutputIteratorType = ImageLinearIteratorWithIndex; + + // for stages after the first + using OutputConstIteratorType = ImageLinearConstIteratorWithIndex; + + using RegionType = ImageRegion; + + typename TInputImage::ConstPointer inputImage(this->GetInput()); + typename TOutputImage::Pointer outputImage(this->GetOutput()); + + // const unsigned int imageDimension = inputImage->GetImageDimension(); + + // outputImage->SetBufferedRegion( outputImage->GetRequestedRegion() ); + // outputImage->Allocate(); + RegionType region = outputRegionForThread; + + InputConstIteratorType inputIterator(inputImage, region); + OutputIteratorType outputIterator(outputImage, region); + OutputConstIteratorType inputIteratorStage2(outputImage, region); + + if (m_Stage == 1) + { + // deal with the first dimension - this should be copied to the + // output if the scale is 0 + if (m_CurrentDimension == 0) + { + if (m_Scale[0] > 0) + { + // Perform as normal + // RealType magnitude = 1.0/(2.0 * m_Scale[0]); + unsigned long LineLength = region.GetSize()[0]; + RealType image_scale = this->GetInput()->GetSpacing()[0]; + + doOneDimension( + inputIterator, + outputIterator, + progress, + LineLength, + 0, + this->m_UseImageSpacing, + image_scale, + this->m_Scale[0], + m_ParabolicAlgorithm); + } + else + { + // copy to output + using InItType = ImageRegionConstIterator; + using OutItType = ImageRegionIterator; + + InItType InIt(inputImage, region); + OutItType OutIt(outputImage, region); + while (!InIt.IsAtEnd()) + { + OutIt.Set(static_cast(InIt.Get())); + ++InIt; + ++OutIt; + } + } + } + else + { + if (m_Scale[m_CurrentDimension] > 0) + { + // now deal with the other dimensions for first stage + unsigned long LineLength = region.GetSize()[m_CurrentDimension]; + RealType image_scale = this->GetInput()->GetSpacing()[m_CurrentDimension]; + + doOneDimension( + inputIteratorStage2, + outputIterator, + progress, + LineLength, + m_CurrentDimension, + this->m_UseImageSpacing, + image_scale, + this->m_Scale[m_CurrentDimension], + m_ParabolicAlgorithm); + } + } + } + else + { + // deal with the other dimensions for second stage + if (m_Scale[m_CurrentDimension] > 0) + { + // RealType magnitude = 1.0/(2.0 * m_Scale[dd]); + unsigned long LineLength = region.GetSize()[m_CurrentDimension]; + RealType image_scale = this->GetInput()->GetSpacing()[m_CurrentDimension]; + + doOneDimension( + inputIteratorStage2, + outputIterator, + progress, + LineLength, + m_CurrentDimension, + this->m_UseImageSpacing, + image_scale, + this->m_Scale[m_CurrentDimension], + m_ParabolicAlgorithm); + } + } +} + +template +void +ParabolicOpenCloseImageFilter::PrintSelf(std::ostream & os, Indent indent) const +{ + Superclass::PrintSelf(os, indent); + if (m_UseImageSpacing) + { + os << "Scale in world units: " << m_Scale << std::endl; + } + else + { + os << "Scale in voxels: " << m_Scale << std::endl; + } +} +} // namespace itk +#endif diff --git a/Modules/Filtering/ParabolicMorphology/include/itkParabolicOpenCloseSafeBorderImageFilter.h b/Modules/Filtering/ParabolicMorphology/include/itkParabolicOpenCloseSafeBorderImageFilter.h new file mode 100644 index 00000000000..869710e5ec3 --- /dev/null +++ b/Modules/Filtering/ParabolicMorphology/include/itkParabolicOpenCloseSafeBorderImageFilter.h @@ -0,0 +1,185 @@ +/*========================================================================= + * + * Copyright NumFOCUS + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * https://www.apache.org/licenses/LICENSE-2.0.txt + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + *=========================================================================*/ +#ifndef itkParabolicOpenCloseSafeBorderImageFilter_h +#define itkParabolicOpenCloseSafeBorderImageFilter_h + +#include "itkParabolicOpenCloseImageFilter.h" +#include "itkCropImageFilter.h" +#include "itkConstantPadImageFilter.h" +#include "itkCastImageFilter.h" +#include "itkStatisticsImageFilter.h" + +/* this class implements padding and cropping, so we don't just + * inherit from the OpenCloseImageFitler */ + +namespace itk +{ +template +class ITK_TEMPLATE_EXPORT ParabolicOpenCloseSafeBorderImageFilter : public ImageToImageFilter +{ +public: + ITK_DISALLOW_COPY_AND_MOVE(ParabolicOpenCloseSafeBorderImageFilter); + + /** Standard class type alias. */ + using Self = ParabolicOpenCloseSafeBorderImageFilter; + using Superclass = ImageToImageFilter; + using Pointer = SmartPointer; + using ConstPointer = SmartPointer; + + /** Method for creation through the object factory. */ + itkNewMacro(Self); + + /** Runtime information support. */ + itkTypeMacro(ParabolicOpenCloseSafeBorderImageFilter, ImageToImageFilter); + + /** Pixel Type of the input image */ + using InputImageType = TInputImage; + using OutputImageType = TOutputImage; + using InputPixelType = typename TInputImage::PixelType; + using RealType = typename NumericTraits::RealType; + using ScalarRealType = typename NumericTraits::ScalarRealType; + using OutputPixelType = typename TOutputImage::PixelType; + + /** Smart pointer type alias support. */ + using InputImagePointer = typename TInputImage::Pointer; + using InputImageConstPointer = typename TInputImage::ConstPointer; + + /** a type to represent the "kernel radius" */ + using RadiusType = typename itk::FixedArray; + + /** Image dimension. */ + static constexpr unsigned int ImageDimension = TInputImage::ImageDimension; + /** Define the image type for internal computations + RealType is usually 'double' in NumericTraits. + Here we prefer float in order to save memory. */ + + // set all of the scales the same + void + SetScale(ScalarRealType scale) + { + RadiusType s = this->GetScale(); + + this->m_MorphFilt->SetScale(scale); + if (s != this->GetScale()) + { + this->Modified(); + } + } + + // different scale for each direction + void + SetScale(RadiusType scale) + { + if (scale != this->GetScale()) + { + this->m_MorphFilt->SetScale(scale); + this->Modified(); + } + } + + // + const RadiusType & + GetScale() const + { + return (this->m_MorphFilt->GetScale()); + } + + void + SetUseImageSpacing(bool B) + { + if (B != this->GetUseImageSpacing()) + { + this->m_MorphFilt->SetUseImageSpacing(B); + this->Modified(); + } + } + + bool + GetUseImageSpacing() const + { + return (this->m_MorphFilt->GetUseImageSpacing()); + } + + itkBooleanMacro(UseImageSpacing); + + itkSetMacro(SafeBorder, bool); + itkGetConstReferenceMacro(SafeBorder, bool); + itkBooleanMacro(SafeBorder); + // should add the Get methods + + enum ParabolicAlgorithm + { + NOCHOICE = 0, // decices based on scale - experimental + CONTACTPOINT = 1, // sometimes faster at low scale + INTERSECTION = 2 // default + }; + /** + * Set/Get the method used. Choices are contact point or + * intersection. Intersection is the default. Contact point can be + * faster at small scales. + */ + + itkSetMacro(ParabolicAlgorithm, int); + itkGetConstReferenceMacro(ParabolicAlgorithm, int); + + /** ParabolicOpenCloseImageFilter must forward the Modified() call to its + internal filters */ + void + Modified() const override; + +protected: + void + GenerateData() override; + + void + PrintSelf(std::ostream & os, Indent indent) const override; + + using MorphFilterType = ParabolicOpenCloseImageFilter; + using PadFilterType = ConstantPadImageFilter; + using CropFilterType = CropImageFilter; + using StatsFilterType = StatisticsImageFilter; + + ParabolicOpenCloseSafeBorderImageFilter() + { + m_MorphFilt = MorphFilterType::New(); + m_PadFilt = PadFilterType::New(); + m_CropFilt = CropFilterType::New(); + m_StatsFilt = StatsFilterType::New(); + m_SafeBorder = true; + m_ParabolicAlgorithm = INTERSECTION; + } + + ~ParabolicOpenCloseSafeBorderImageFilter() override = default; + int m_ParabolicAlgorithm; + +private: + typename MorphFilterType::Pointer m_MorphFilt; + typename PadFilterType::Pointer m_PadFilt; + typename CropFilterType::Pointer m_CropFilt; + typename StatsFilterType::Pointer m_StatsFilt; + + bool m_SafeBorder; + bool m_UseContactPoint; + bool m_UseIntersection; +}; +} // end namespace itk +#ifndef ITK_MANUAL_INSTANTIATION +# include "itkParabolicOpenCloseSafeBorderImageFilter.hxx" +#endif + +#endif diff --git a/Modules/Filtering/ParabolicMorphology/include/itkParabolicOpenCloseSafeBorderImageFilter.hxx b/Modules/Filtering/ParabolicMorphology/include/itkParabolicOpenCloseSafeBorderImageFilter.hxx new file mode 100644 index 00000000000..93d4bc40c31 --- /dev/null +++ b/Modules/Filtering/ParabolicMorphology/include/itkParabolicOpenCloseSafeBorderImageFilter.hxx @@ -0,0 +1,144 @@ +/*========================================================================= + * + * Copyright NumFOCUS + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * https://www.apache.org/licenses/LICENSE-2.0.txt + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + *=========================================================================*/ +#ifndef itkParabolicOpenCloseSafeBorderImageFilter_hxx +#define itkParabolicOpenCloseSafeBorderImageFilter_hxx + +#include "itkProgressAccumulator.h" + +namespace itk +{ +template +void +ParabolicOpenCloseSafeBorderImageFilter::GenerateData() +{ + ProgressAccumulator::Pointer progress = ProgressAccumulator::New(); + + progress->SetMiniPipelineFilter(this); + + // Allocate the output + this->AllocateOutputs(); + InputImageConstPointer inputImage; + typename PadFilterType::SizeType Bounds; + typename PadFilterType::SizeType BoundsSize; + if (this->m_SafeBorder) + { + auto localInput = TInputImage::New(); + localInput->Graft(this->GetInput()); + + // need to compute some image statistics and determine the padding + // extent. This will almost certainly be an over estimate + m_StatsFilt->SetInput(localInput); + m_StatsFilt->Update(); + InputPixelType range = m_StatsFilt->GetMaximum() - m_StatsFilt->GetMinimum(); + typename MorphFilterType::RadiusType Sigma = m_MorphFilt->GetScale(); + typename TInputImage::SpacingType spcing = localInput->GetSpacing(); + for (unsigned s = 0; s < ImageDimension; s++) + { + if (m_MorphFilt->GetUseImageSpacing()) + { + RealType image_scale = spcing[s]; + Bounds[s] = (unsigned long)ceil(sqrt(2 * (Sigma[s] / (image_scale * image_scale)) * range)); + BoundsSize[s] = Bounds[s]; + } + else + { + Bounds[s] = (unsigned long)ceil(sqrt(2 * Sigma[s] * range)); + BoundsSize[s] = Bounds[s]; + } + } + m_PadFilt->SetPadLowerBound(Bounds); + m_PadFilt->SetPadUpperBound(Bounds); + + // need to select between opening and closing here + if (DoOpen) + { + // m_PadFilt->SetConstant(NumericTraits::max()); + m_PadFilt->SetConstant(m_StatsFilt->GetMaximum()); + } + else + { + // m_PadFilt->SetConstant(NumericTraits::NonpositiveMin()); + m_PadFilt->SetConstant(m_StatsFilt->GetMinimum()); + } + m_PadFilt->SetInput(localInput); + progress->RegisterInternalFilter(m_PadFilt, 0.1f); + inputImage = m_PadFilt->GetOutput(); + } + else + { + auto localInput = TInputImage::New(); + localInput->Graft(this->GetInput()); + + inputImage = localInput; + } + + m_MorphFilt->SetInput(inputImage); + m_MorphFilt->SetParabolicAlgorithm(m_ParabolicAlgorithm); + + progress->RegisterInternalFilter(m_MorphFilt, 0.8f); + + if (this->m_SafeBorder) + { + // crop + m_CropFilt->SetInput(m_MorphFilt->GetOutput()); + m_CropFilt->SetUpperBoundaryCropSize(BoundsSize); + m_CropFilt->SetLowerBoundaryCropSize(BoundsSize); + progress->RegisterInternalFilter(m_CropFilt, 0.1f); + m_CropFilt->GraftOutput(this->GetOutput()); + m_CropFilt->Update(); + this->GraftOutput(m_CropFilt->GetOutput()); + } + else + { + m_MorphFilt->GraftOutput(this->GetOutput()); + m_MorphFilt->Update(); + this->GraftOutput(m_MorphFilt->GetOutput()); + // std::cout << "Finished grafting" << std::endl; + } +} + +template +void +ParabolicOpenCloseSafeBorderImageFilter::Modified() const +{ + Superclass::Modified(); + m_MorphFilt->Modified(); + m_PadFilt->Modified(); + m_CropFilt->Modified(); + m_StatsFilt->Modified(); +} + +/////////////////////////////////// +template +void +ParabolicOpenCloseSafeBorderImageFilter::PrintSelf(std::ostream & os, + Indent indent) const +{ + os << indent << "SafeBorder: " << m_SafeBorder << std::endl; + if (this->GetUseImageSpacing()) + { + os << "Scale in world units: " << this->GetScale() << std::endl; + } + else + { + os << "Scale in voxels: " << this->GetScale() << std::endl; + } +} +} // namespace itk + +#endif diff --git a/Modules/Filtering/ParabolicMorphology/include/itkParabolicOpenImageFilter.h b/Modules/Filtering/ParabolicMorphology/include/itkParabolicOpenImageFilter.h new file mode 100644 index 00000000000..f3fbb5d2d17 --- /dev/null +++ b/Modules/Filtering/ParabolicMorphology/include/itkParabolicOpenImageFilter.h @@ -0,0 +1,91 @@ +/*========================================================================= + * + * Copyright NumFOCUS + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * https://www.apache.org/licenses/LICENSE-2.0.txt + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + *=========================================================================*/ +#ifndef itkParabolicOpenImageFilter_h +#define itkParabolicOpenImageFilter_h + +#include "itkParabolicOpenCloseSafeBorderImageFilter.h" +#include "itkNumericTraits.h" + +namespace itk +{ +/** + * \class ParabolicOpenImageFilter + * \brief Class for morphological opening + * operations with parabolic structuring elements. + * + * This filter provides options for padded borders + * + * This filter is threaded. + * + * Core methods described in the InsightJournal article: + * "Morphology with parabolic structuring elements" + * + * https://hdl.handle.net/1926/1370 + * + * \sa itkParabolicOpenCloseImageFilter + * + * \ingroup ParabolicMorphology + * + * \author Richard Beare, Department of Medicine, Monash University, + * Australia. + **/ + +template +class ITK_TEMPLATE_EXPORT ParabolicOpenImageFilter + : public ParabolicOpenCloseSafeBorderImageFilter +{ +public: + ITK_DISALLOW_COPY_AND_MOVE(ParabolicOpenImageFilter); + + /** Standard class type alias. */ + using Self = ParabolicOpenImageFilter; + using Superclass = ParabolicOpenCloseSafeBorderImageFilter; + using Pointer = SmartPointer; + using ConstPointer = SmartPointer; + + /** Method for creation through the object factory. */ + itkNewMacro(Self); + + /** Pixel Type of the input image */ + using InputImageType = TInputImage; + using OutputImageType = TOutputImage; + using PixelType = typename TInputImage::PixelType; + using RealType = typename NumericTraits::RealType; + using ScalarRealType = typename NumericTraits::ScalarRealType; + using OutputPixelType = typename TOutputImage::PixelType; + + /** Smart pointer type alias support. */ + using InputImagePointer = typename TInputImage::Pointer; + using InputImageConstPointer = typename TInputImage::ConstPointer; + + /** a type to represent the "kernel radius" */ + using RadiusType = typename itk::FixedArray; + + /** Image dimension. */ + static constexpr unsigned int ImageDimension = TInputImage::ImageDimension; + /** Define the image type for internal computations + RealType is usually 'double' in NumericTraits. + Here we prefer float in order to save memory. */ +protected: + ParabolicOpenImageFilter() = default; + ~ParabolicOpenImageFilter() override = default; + // void PrintSelf(std::ostream& os, Indent indent) const; +}; +} // end namespace itk + +#endif diff --git a/Modules/Filtering/ParabolicMorphology/include/itkSharpenOpImageFilter.h b/Modules/Filtering/ParabolicMorphology/include/itkSharpenOpImageFilter.h new file mode 100644 index 00000000000..9412662bd0d --- /dev/null +++ b/Modules/Filtering/ParabolicMorphology/include/itkSharpenOpImageFilter.h @@ -0,0 +1,126 @@ +/*========================================================================= + * + * Copyright NumFOCUS + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * https://www.apache.org/licenses/LICENSE-2.0.txt + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + *=========================================================================*/ +#ifndef itkSharpenOpImageFilter_h +#define itkSharpenOpImageFilter_h + +#include "itkTernaryFunctorImageFilter.h" + +namespace itk +{ +/** \class SharpenOpImageFilter + * \brief Implements the sharpening operation. The inputs are the + * dilated, eroded and original images. + * + * This class is parametrized over the types of the three + * input images and the type of the output image. + * Numeric conversions (castings) are done by the C++ defaults. + * + * In reality the input and output types of this filter are expected + * to be the same. + * + * Core methods described in the InsightJournal article: + * "Morphology with parabolic structuring elements" + * + * https://hdl.handle.net/1926/1370 + * \ingroup ParabolicMorphology + * + * \author Richard Beare, Department of Medicine, Monash University, + * Australia. + * + */ +namespace Function +{ +template +class SharpM +{ +public: + SharpM() = default; + ~SharpM() = default; + bool + operator!=(const SharpM &) const + { + return false; + } + + bool + operator==(const SharpM & other) const + { + return !(*this != other); + } + + inline TOutput + operator()(const TInput1 & A, const TInput2 & B, const TInput3 & C) + { + // the sharpening operator. A is the dilation, B the original, C + // the erosion + TInput2 diff1 = A - B; + TInput2 diff2 = B - C; + + if (diff1 < diff2) + { + return (TOutput)A; + } + if (diff2 < diff1) + { + return (TOutput)C; + } + return ((TOutput)B); + } +}; +} // namespace Function + +template +class ITK_TEMPLATE_EXPORT SharpenOpImageFilter + : public TernaryFunctorImageFilter> +{ +public: + ITK_DISALLOW_COPY_AND_MOVE(SharpenOpImageFilter); + + /** Standard class type alias. */ + using Self = SharpenOpImageFilter; + using Superclass = TernaryFunctorImageFilter>; + using Pointer = SmartPointer; + using ConstPointer = SmartPointer; + + /** Method for creation through the object factory. */ + itkNewMacro(Self); + + /** Runtime information support. */ + itkTypeMacro(SharpenOpImageFilter, TernaryFunctorImageFilter); + +protected: + SharpenOpImageFilter() = default; + ~SharpenOpImageFilter() override = default; +}; +} // end namespace itk + +#endif diff --git a/Modules/Filtering/ParabolicMorphology/itk-module.cmake b/Modules/Filtering/ParabolicMorphology/itk-module.cmake new file mode 100644 index 00000000000..969625a2ae9 --- /dev/null +++ b/Modules/Filtering/ParabolicMorphology/itk-module.cmake @@ -0,0 +1,15 @@ +# Maintainer: Richard Beare +itk_module( + ParabolicMorphology + DEPENDS + ITKIOImageBase + ITKThresholding + TEST_DEPENDS + ITKImageGrid + ITKTestKernel + ITKMathematicalMorphology + ITKSmoothing + EXCLUDE_FROM_DEFAULT + DESCRIPTION + "Classes performing morphology using parabolic functions. Fast distance transforms and binary erosions/dilations/openings/closings by spheres, sharpenings and grayscale operations. https://doi.org/10.54294/aq68pt" +) diff --git a/Modules/Filtering/ParabolicMorphology/test/CMakeLists.txt b/Modules/Filtering/ParabolicMorphology/test/CMakeLists.txt new file mode 100644 index 00000000000..ba83d9ce1ea --- /dev/null +++ b/Modules/Filtering/ParabolicMorphology/test/CMakeLists.txt @@ -0,0 +1,292 @@ +itk_module_test() +set( + ParabolicMorphologyTests + itkParaErodeTest.cxx + itkParaDilateTest.cxx + itkParaOpenTest.cxx + itkParaSpacingTest.cxx + itkParaSharpenTest.cxx + itkParaDTTest.cxx + itkBinaryDilateParabolicTest.cxx + itkBinaryErodeParabolicTest.cxx + itkBinaryOpenParabolicTest.cxx + itkBinaryCloseParabolicTest.cxx +) + +set(INPUT_IMAGE ${CMAKE_CURRENT_SOURCE_DIR}/images/cthead1.png) + +set(ITK_TEST_DRIVER itkTestDriver) + +createtestdriver(ParabolicMorphology "${ParabolicMorphology-Test_LIBRARIES}" "${ParabolicMorphologyTests}") + +## both intersection and contact point erosion +## default scale +itk_add_test( + NAME itkParaErodeTest2D_1 + COMMAND + ParabolicMorphologyTestDriver + --compare + outEInta.png + outECPa.png + --compare + outEInta.png + ${CMAKE_CURRENT_SOURCE_DIR}/baseline/outECPa.png + itkParaErodeTest + ${INPUT_IMAGE} + outEInta.png + outECPa.png +) + +## small scale +itk_add_test( + NAME itkParaErodeTest2D_2 + COMMAND + ParabolicMorphologyTestDriver + --compare + outEIntb.png + outECPb.png + --compare + outEIntb.png + ${CMAKE_CURRENT_SOURCE_DIR}/baseline/outECPb.png + itkParaErodeTest + ${INPUT_IMAGE} + outEIntb.png + outECPb.png + 0.2 +) + +## large scale +itk_add_test( + NAME itkParaErodeTest2D_3 + COMMAND + ParabolicMorphologyTestDriver + --compare + outEIntc.png + outECPc.png + --compare + outEIntc.png + ${CMAKE_CURRENT_SOURCE_DIR}/baseline/outECPc.png + itkParaErodeTest + ${INPUT_IMAGE} + outEIntc.png + outECPc.png + 5 +) + +## default scale +itk_add_test( + NAME itkParaDilateTest2D_1 + COMMAND + ParabolicMorphologyTestDriver + --compare + outDInta.png + outDCPa.png + --compare + outDInta.png + ${CMAKE_CURRENT_SOURCE_DIR}/baseline/outDCPa.png + itkParaDilateTest + ${INPUT_IMAGE} + outDInta.png + outDCPa.png +) + +## small scale +itk_add_test( + NAME itkParaDilateTest2D_2 + COMMAND + ParabolicMorphologyTestDriver + --compare + outDIntb.png + outDCPb.png + --compare + outDIntb.png + ${CMAKE_CURRENT_SOURCE_DIR}/baseline/outDCPb.png + itkParaDilateTest + ${INPUT_IMAGE} + outDIntb.png + outDCPb.png + 0.2 +) + +## large scale +itk_add_test( + NAME itkParaDilateTest2D_3 + COMMAND + ParabolicMorphologyTestDriver + --compare + outDIntc.png + outDCPc.png + --compare + outDIntc.png + ${CMAKE_CURRENT_SOURCE_DIR}/baseline/outDCPc.png + itkParaDilateTest + ${INPUT_IMAGE} + outDIntc.png + outDCPc.png + 5 +) + +## Opening +itk_add_test( + NAME itkParaOpenTest2D_1 + COMMAND + ParabolicMorphologyTestDriver + --compare + openCP.png + openInt.png + --compare + openInt.png + ${CMAKE_CURRENT_SOURCE_DIR}/baseline/openInt.png + itkParaOpenTest + ${INPUT_IMAGE} + openInt.png + openCP.png +) + +itk_add_test( + NAME itkParaSpacingTest2D_1 + COMMAND + ParabolicMorphologyTestDriver + --compare + osp1.png + osp2.png + itkParaSpacingTest + ${INPUT_IMAGE} + osp1.png + osp2.png +) + +## sharpening +itk_add_test( + NAME itkParaSharpenTest2D_1 + COMMAND + ParabolicMorphologyTestDriver + --compare + sharp1.mha + ${CMAKE_CURRENT_SOURCE_DIR}/baseline/sharp1.mha + itkParaSharpenTest + 1 + sharp1.mha +) + +itk_add_test( + NAME itkParaSharpenTest2D_3 + COMMAND + ParabolicMorphologyTestDriver + --compare + sharp3.mha + ${CMAKE_CURRENT_SOURCE_DIR}/baseline/sharp3.mha + itkParaSharpenTest + 3 + sharp3.mha +) + +itk_add_test( + NAME itkParaSharpenTest2D_10 + COMMAND + ParabolicMorphologyTestDriver + --compare + sharp10.mha + ${CMAKE_CURRENT_SOURCE_DIR}/baseline/sharp10.mha + itkParaSharpenTest + 10 + sharp10.mha +) + +itk_add_test( + NAME itkParaSharpenTest2D_100 + COMMAND + ParabolicMorphologyTestDriver + --compare + sharp100.mha + ${CMAKE_CURRENT_SOURCE_DIR}/baseline/sharp100.mha + itkParaSharpenTest + 100 + sharp100.mha +) + +## distance transform +itk_add_test( + NAME itkParaDTTest2D_0 + COMMAND + ParabolicMorphologyTestDriver + --compare + dist1.mha + ${CMAKE_CURRENT_SOURCE_DIR}/baseline/dist1.mha + itkParaDTTest + ${INPUT_IMAGE} + 100 + 0 + dist1.mha +) + +itk_add_test( + NAME itkParaDTTest2D_255 + COMMAND + ParabolicMorphologyTestDriver + --compare + dist255.mha + ${CMAKE_CURRENT_SOURCE_DIR}/baseline/dist255.mha + itkParaDTTest + ${INPUT_IMAGE} + 100 + 255 + dist255.mha +) + +## Binary morphology +itk_add_test( + NAME itkBinaryDilateParabolic2D_10 + COMMAND + ParabolicMorphologyTestDriver + --compare + dilatebinary10.mha + ${CMAKE_CURRENT_SOURCE_DIR}/baseline/dilatebinary10.mha + itkBinaryDilateParabolicTest + ${INPUT_IMAGE} + 150 + 10 + dilatebinary10.mha +) + +itk_add_test( + NAME itkBinaryErodeParabolic2D_10 + COMMAND + ParabolicMorphologyTestDriver + --compare + erodebinary10.mha + ${CMAKE_CURRENT_SOURCE_DIR}/baseline/erodebinary10.mha + itkBinaryErodeParabolicTest + ${INPUT_IMAGE} + 120 + 10 + erodebinary10.mha +) + +itk_add_test( + NAME itkBinaryOpenParabolic2D_10 + COMMAND + ParabolicMorphologyTestDriver + --compare + openbinary10.mha + ${CMAKE_CURRENT_SOURCE_DIR}/baseline/openbinary10.mha + itkBinaryOpenParabolicTest + ${INPUT_IMAGE} + 130 + 10 + openbinary10.mha +) + +itk_add_test( + NAME itkBinaryCloseParabolic2D_10 + COMMAND + ParabolicMorphologyTestDriver + --compare + closebinary10.mha + ${CMAKE_CURRENT_SOURCE_DIR}/baseline/closebinary10.mha + itkBinaryCloseParabolicTest + ${INPUT_IMAGE} + 150 + 10 + closebinary10.mha +) diff --git a/Modules/Filtering/ParabolicMorphology/test/baseline/closebinary10.mha b/Modules/Filtering/ParabolicMorphology/test/baseline/closebinary10.mha new file mode 100644 index 00000000000..9bbb6100512 Binary files /dev/null and b/Modules/Filtering/ParabolicMorphology/test/baseline/closebinary10.mha differ diff --git a/Modules/Filtering/ParabolicMorphology/test/baseline/dilatebinary10.mha b/Modules/Filtering/ParabolicMorphology/test/baseline/dilatebinary10.mha new file mode 100644 index 00000000000..fceb59521dd Binary files /dev/null and b/Modules/Filtering/ParabolicMorphology/test/baseline/dilatebinary10.mha differ diff --git a/Modules/Filtering/ParabolicMorphology/test/baseline/dist1.mha b/Modules/Filtering/ParabolicMorphology/test/baseline/dist1.mha new file mode 100644 index 00000000000..845174e7e7a Binary files /dev/null and b/Modules/Filtering/ParabolicMorphology/test/baseline/dist1.mha differ diff --git a/Modules/Filtering/ParabolicMorphology/test/baseline/dist255.mha b/Modules/Filtering/ParabolicMorphology/test/baseline/dist255.mha new file mode 100644 index 00000000000..39f829ce09f Binary files /dev/null and b/Modules/Filtering/ParabolicMorphology/test/baseline/dist255.mha differ diff --git a/Modules/Filtering/ParabolicMorphology/test/baseline/erodebinary10.mha b/Modules/Filtering/ParabolicMorphology/test/baseline/erodebinary10.mha new file mode 100644 index 00000000000..6781879d28b Binary files /dev/null and b/Modules/Filtering/ParabolicMorphology/test/baseline/erodebinary10.mha differ diff --git a/Modules/Filtering/ParabolicMorphology/test/baseline/openInt.png b/Modules/Filtering/ParabolicMorphology/test/baseline/openInt.png new file mode 100644 index 00000000000..3c4d3e8765d Binary files /dev/null and b/Modules/Filtering/ParabolicMorphology/test/baseline/openInt.png differ diff --git a/Modules/Filtering/ParabolicMorphology/test/baseline/openbinary10.mha b/Modules/Filtering/ParabolicMorphology/test/baseline/openbinary10.mha new file mode 100644 index 00000000000..48e3301cbb0 Binary files /dev/null and b/Modules/Filtering/ParabolicMorphology/test/baseline/openbinary10.mha differ diff --git a/Modules/Filtering/ParabolicMorphology/test/baseline/outDCPa.png b/Modules/Filtering/ParabolicMorphology/test/baseline/outDCPa.png new file mode 100644 index 00000000000..e8fbad578a0 Binary files /dev/null and b/Modules/Filtering/ParabolicMorphology/test/baseline/outDCPa.png differ diff --git a/Modules/Filtering/ParabolicMorphology/test/baseline/outDCPb.png b/Modules/Filtering/ParabolicMorphology/test/baseline/outDCPb.png new file mode 100644 index 00000000000..3be55bf5ed3 Binary files /dev/null and b/Modules/Filtering/ParabolicMorphology/test/baseline/outDCPb.png differ diff --git a/Modules/Filtering/ParabolicMorphology/test/baseline/outDCPc.png b/Modules/Filtering/ParabolicMorphology/test/baseline/outDCPc.png new file mode 100644 index 00000000000..04e76fcf7e9 Binary files /dev/null and b/Modules/Filtering/ParabolicMorphology/test/baseline/outDCPc.png differ diff --git a/Modules/Filtering/ParabolicMorphology/test/baseline/outDInta.png b/Modules/Filtering/ParabolicMorphology/test/baseline/outDInta.png new file mode 100644 index 00000000000..e8fbad578a0 Binary files /dev/null and b/Modules/Filtering/ParabolicMorphology/test/baseline/outDInta.png differ diff --git a/Modules/Filtering/ParabolicMorphology/test/baseline/outDIntb.png b/Modules/Filtering/ParabolicMorphology/test/baseline/outDIntb.png new file mode 100644 index 00000000000..3be55bf5ed3 Binary files /dev/null and b/Modules/Filtering/ParabolicMorphology/test/baseline/outDIntb.png differ diff --git a/Modules/Filtering/ParabolicMorphology/test/baseline/outDIntc.png b/Modules/Filtering/ParabolicMorphology/test/baseline/outDIntc.png new file mode 100644 index 00000000000..04e76fcf7e9 Binary files /dev/null and b/Modules/Filtering/ParabolicMorphology/test/baseline/outDIntc.png differ diff --git a/Modules/Filtering/ParabolicMorphology/test/baseline/outECPa.png b/Modules/Filtering/ParabolicMorphology/test/baseline/outECPa.png new file mode 100644 index 00000000000..3f2b8e5395e Binary files /dev/null and b/Modules/Filtering/ParabolicMorphology/test/baseline/outECPa.png differ diff --git a/Modules/Filtering/ParabolicMorphology/test/baseline/outECPb.png b/Modules/Filtering/ParabolicMorphology/test/baseline/outECPb.png new file mode 100644 index 00000000000..3713937daf7 Binary files /dev/null and b/Modules/Filtering/ParabolicMorphology/test/baseline/outECPb.png differ diff --git a/Modules/Filtering/ParabolicMorphology/test/baseline/outECPc.png b/Modules/Filtering/ParabolicMorphology/test/baseline/outECPc.png new file mode 100644 index 00000000000..172baa3734b Binary files /dev/null and b/Modules/Filtering/ParabolicMorphology/test/baseline/outECPc.png differ diff --git a/Modules/Filtering/ParabolicMorphology/test/baseline/outEInta.png b/Modules/Filtering/ParabolicMorphology/test/baseline/outEInta.png new file mode 100644 index 00000000000..3f2b8e5395e Binary files /dev/null and b/Modules/Filtering/ParabolicMorphology/test/baseline/outEInta.png differ diff --git a/Modules/Filtering/ParabolicMorphology/test/baseline/outEIntb.png b/Modules/Filtering/ParabolicMorphology/test/baseline/outEIntb.png new file mode 100644 index 00000000000..3713937daf7 Binary files /dev/null and b/Modules/Filtering/ParabolicMorphology/test/baseline/outEIntb.png differ diff --git a/Modules/Filtering/ParabolicMorphology/test/baseline/outEIntc.png b/Modules/Filtering/ParabolicMorphology/test/baseline/outEIntc.png new file mode 100644 index 00000000000..172baa3734b Binary files /dev/null and b/Modules/Filtering/ParabolicMorphology/test/baseline/outEIntc.png differ diff --git a/Modules/Filtering/ParabolicMorphology/test/baseline/sharp1.mha b/Modules/Filtering/ParabolicMorphology/test/baseline/sharp1.mha new file mode 100644 index 00000000000..61e40dca387 Binary files /dev/null and b/Modules/Filtering/ParabolicMorphology/test/baseline/sharp1.mha differ diff --git a/Modules/Filtering/ParabolicMorphology/test/baseline/sharp10.mha b/Modules/Filtering/ParabolicMorphology/test/baseline/sharp10.mha new file mode 100644 index 00000000000..1b93d8d631f Binary files /dev/null and b/Modules/Filtering/ParabolicMorphology/test/baseline/sharp10.mha differ diff --git a/Modules/Filtering/ParabolicMorphology/test/baseline/sharp100.mha b/Modules/Filtering/ParabolicMorphology/test/baseline/sharp100.mha new file mode 100644 index 00000000000..1b93d8d631f Binary files /dev/null and b/Modules/Filtering/ParabolicMorphology/test/baseline/sharp100.mha differ diff --git a/Modules/Filtering/ParabolicMorphology/test/baseline/sharp3.mha b/Modules/Filtering/ParabolicMorphology/test/baseline/sharp3.mha new file mode 100644 index 00000000000..bf312a664e8 Binary files /dev/null and b/Modules/Filtering/ParabolicMorphology/test/baseline/sharp3.mha differ diff --git a/Modules/Filtering/ParabolicMorphology/test/images/blurredprof.txt b/Modules/Filtering/ParabolicMorphology/test/images/blurredprof.txt new file mode 100644 index 00000000000..2872700d239 --- /dev/null +++ b/Modules/Filtering/ParabolicMorphology/test/images/blurredprof.txt @@ -0,0 +1,100 @@ +0 0 +1 0 +2 0 +3 0 +4 0 +5 0 +6 0 +7 0 +8 0 +9 0 +10 0 +11 0 +12 0 +13 0 +14 0 +15 0 +16 0 +17 0 +18 0 +19 0 +20 0 +21 0 +22 0 +23 0 +24 0 +25 0 +26 0 +27 0 +28 0 +29 0 +30 0 +31 0 +32 0 +33 0 +34 0 +35 0 +36 2 +37 5 +38 12 +39 23 +40 40 +41 63 +42 92 +43 124 +44 156 +45 184 +46 208 +47 225 +48 237 +49 243 +50 245 +51 243 +52 237 +53 225 +54 208 +55 184 +56 156 +57 124 +58 92 +59 63 +60 40 +61 23 +62 12 +63 5 +64 2 +65 0 +66 0 +67 0 +68 0 +69 0 +70 0 +71 0 +72 0 +73 0 +74 0 +75 0 +76 0 +77 0 +78 0 +79 0 +80 0 +81 0 +82 0 +83 0 +84 0 +85 0 +86 0 +87 0 +88 0 +89 0 +90 0 +91 0 +92 0 +93 0 +94 0 +95 0 +96 0 +97 0 +98 0 +99 0 diff --git a/Modules/Filtering/ParabolicMorphology/test/images/bunnyPadded.nrrd b/Modules/Filtering/ParabolicMorphology/test/images/bunnyPadded.nrrd new file mode 100644 index 00000000000..03e2ad5638f Binary files /dev/null and b/Modules/Filtering/ParabolicMorphology/test/images/bunnyPadded.nrrd differ diff --git a/Modules/Filtering/ParabolicMorphology/test/images/cthead1.png b/Modules/Filtering/ParabolicMorphology/test/images/cthead1.png new file mode 100644 index 00000000000..8c40498a770 Binary files /dev/null and b/Modules/Filtering/ParabolicMorphology/test/images/cthead1.png differ diff --git a/Modules/Filtering/ParabolicMorphology/test/images/dilate_parabolic.fig b/Modules/Filtering/ParabolicMorphology/test/images/dilate_parabolic.fig new file mode 100644 index 00000000000..423e484c709 --- /dev/null +++ b/Modules/Filtering/ParabolicMorphology/test/images/dilate_parabolic.fig @@ -0,0 +1,74 @@ +#FIG 3.2 +Landscape +Center +Inches +Letter +100.00 +Single +-2 +1200 2 +6 2850 300 8850 3300 +3 2 0 1 0 7 50 -1 -1 0.000 0 0 0 6 + 5850 3300 6450 3150 7050 2850 7650 2250 8250 1350 8850 300 + 0.000 -1.000 -1.000 -1.000 -1.000 0.000 +3 2 0 1 0 7 50 -1 -1 0.000 0 0 0 6 + 5850 3300 5250 3150 4650 2850 4050 2250 3450 1350 2850 300 + 0.000 -1.000 -1.000 -1.000 -1.000 0.000 +-6 +2 2 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 5 + 2100 1800 2400 1800 2400 7200 2100 7200 2100 1800 +2 2 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 5 + 2400 2400 2700 2400 2700 7200 2400 7200 2400 2400 +2 2 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 5 + 2700 2700 3000 2700 3000 7200 2700 7200 2700 2700 +2 2 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 5 + 3000 2400 3300 2400 3300 7200 3000 7200 3000 2400 +2 2 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 5 + 3600 7200 3300 7200 3300 2400 3600 2400 3600 7200 +2 2 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 5 + 3900 7200 3600 7200 3600 2100 3900 2100 3900 7200 +2 2 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 5 + 3900 3000 4200 3000 4200 7200 3900 7200 3900 3000 +2 2 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 5 + 4200 5325 4500 5325 4500 7200 4200 7200 4200 5325 +2 2 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 5 + 4500 5100 4800 5100 4800 7200 4500 7200 4500 5100 +2 2 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 5 + 4800 5025 5100 5025 5100 7200 4800 7200 4800 5025 +2 2 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 5 + 5100 4800 5400 4800 5400 7200 5100 7200 5100 4800 +2 2 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 5 + 5400 5100 5700 5100 5700 7200 5400 7200 5400 5100 +2 2 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 5 + 5700 4500 6000 4500 6000 7200 5700 7200 5700 4500 +2 2 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 5 + 6300 7200 6000 7200 6000 5475 6300 5475 6300 7200 +2 2 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 5 + 6300 5100 6600 5100 6600 7200 6300 7200 6300 5100 +2 2 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 5 + 6600 3675 6900 3675 6900 7200 6600 7200 6600 3675 +2 2 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 5 + 6900 3900 7200 3900 7200 7200 6900 7200 6900 3900 +2 2 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 5 + 7200 3675 7500 3675 7500 7200 7200 7200 7200 3675 +2 2 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 5 + 7500 3000 7800 3000 7800 7200 7500 7200 7500 3000 +2 2 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 5 + 7800 2925 8100 2925 8100 7200 7800 7200 7800 2925 +2 2 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 5 + 8100 2775 8400 2775 8400 7200 8100 7200 8100 2775 +2 2 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 5 + 8400 2700 8700 2700 8700 7200 8400 7200 8400 2700 +2 2 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 5 + 8700 2700 9000 2700 9000 7200 8700 7200 8700 2700 +2 2 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 5 + 9000 2700 9300 2700 9300 7200 9000 7200 9000 2700 +2 1 0 3 0 7 50 -1 -1 0.000 0 0 -1 1 0 2 + 2 1 3.00 90.00 150.00 + 5850 1350 5850 3300 +4 0 0 50 -1 0 16 0.0000 4 225 2805 8775 900 parabolic tructuring element\001 +4 0 0 50 -1 0 12 0.0000 4 180 1500 5475 825 Dilated signal value\001 +4 0 0 50 -1 0 12 0.0000 4 135 135 5775 7500 A\001 +4 0 0 50 -1 0 12 0.0000 4 135 120 3975 1950 B\001 +4 0 0 50 -1 0 12 0.0000 4 135 120 5850 3525 C\001 +4 0 0 50 -1 0 16 0.0000 4 225 630 2100 7500 Signal\001 diff --git a/Modules/Filtering/ParabolicMorphology/test/images/inputprof.txt b/Modules/Filtering/ParabolicMorphology/test/images/inputprof.txt new file mode 100644 index 00000000000..9c046c18d08 --- /dev/null +++ b/Modules/Filtering/ParabolicMorphology/test/images/inputprof.txt @@ -0,0 +1,100 @@ +0 0 +1 0 +2 0 +3 0 +4 0 +5 0 +6 0 +7 0 +8 0 +9 0 +10 0 +11 0 +12 0 +13 0 +14 0 +15 0 +16 0 +17 0 +18 0 +19 0 +20 0 +21 0 +22 0 +23 0 +24 0 +25 0 +26 0 +27 0 +28 0 +29 0 +30 0 +31 0 +32 0 +33 0 +34 0 +35 0 +36 0 +37 0 +38 0 +39 0 +40 0 +41 0 +42 0 +43 255 +44 255 +45 255 +46 255 +47 255 +48 255 +49 255 +50 255 +51 255 +52 255 +53 255 +54 255 +55 255 +56 255 +57 255 +58 0 +59 0 +60 0 +61 0 +62 0 +63 0 +64 0 +65 0 +66 0 +67 0 +68 0 +69 0 +70 0 +71 0 +72 0 +73 0 +74 0 +75 0 +76 0 +77 0 +78 0 +79 0 +80 0 +81 0 +82 0 +83 0 +84 0 +85 0 +86 0 +87 0 +88 0 +89 0 +90 0 +91 0 +92 0 +93 0 +94 0 +95 0 +96 0 +97 0 +98 0 +99 0 diff --git a/Modules/Filtering/ParabolicMorphology/test/images/open.png b/Modules/Filtering/ParabolicMorphology/test/images/open.png new file mode 100644 index 00000000000..1b6e75a60c2 Binary files /dev/null and b/Modules/Filtering/ParabolicMorphology/test/images/open.png differ diff --git a/Modules/Filtering/ParabolicMorphology/test/images/sharp.gnuplot b/Modules/Filtering/ParabolicMorphology/test/images/sharp.gnuplot new file mode 100644 index 00000000000..0e7ac70040a --- /dev/null +++ b/Modules/Filtering/ParabolicMorphology/test/images/sharp.gnuplot @@ -0,0 +1,10 @@ +set term postscript eps enhanced color +set xlabel "image index" +set ylabel "gray level" +plot[0:100][0:255] "../images/inputprof.txt" using 1:2 with linespoints title "Raw input",\ +"../images/blurredprof.txt" using 1:2 with lines title "Blurred input",\ +"../images/sharp1.txt" using 1:2 with lines title "Sharpening 1 iteration",\ +"../images/sharp2.txt" using 1:2 with lines title "Sharpening 2 iterations",\ +"../images/sharp3.txt" using 1:2 with lines title "Sharpening 3 iterations",\ +"../images/sharp10.txt" using 1:2 with lines title "Sharpening 10 iterations",\ +"../images/sharp100.txt" using 1:2 with lines title "Sharpening 100 iterations" diff --git a/Modules/Filtering/ParabolicMorphology/test/images/sharp1.txt b/Modules/Filtering/ParabolicMorphology/test/images/sharp1.txt new file mode 100644 index 00000000000..453ae9bd02d --- /dev/null +++ b/Modules/Filtering/ParabolicMorphology/test/images/sharp1.txt @@ -0,0 +1,100 @@ +0 0 +1 0 +2 0 +3 0 +4 0 +5 0 +6 0 +7 0 +8 0 +9 0 +10 0 +11 0 +12 0 +13 0 +14 0 +15 0 +16 0 +17 0 +18 0 +19 0 +20 0 +21 0 +22 0 +23 0 +24 0 +25 0 +26 0 +27 0 +28 0 +29 0 +30 0 +31 0 +32 0 +33 0 +34 0 +35 0 +36 0.5 +37 2 +38 4 +39 6.5 +40 9.5 +41 13 +42 17.5 +43 124 +44 230.5 +45 235 +46 238.5 +47 241 +48 243 +49 244.5 +50 245 +51 244.5 +52 243 +53 241 +54 238.5 +55 235 +56 230.5 +57 124 +58 17.5 +59 13 +60 9.5 +61 6.5 +62 4 +63 2 +64 0.5 +65 0 +66 0 +67 0 +68 0 +69 0 +70 0 +71 0 +72 0 +73 0 +74 0 +75 0 +76 0 +77 0 +78 0 +79 0 +80 0 +81 0 +82 0 +83 0 +84 0 +85 0 +86 0 +87 0 +88 0 +89 0 +90 0 +91 0 +92 0 +93 0 +94 0 +95 0 +96 0 +97 0 +98 0 +99 0 diff --git a/Modules/Filtering/ParabolicMorphology/test/images/sharp10.txt b/Modules/Filtering/ParabolicMorphology/test/images/sharp10.txt new file mode 100644 index 00000000000..f4fc288a648 --- /dev/null +++ b/Modules/Filtering/ParabolicMorphology/test/images/sharp10.txt @@ -0,0 +1,100 @@ +0 0 +1 0 +2 0 +3 0 +4 0 +5 0 +6 0 +7 0 +8 0 +9 0 +10 0 +11 0 +12 0 +13 0 +14 0 +15 0 +16 0 +17 0 +18 0 +19 0 +20 0 +21 0 +22 0 +23 0 +24 0 +25 0 +26 0 +27 0 +28 0 +29 0 +30 0 +31 0 +32 0 +33 0 +34 0 +35 0 +36 0.5 +37 1 +38 1.5 +39 2 +40 2.5 +41 3 +42 3.5 +43 241.5 +44 242 +45 242.5 +46 243 +47 243.5 +48 244 +49 244.5 +50 245 +51 244.5 +52 244 +53 243.5 +54 243 +55 242.5 +56 242 +57 241.5 +58 3.5 +59 3 +60 2.5 +61 2 +62 1.5 +63 1 +64 0.5 +65 0 +66 0 +67 0 +68 0 +69 0 +70 0 +71 0 +72 0 +73 0 +74 0 +75 0 +76 0 +77 0 +78 0 +79 0 +80 0 +81 0 +82 0 +83 0 +84 0 +85 0 +86 0 +87 0 +88 0 +89 0 +90 0 +91 0 +92 0 +93 0 +94 0 +95 0 +96 0 +97 0 +98 0 +99 0 diff --git a/Modules/Filtering/ParabolicMorphology/test/images/sharp100.txt b/Modules/Filtering/ParabolicMorphology/test/images/sharp100.txt new file mode 100644 index 00000000000..f4fc288a648 --- /dev/null +++ b/Modules/Filtering/ParabolicMorphology/test/images/sharp100.txt @@ -0,0 +1,100 @@ +0 0 +1 0 +2 0 +3 0 +4 0 +5 0 +6 0 +7 0 +8 0 +9 0 +10 0 +11 0 +12 0 +13 0 +14 0 +15 0 +16 0 +17 0 +18 0 +19 0 +20 0 +21 0 +22 0 +23 0 +24 0 +25 0 +26 0 +27 0 +28 0 +29 0 +30 0 +31 0 +32 0 +33 0 +34 0 +35 0 +36 0.5 +37 1 +38 1.5 +39 2 +40 2.5 +41 3 +42 3.5 +43 241.5 +44 242 +45 242.5 +46 243 +47 243.5 +48 244 +49 244.5 +50 245 +51 244.5 +52 244 +53 243.5 +54 243 +55 242.5 +56 242 +57 241.5 +58 3.5 +59 3 +60 2.5 +61 2 +62 1.5 +63 1 +64 0.5 +65 0 +66 0 +67 0 +68 0 +69 0 +70 0 +71 0 +72 0 +73 0 +74 0 +75 0 +76 0 +77 0 +78 0 +79 0 +80 0 +81 0 +82 0 +83 0 +84 0 +85 0 +86 0 +87 0 +88 0 +89 0 +90 0 +91 0 +92 0 +93 0 +94 0 +95 0 +96 0 +97 0 +98 0 +99 0 diff --git a/Modules/Filtering/ParabolicMorphology/test/images/sharp2.txt b/Modules/Filtering/ParabolicMorphology/test/images/sharp2.txt new file mode 100644 index 00000000000..feddb667388 --- /dev/null +++ b/Modules/Filtering/ParabolicMorphology/test/images/sharp2.txt @@ -0,0 +1,100 @@ +0 0 +1 0 +2 0 +3 0 +4 0 +5 0 +6 0 +7 0 +8 0 +9 0 +10 0 +11 0 +12 0 +13 0 +14 0 +15 0 +16 0 +17 0 +18 0 +19 0 +20 0 +21 0 +22 0 +23 0 +24 0 +25 0 +26 0 +27 0 +28 0 +29 0 +30 0 +31 0 +32 0 +33 0 +34 0 +35 0 +36 0.5 +37 1 +38 2.5 +39 4 +40 6 +41 8.5 +42 11 +43 124 +44 236.5 +45 239 +46 241 +47 242.5 +48 244 +49 244.5 +50 245 +51 244.5 +52 244 +53 242.5 +54 241 +55 239 +56 236.5 +57 124 +58 11 +59 8.5 +60 6 +61 4 +62 2.5 +63 1 +64 0.5 +65 0 +66 0 +67 0 +68 0 +69 0 +70 0 +71 0 +72 0 +73 0 +74 0 +75 0 +76 0 +77 0 +78 0 +79 0 +80 0 +81 0 +82 0 +83 0 +84 0 +85 0 +86 0 +87 0 +88 0 +89 0 +90 0 +91 0 +92 0 +93 0 +94 0 +95 0 +96 0 +97 0 +98 0 +99 0 diff --git a/Modules/Filtering/ParabolicMorphology/test/images/sharp3.txt b/Modules/Filtering/ParabolicMorphology/test/images/sharp3.txt new file mode 100644 index 00000000000..1f9102c675a --- /dev/null +++ b/Modules/Filtering/ParabolicMorphology/test/images/sharp3.txt @@ -0,0 +1,100 @@ +0 0 +1 0 +2 0 +3 0 +4 0 +5 0 +6 0 +7 0 +8 0 +9 0 +10 0 +11 0 +12 0 +13 0 +14 0 +15 0 +16 0 +17 0 +18 0 +19 0 +20 0 +21 0 +22 0 +23 0 +24 0 +25 0 +26 0 +27 0 +28 0 +29 0 +30 0 +31 0 +32 0 +33 0 +34 0 +35 0 +36 0.5 +37 1 +38 1.5 +39 3 +40 4.5 +41 6 +42 8 +43 237 +44 239 +45 240.5 +46 242 +47 243.5 +48 244 +49 244.5 +50 245 +51 244.5 +52 244 +53 243.5 +54 242 +55 240.5 +56 239 +57 237 +58 8 +59 6 +60 4.5 +61 3 +62 1.5 +63 1 +64 0.5 +65 0 +66 0 +67 0 +68 0 +69 0 +70 0 +71 0 +72 0 +73 0 +74 0 +75 0 +76 0 +77 0 +78 0 +79 0 +80 0 +81 0 +82 0 +83 0 +84 0 +85 0 +86 0 +87 0 +88 0 +89 0 +90 0 +91 0 +92 0 +93 0 +94 0 +95 0 +96 0 +97 0 +98 0 +99 0 diff --git a/Modules/Filtering/ParabolicMorphology/test/images/test.png b/Modules/Filtering/ParabolicMorphology/test/images/test.png new file mode 100644 index 00000000000..11ba5918d9e Binary files /dev/null and b/Modules/Filtering/ParabolicMorphology/test/images/test.png differ diff --git a/Modules/Filtering/ParabolicMorphology/test/itkBinaryCloseParabolicTest.cxx b/Modules/Filtering/ParabolicMorphology/test/itkBinaryCloseParabolicTest.cxx new file mode 100644 index 00000000000..899a37c35b2 --- /dev/null +++ b/Modules/Filtering/ParabolicMorphology/test/itkBinaryCloseParabolicTest.cxx @@ -0,0 +1,88 @@ +/*========================================================================= + * + * Copyright NumFOCUS + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * https://www.apache.org/licenses/LICENSE-2.0.txt + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + *=========================================================================*/ +#include +#include "itkImageFileReader.h" +#include "itkImageFileWriter.h" +#include "itkCommand.h" +#include "itkSimpleFilterWatcher.h" + +#include "itkTimeProbe.h" +#include "itkMultiThreaderBase.h" +#include +#include "itkBinaryCloseParabolicImageFilter.h" + +int +itkBinaryCloseParabolicTest(int argc, char * argv[]) +{ + if (argc != 5) + { + std::cerr << "Usage: " << argv[0] << " inputimage threshold size outim " << std::endl; + return (EXIT_FAILURE); + } + + itk::MultiThreaderBase::SetGlobalMaximumNumberOfThreads(1); + constexpr int dim = 2; + + using PType = unsigned char; + using IType = itk::Image; + + using ReaderType = itk::ImageFileReader; + ReaderType::Pointer reader = ReaderType::New(); + reader->SetFileName(argv[1]); + try + { + reader->Update(); + } + catch (itk::ExceptionObject & excp) + { + std::cerr << excp << std::endl; + return EXIT_FAILURE; + } + // threshold the input to create a mask + using ThreshType = itk::BinaryThresholdImageFilter; + ThreshType::Pointer thresh = ThreshType::New(); + thresh->SetInput(reader->GetOutput()); + + thresh->SetUpperThreshold(std::stoi(argv[2])); + thresh->SetInsideValue(0); + thresh->SetOutsideValue(1); + // now to apply the erosion + using FilterType = itk::BinaryCloseParabolicImageFilter; + + FilterType::Pointer filter = FilterType::New(); + + filter->SetInput(thresh->GetOutput()); + filter->SetUseImageSpacing(true); + filter->SetRadius(std::stod(argv[3])); + + using WriterType = itk::ImageFileWriter; + WriterType::Pointer writer = WriterType::New(); + writer->SetInput(filter->GetOutput()); + writer->SetFileName(argv[4]); + try + { + writer->Update(); + } + catch (itk::ExceptionObject & excp) + { + std::cerr << excp << std::endl; + return EXIT_FAILURE; + } + + return EXIT_SUCCESS; +} diff --git a/Modules/Filtering/ParabolicMorphology/test/itkBinaryDilateParabolicTest.cxx b/Modules/Filtering/ParabolicMorphology/test/itkBinaryDilateParabolicTest.cxx new file mode 100644 index 00000000000..6993cada922 --- /dev/null +++ b/Modules/Filtering/ParabolicMorphology/test/itkBinaryDilateParabolicTest.cxx @@ -0,0 +1,88 @@ +/*========================================================================= + * + * Copyright NumFOCUS + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * https://www.apache.org/licenses/LICENSE-2.0.txt + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + *=========================================================================*/ +#include +#include "itkImageFileReader.h" +#include "itkImageFileWriter.h" +#include "itkCommand.h" +#include "itkSimpleFilterWatcher.h" + +#include "itkTimeProbe.h" +#include "itkMultiThreaderBase.h" +#include +#include "itkBinaryDilateParabolicImageFilter.h" + +int +itkBinaryDilateParabolicTest(int argc, char * argv[]) +{ + if (argc != 5) + { + std::cerr << "Usage: " << argv[0] << " inputimage threshold size outim " << std::endl; + return (EXIT_FAILURE); + } + + itk::MultiThreaderBase::SetGlobalMaximumNumberOfThreads(1); + constexpr int dim = 2; + + using PType = unsigned char; + using IType = itk::Image; + + using ReaderType = itk::ImageFileReader; + ReaderType::Pointer reader = ReaderType::New(); + reader->SetFileName(argv[1]); + try + { + reader->Update(); + } + catch (itk::ExceptionObject & excp) + { + std::cerr << excp << std::endl; + return EXIT_FAILURE; + } + // threshold the input to create a mask + using ThreshType = itk::BinaryThresholdImageFilter; + ThreshType::Pointer thresh = ThreshType::New(); + thresh->SetInput(reader->GetOutput()); + + thresh->SetUpperThreshold(std::stoi(argv[2])); + thresh->SetInsideValue(0); + thresh->SetOutsideValue(1); + // now to apply the erosion + using FilterType = itk::BinaryDilateParabolicImageFilter; + + FilterType::Pointer filter = FilterType::New(); + + filter->SetInput(thresh->GetOutput()); + filter->SetUseImageSpacing(true); + filter->SetRadius(std::stod(argv[3])); + + using WriterType = itk::ImageFileWriter; + WriterType::Pointer writer = WriterType::New(); + writer->SetInput(filter->GetOutput()); + writer->SetFileName(argv[4]); + try + { + writer->Update(); + } + catch (itk::ExceptionObject & excp) + { + std::cerr << excp << std::endl; + return EXIT_FAILURE; + } + + return EXIT_SUCCESS; +} diff --git a/Modules/Filtering/ParabolicMorphology/test/itkBinaryErodeParabolicTest.cxx b/Modules/Filtering/ParabolicMorphology/test/itkBinaryErodeParabolicTest.cxx new file mode 100644 index 00000000000..5da07c94e8a --- /dev/null +++ b/Modules/Filtering/ParabolicMorphology/test/itkBinaryErodeParabolicTest.cxx @@ -0,0 +1,88 @@ +/*========================================================================= + * + * Copyright NumFOCUS + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * https://www.apache.org/licenses/LICENSE-2.0.txt + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + *=========================================================================*/ +#include +#include "itkImageFileReader.h" +#include "itkImageFileWriter.h" +#include "itkCommand.h" +#include "itkSimpleFilterWatcher.h" + +#include "itkTimeProbe.h" +#include "itkMultiThreaderBase.h" +#include +#include "itkBinaryErodeParabolicImageFilter.h" + +int +itkBinaryErodeParabolicTest(int argc, char * argv[]) +{ + if (argc != 5) + { + std::cerr << "Usage: " << argv[0] << " inputimage threshold size outim " << std::endl; + return (EXIT_FAILURE); + } + + itk::MultiThreaderBase::SetGlobalMaximumNumberOfThreads(1); + constexpr int dim = 2; + + using PType = unsigned char; + using IType = itk::Image; + + using ReaderType = itk::ImageFileReader; + ReaderType::Pointer reader = ReaderType::New(); + reader->SetFileName(argv[1]); + try + { + reader->Update(); + } + catch (itk::ExceptionObject & excp) + { + std::cerr << excp << std::endl; + return EXIT_FAILURE; + } + // threshold the input to create a mask + using ThreshType = itk::BinaryThresholdImageFilter; + ThreshType::Pointer thresh = ThreshType::New(); + thresh->SetInput(reader->GetOutput()); + + thresh->SetUpperThreshold(std::stoi(argv[2])); + thresh->SetInsideValue(0); + thresh->SetOutsideValue(1); + // now to apply the erosion + using FilterType = itk::BinaryErodeParabolicImageFilter; + + FilterType::Pointer filter = FilterType::New(); + + filter->SetInput(thresh->GetOutput()); + filter->SetUseImageSpacing(true); + filter->SetRadius(std::stod(argv[3])); + + using WriterType = itk::ImageFileWriter; + WriterType::Pointer writer = WriterType::New(); + writer->SetInput(filter->GetOutput()); + writer->SetFileName(argv[4]); + try + { + writer->Update(); + } + catch (itk::ExceptionObject & excp) + { + std::cerr << excp << std::endl; + return EXIT_FAILURE; + } + + return EXIT_SUCCESS; +} diff --git a/Modules/Filtering/ParabolicMorphology/test/itkBinaryOpenParabolicTest.cxx b/Modules/Filtering/ParabolicMorphology/test/itkBinaryOpenParabolicTest.cxx new file mode 100644 index 00000000000..f07d6771eeb --- /dev/null +++ b/Modules/Filtering/ParabolicMorphology/test/itkBinaryOpenParabolicTest.cxx @@ -0,0 +1,88 @@ +/*========================================================================= + * + * Copyright NumFOCUS + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * https://www.apache.org/licenses/LICENSE-2.0.txt + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + *=========================================================================*/ +#include +#include "itkImageFileReader.h" +#include "itkImageFileWriter.h" +#include "itkCommand.h" +#include "itkSimpleFilterWatcher.h" + +#include "itkTimeProbe.h" +#include "itkMultiThreaderBase.h" +#include +#include "itkBinaryOpenParabolicImageFilter.h" + +int +itkBinaryOpenParabolicTest(int argc, char * argv[]) +{ + if (argc != 5) + { + std::cerr << "Usage: " << argv[0] << " inputimage threshold size outim " << std::endl; + return (EXIT_FAILURE); + } + + itk::MultiThreaderBase::SetGlobalMaximumNumberOfThreads(1); + constexpr int dim = 2; + + using PType = unsigned char; + using IType = itk::Image; + + using ReaderType = itk::ImageFileReader; + ReaderType::Pointer reader = ReaderType::New(); + reader->SetFileName(argv[1]); + try + { + reader->Update(); + } + catch (itk::ExceptionObject & excp) + { + std::cerr << excp << std::endl; + return EXIT_FAILURE; + } + // threshold the input to create a mask + using ThreshType = itk::BinaryThresholdImageFilter; + ThreshType::Pointer thresh = ThreshType::New(); + thresh->SetInput(reader->GetOutput()); + + thresh->SetUpperThreshold(std::stoi(argv[2])); + thresh->SetInsideValue(0); + thresh->SetOutsideValue(1); + // now to apply the erosion + using FilterType = itk::BinaryOpenParabolicImageFilter; + + FilterType::Pointer filter = FilterType::New(); + + filter->SetInput(thresh->GetOutput()); + filter->SetUseImageSpacing(true); + filter->SetRadius(std::stod(argv[3])); + + using WriterType = itk::ImageFileWriter; + WriterType::Pointer writer = WriterType::New(); + writer->SetInput(filter->GetOutput()); + writer->SetFileName(argv[4]); + try + { + writer->Update(); + } + catch (itk::ExceptionObject & excp) + { + std::cerr << excp << std::endl; + return EXIT_FAILURE; + } + + return EXIT_SUCCESS; +} diff --git a/Modules/Filtering/ParabolicMorphology/test/itkParaDTTest.cxx b/Modules/Filtering/ParabolicMorphology/test/itkParaDTTest.cxx new file mode 100644 index 00000000000..8718e6eaa79 --- /dev/null +++ b/Modules/Filtering/ParabolicMorphology/test/itkParaDTTest.cxx @@ -0,0 +1,103 @@ +/*========================================================================= + * + * Copyright NumFOCUS + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * https://www.apache.org/licenses/LICENSE-2.0.txt + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + *=========================================================================*/ +#include +#include "itkImageFileReader.h" +#include "itkImageFileWriter.h" +#include "itkCommand.h" +#include "itkSimpleFilterWatcher.h" + +#include "itkBinaryThresholdImageFilter.h" +#include "itkMorphologicalDistanceTransformImageFilter.h" +#include "itkTimeProbe.h" +#include "itkMultiThreaderBase.h" + +int +itkParaDTTest(int argc, char * argv[]) +{ + // int iterations = 1; + + if (argc != 5) + { + std::cerr << "Usage: " << argv[0] << " inputimage threshold outsideval outim1" << std::endl; + return (EXIT_FAILURE); + } + + itk::MultiThreaderBase::SetGlobalMaximumNumberOfThreads(1); + constexpr int dim = 2; + + using PType = unsigned char; + using IType = itk::Image; + using FType = itk::Image; + + using ReaderType = itk::ImageFileReader; + ReaderType::Pointer reader = ReaderType::New(); + reader->SetFileName(argv[1]); + try + { + reader->Update(); + } + catch (itk::ExceptionObject & excp) + { + std::cerr << excp << std::endl; + return EXIT_FAILURE; + } + + IType::Pointer input = reader->GetOutput(); + + // threshold the input to create a mask + using ThreshType = itk::BinaryThresholdImageFilter; + ThreshType::Pointer thresh = ThreshType::New(); + thresh->SetInput(input); + + thresh->SetUpperThreshold(std::stoi(argv[2])); + thresh->SetInsideValue(0); + thresh->SetOutsideValue(255); + + // now to apply the distance transform + using FilterType = itk::MorphologicalDistanceTransformImageFilter; + + FilterType::Pointer filter = FilterType::New(); + + filter->SetInput(thresh->GetOutput()); + filter->SetOutsideValue(std::stoi(argv[3])); + try + { + filter->Update(); + } + catch (itk::ExceptionObject & excp) + { + std::cerr << excp << std::endl; + return EXIT_FAILURE; + } + + using WriterType = itk::ImageFileWriter; + WriterType::Pointer writer = WriterType::New(); + writer->SetInput(filter->GetOutput()); + writer->SetFileName(argv[4]); + try + { + writer->Update(); + } + catch (itk::ExceptionObject & excp) + { + std::cerr << excp << std::endl; + return EXIT_FAILURE; + } + + return EXIT_SUCCESS; +} diff --git a/Modules/Filtering/ParabolicMorphology/test/itkParaDilateTest.cxx b/Modules/Filtering/ParabolicMorphology/test/itkParaDilateTest.cxx new file mode 100644 index 00000000000..020d0f425e7 --- /dev/null +++ b/Modules/Filtering/ParabolicMorphology/test/itkParaDilateTest.cxx @@ -0,0 +1,110 @@ +/*========================================================================= + * + * Copyright NumFOCUS + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * https://www.apache.org/licenses/LICENSE-2.0.txt + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + *=========================================================================*/ +#include +#include "itkImageFileReader.h" +#include "itkImageFileWriter.h" +#include "itkChangeInformationImageFilter.h" +#include "itkCommand.h" +#include "itkSimpleFilterWatcher.h" + +#include "itkParabolicDilateImageFilter.h" +#include "itkTimeProbe.h" +#include "itkMultiThreaderBase.h" + +// sanity check of the image spacing option + +int +itkParaDilateTest(int argc, char * argv[]) +{ + itk::MultiThreaderBase::SetGlobalMaximumNumberOfThreads(1); + constexpr int dim = 2; + + using PType = unsigned char; + using IType = itk::Image; + + float scale(1.0); + if (argc > 4) + { + scale = std::stod(argv[4]); + } + + using ReaderType = itk::ImageFileReader; + ReaderType::Pointer reader = ReaderType::New(); + reader->SetFileName(argv[1]); + try + { + reader->Update(); + } + catch (itk::ExceptionObject & excp) + { + std::cerr << excp << std::endl; + return EXIT_FAILURE; + } + + using FilterType = itk::ParabolicDilateImageFilter; + + FilterType::Pointer filter = FilterType::New(); + + filter->SetInput(reader->GetOutput()); + + filter->SetScale(scale); + filter->SetUseImageSpacing(true); + filter->SetParabolicAlgorithm(FilterType::INTERSECTION); + try + { + filter->Update(); + } + catch (itk::ExceptionObject & excp) + { + std::cerr << excp << std::endl; + return EXIT_FAILURE; + } + + using WriterType = itk::ImageFileWriter; + WriterType::Pointer writer = WriterType::New(); + writer->SetInput(filter->GetOutput()); + writer->SetFileName(argv[2]); + try + { + writer->Update(); + } + catch (itk::ExceptionObject & excp) + { + std::cerr << excp << std::endl; + return EXIT_FAILURE; + } + + filter->SetScale(scale); + filter->SetUseImageSpacing(true); + filter->SetParabolicAlgorithm(FilterType::CONTACTPOINT); + filter->Update(); + + writer->SetInput(filter->GetOutput()); + writer->SetFileName(argv[3]); + try + { + writer->Update(); + } + catch (itk::ExceptionObject & excp) + { + std::cerr << excp << std::endl; + return EXIT_FAILURE; + } + + return EXIT_SUCCESS; +} diff --git a/Modules/Filtering/ParabolicMorphology/test/itkParaErodeTest.cxx b/Modules/Filtering/ParabolicMorphology/test/itkParaErodeTest.cxx new file mode 100644 index 00000000000..01006387c41 --- /dev/null +++ b/Modules/Filtering/ParabolicMorphology/test/itkParaErodeTest.cxx @@ -0,0 +1,110 @@ +/*========================================================================= + * + * Copyright NumFOCUS + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * https://www.apache.org/licenses/LICENSE-2.0.txt + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + *=========================================================================*/ +#include +#include "itkImageFileReader.h" +#include "itkImageFileWriter.h" +#include "itkChangeInformationImageFilter.h" +#include "itkCommand.h" +#include "itkSimpleFilterWatcher.h" + +#include "itkParabolicErodeImageFilter.h" +#include "itkTimeProbe.h" +#include "itkMultiThreaderBase.h" + +// sanity check of the image spacing option + +int +itkParaErodeTest(int argc, char * argv[]) +{ + itk::MultiThreaderBase::SetGlobalMaximumNumberOfThreads(1); + constexpr int dim = 2; + + using PType = unsigned char; + using IType = itk::Image; + + float scale(1.0); + if (argc > 4) + { + scale = std::stod(argv[4]); + } + + using ReaderType = itk::ImageFileReader; + ReaderType::Pointer reader = ReaderType::New(); + reader->SetFileName(argv[1]); + try + { + reader->Update(); + } + catch (itk::ExceptionObject & excp) + { + std::cerr << excp << std::endl; + return EXIT_FAILURE; + } + + using FilterType = itk::ParabolicErodeImageFilter; + + FilterType::Pointer filter = FilterType::New(); + + filter->SetInput(reader->GetOutput()); + + filter->SetScale(scale); + filter->SetUseImageSpacing(true); + filter->SetParabolicAlgorithm(FilterType::INTERSECTION); + try + { + filter->Update(); + } + catch (itk::ExceptionObject & excp) + { + std::cerr << excp << std::endl; + return EXIT_FAILURE; + } + + using WriterType = itk::ImageFileWriter; + WriterType::Pointer writer = WriterType::New(); + writer->SetInput(filter->GetOutput()); + writer->SetFileName(argv[2]); + try + { + writer->Update(); + } + catch (itk::ExceptionObject & excp) + { + std::cerr << excp << std::endl; + return EXIT_FAILURE; + } + + filter->SetScale(scale); + filter->SetUseImageSpacing(true); + filter->SetParabolicAlgorithm(FilterType::CONTACTPOINT); + filter->Update(); + + writer->SetInput(filter->GetOutput()); + writer->SetFileName(argv[3]); + try + { + writer->Update(); + } + catch (itk::ExceptionObject & excp) + { + std::cerr << excp << std::endl; + return EXIT_FAILURE; + } + + return EXIT_SUCCESS; +} diff --git a/Modules/Filtering/ParabolicMorphology/test/itkParaOpenTest.cxx b/Modules/Filtering/ParabolicMorphology/test/itkParaOpenTest.cxx new file mode 100644 index 00000000000..ea5afdc2589 --- /dev/null +++ b/Modules/Filtering/ParabolicMorphology/test/itkParaOpenTest.cxx @@ -0,0 +1,112 @@ +/*========================================================================= + * + * Copyright NumFOCUS + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * https://www.apache.org/licenses/LICENSE-2.0.txt + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + *=========================================================================*/ +#include +#include "itkImageFileReader.h" +#include "itkImageFileWriter.h" +#include "itkCommand.h" +#include "itkSimpleFilterWatcher.h" + +#include "itkParabolicOpenImageFilter.h" +#include "itkTimeProbe.h" +#include "itkMultiThreaderBase.h" + +int +itkParaOpenTest(int argc, char * argv[]) +{ + itk::MultiThreaderBase::SetGlobalMaximumNumberOfThreads(1); + constexpr int dim = 2; + + using PType = unsigned char; + using IType = itk::Image; + + using ReaderType = itk::ImageFileReader; + ReaderType::Pointer reader = ReaderType::New(); + reader->SetFileName(argv[1]); + try + { + reader->Update(); + } + catch (itk::ExceptionObject & excp) + { + std::cerr << excp << std::endl; + return EXIT_FAILURE; + } + + using FilterType = itk::ParabolicOpenImageFilter; + + FilterType::Pointer filter = FilterType::New(); + + filter->SetInput(reader->GetOutput()); + filter->SetSafeBorder(true); + FilterType::RadiusType scale; + scale[0] = 1; + scale[1] = 0.5; + filter->SetScale(scale); + + filter->SetParabolicAlgorithm(FilterType::INTERSECTION); + try + { + filter->Update(); + } + catch (itk::ExceptionObject & excp) + { + std::cerr << excp << std::endl; + return EXIT_FAILURE; + } + + using WriterType = itk::ImageFileWriter; + WriterType::Pointer writer = WriterType::New(); + writer->SetInput(filter->GetOutput()); + writer->SetFileName(argv[2]); + try + { + writer->Update(); + } + catch (itk::ExceptionObject & excp) + { + std::cerr << excp << std::endl; + return EXIT_FAILURE; + } + + if (argc > 3) + { + // testing equivalence + filter->SetParabolicAlgorithm(FilterType::CONTACTPOINT); + try + { + filter->Update(); + } + catch (itk::ExceptionObject & excp) + { + std::cerr << excp << std::endl; + return EXIT_FAILURE; + } + writer->SetFileName(argv[3]); + try + { + writer->Update(); + } + catch (itk::ExceptionObject & excp) + { + std::cerr << excp << std::endl; + return EXIT_FAILURE; + } + } + + return EXIT_SUCCESS; +} diff --git a/Modules/Filtering/ParabolicMorphology/test/itkParaSharpenTest.cxx b/Modules/Filtering/ParabolicMorphology/test/itkParaSharpenTest.cxx new file mode 100644 index 00000000000..71150722958 --- /dev/null +++ b/Modules/Filtering/ParabolicMorphology/test/itkParaSharpenTest.cxx @@ -0,0 +1,170 @@ +/*========================================================================= + * + * Copyright NumFOCUS + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * https://www.apache.org/licenses/LICENSE-2.0.txt + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + *=========================================================================*/ +#include +#include "itkImageFileReader.h" +#include "itkImageFileWriter.h" +#include "itkCommand.h" +#include "itkSimpleFilterWatcher.h" + +// #include "plotutils.h" + +#include "itkGrayscaleDilateImageFilter.h" +#include "itkBinaryBallStructuringElement.h" + +#include "itkSmoothingRecursiveGaussianImageFilter.h" +#include "itkBinaryThresholdImageFilter.h" +#include "itkMorphologicalSharpeningImageFilter.h" +#include "itkTimeProbe.h" +#include "itkMultiThreaderBase.h" + +int +itkParaSharpenTest(int argc, char * argv[]) +{ + int iterations = 1; + + if (argc != 3) + { + std::cerr << "Usage: " << argv[0] << " iterations outputimage " << std::endl; + return (EXIT_FAILURE); + } + + iterations = std::stoi(argv[1]); + + // itk::MultiThreaderBase::SetGlobalMaximumNumberOfThreads( 1 ); + constexpr int dim = 2; + + using PType = unsigned char; + using IType = itk::Image; + using FType = itk::Image; + + using WriterType = itk::ImageFileWriter; + WriterType::Pointer writer = WriterType::New(); + + // create the input image - we will blur a dot, threshold then blur again + IType::Pointer input = IType::New(); + IType::SizeType size; + IType::IndexType index; + IType::RegionType region; + IType::SpacingType spacing; + + size.Fill(100); + spacing.Fill(1); + region.SetSize(size); + + index.Fill(50); + + input->SetRegions(region); + input->SetSpacing(spacing); + input->Allocate(); + + input->FillBuffer(0); + input->SetPixel(index, 255); + + using SRType = itk::BinaryBallStructuringElement; + using DilateType = itk::GrayscaleDilateImageFilter; + DilateType::Pointer smallDilate = DilateType::New(); + SRType smallkernel; + SRType::RadiusType smallrad = smallkernel.GetRadius(); + smallrad.Fill(7); + smallkernel.SetRadius(smallrad); + smallkernel.CreateStructuringElement(); + smallDilate->SetKernel(smallkernel); + + using SmootherType = itk::SmoothingRecursiveGaussianImageFilter; + + SmootherType::Pointer smoother = SmootherType::New(); + + smallDilate->SetInput(input); + + smoother->SetInput(smallDilate->GetOutput()); + smoother->SetSigma(3); + + writer->SetInput(smallDilate->GetOutput()); + writer->SetFileName("input.tif"); + try + { + writer->Update(); + } + catch (itk::ExceptionObject & excp) + { + std::cerr << excp << std::endl; + return EXIT_FAILURE; + } + + writer->SetInput(smoother->GetOutput()); + writer->SetFileName("blurrredinput.tif"); + try + { + writer->Update(); + } + catch (itk::ExceptionObject & excp) + { + std::cerr << excp << std::endl; + return EXIT_FAILURE; + } + + // now to apply the sharpening + + using FilterType = itk::MorphologicalSharpeningImageFilter; + + FilterType::Pointer filter = FilterType::New(); + + filter->SetInput(smoother->GetOutput()); + filter->SetScale(1); + filter->SetIterations(iterations); + try + { + filter->Update(); + } + catch (itk::ExceptionObject & excp) + { + std::cerr << excp << std::endl; + return EXIT_FAILURE; + } + + using FlWriterType = itk::ImageFileWriter; + FlWriterType::Pointer flwriter = FlWriterType::New(); + + flwriter->SetInput(filter->GetOutput()); + flwriter->SetFileName(argv[2]); + try + { + flwriter->Update(); + } + catch (itk::ExceptionObject & excp) + { + std::cerr << excp << std::endl; + return EXIT_FAILURE; + } + +#if 0 + // this stuff as for the article + // write out profiles + // input + IType::IndexType first, last; + first[0] = 50; + first[1] = 0; + last[0] = 50; + last[1] = 99; + + extractProfile< IType >(smallDilate->GetOutput(), first, last, "inputprof.txt"); + extractProfile< IType >(smoother->GetOutput(), first, last, "blurredprof.txt"); + extractProfile< FType >(filter->GetOutput(), first, last, argv[3]); +#endif + return EXIT_SUCCESS; +} diff --git a/Modules/Filtering/ParabolicMorphology/test/itkParaSpacingTest.cxx b/Modules/Filtering/ParabolicMorphology/test/itkParaSpacingTest.cxx new file mode 100644 index 00000000000..7a52d35d709 --- /dev/null +++ b/Modules/Filtering/ParabolicMorphology/test/itkParaSpacingTest.cxx @@ -0,0 +1,135 @@ +/*========================================================================= + * + * Copyright NumFOCUS + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * https://www.apache.org/licenses/LICENSE-2.0.txt + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + *=========================================================================*/ +#include +#include "itkImageFileReader.h" +#include "itkImageFileWriter.h" +#include "itkChangeInformationImageFilter.h" +#include "itkCommand.h" +#include "itkSimpleFilterWatcher.h" + +#include "itkParabolicOpenImageFilter.h" +#include "itkTimeProbe.h" +#include "itkMultiThreaderBase.h" + +// sanity check of the image spacing option + +int +itkParaSpacingTest(int, char * argv[]) +{ + itk::MultiThreaderBase::SetGlobalMaximumNumberOfThreads(1); + constexpr int dim = 2; + + using PType = unsigned char; + using IType = itk::Image; + + using ReaderType = itk::ImageFileReader; + ReaderType::Pointer reader = ReaderType::New(); + reader->SetFileName(argv[1]); + try + { + reader->Update(); + } + catch (itk::ExceptionObject & excp) + { + std::cerr << excp << std::endl; + return EXIT_FAILURE; + } + + using FilterType = itk::ParabolicOpenImageFilter; + + FilterType::Pointer filter = FilterType::New(); + + filter->SetInput(reader->GetOutput()); + filter->SetSafeBorder(true); + FilterType::RadiusType scale; + scale[0] = 1; + scale[1] = 0.5; + + filter->SetScale(scale); + filter->SetUseImageSpacing(false); + try + { + filter->Update(); + } + catch (itk::ExceptionObject & excp) + { + std::cerr << excp << std::endl; + return EXIT_FAILURE; + } + + using WriterType = itk::ImageFileWriter; + WriterType::Pointer writer = WriterType::New(); + writer->SetInput(filter->GetOutput()); + writer->SetFileName(argv[2]); + try + { + writer->Update(); + } + catch (itk::ExceptionObject & excp) + { + std::cerr << excp << std::endl; + return EXIT_FAILURE; + } + + // now we'll change the image spacing and see if we can reproduce + // the result + using ChangeType = itk::ChangeInformationImageFilter; + ChangeType::Pointer changer = ChangeType::New(); + changer->SetInput(reader->GetOutput()); + ChangeType::SpacingType oldspacing, newspacing; + + oldspacing = filter->GetOutput()->GetSpacing(); + newspacing[0] = 1 / sqrt((float)1); + newspacing[1] = 1 / sqrt((float)0.5); + + changer->SetOutputSpacing(newspacing); + changer->ChangeSpacingOn(); + // set scales to deliver the same result + scale[0] = 1; + scale[1] = 1; + filter->SetInput(changer->GetOutput()); + filter->SetScale(scale); + filter->SetUseImageSpacing(true); + // change the spacing back to original to allow comparison + ChangeType::Pointer changerback = ChangeType::New(); + changerback->SetInput(filter->GetOutput()); + changerback->SetOutputSpacing(oldspacing); + changerback->ChangeSpacingOn(); + + try + { + changerback->Update(); + } + catch (itk::ExceptionObject & excp) + { + std::cerr << excp << std::endl; + return EXIT_FAILURE; + } + writer->SetInput(changerback->GetOutput()); + writer->SetFileName(argv[3]); + try + { + writer->Update(); + } + catch (itk::ExceptionObject & excp) + { + std::cerr << excp << std::endl; + return EXIT_FAILURE; + } + return EXIT_SUCCESS; +} diff --git a/Modules/Filtering/ParabolicMorphology/wrapping/CMakeLists.txt b/Modules/Filtering/ParabolicMorphology/wrapping/CMakeLists.txt new file mode 100644 index 00000000000..c110d759119 --- /dev/null +++ b/Modules/Filtering/ParabolicMorphology/wrapping/CMakeLists.txt @@ -0,0 +1,3 @@ +itk_wrap_module(ParabolicMorphology) +itk_auto_load_submodules() +itk_end_wrap_module() diff --git a/Modules/Filtering/ParabolicMorphology/wrapping/itkBinaryCloseParabolicImageFilter.wrap b/Modules/Filtering/ParabolicMorphology/wrapping/itkBinaryCloseParabolicImageFilter.wrap new file mode 100644 index 00000000000..4c9152b989a --- /dev/null +++ b/Modules/Filtering/ParabolicMorphology/wrapping/itkBinaryCloseParabolicImageFilter.wrap @@ -0,0 +1,3 @@ +itk_wrap_class("itk::BinaryCloseParabolicImageFilter" POINTER_WITH_SUPERCLASS) +itk_wrap_image_filter("${WRAP_ITK_SCALAR}" 2 2+) +itk_end_wrap_class() diff --git a/Modules/Filtering/ParabolicMorphology/wrapping/itkBinaryDilateParabolicImageFilter.wrap b/Modules/Filtering/ParabolicMorphology/wrapping/itkBinaryDilateParabolicImageFilter.wrap new file mode 100644 index 00000000000..8e248d520a8 --- /dev/null +++ b/Modules/Filtering/ParabolicMorphology/wrapping/itkBinaryDilateParabolicImageFilter.wrap @@ -0,0 +1,3 @@ +itk_wrap_class("itk::BinaryDilateParabolicImageFilter" POINTER_WITH_SUPERCLASS) +itk_wrap_image_filter("${WRAP_ITK_SCALAR}" 2 2+) +itk_end_wrap_class() diff --git a/Modules/Filtering/ParabolicMorphology/wrapping/itkBinaryErodeParabolicImageFilter.wrap b/Modules/Filtering/ParabolicMorphology/wrapping/itkBinaryErodeParabolicImageFilter.wrap new file mode 100644 index 00000000000..486aabd9631 --- /dev/null +++ b/Modules/Filtering/ParabolicMorphology/wrapping/itkBinaryErodeParabolicImageFilter.wrap @@ -0,0 +1,3 @@ +itk_wrap_class("itk::BinaryErodeParabolicImageFilter" POINTER_WITH_SUPERCLASS) +itk_wrap_image_filter("${WRAP_ITK_SCALAR}" 2 2+) +itk_end_wrap_class() diff --git a/Modules/Filtering/ParabolicMorphology/wrapping/itkBinaryOpenParabolicImageFilter.wrap b/Modules/Filtering/ParabolicMorphology/wrapping/itkBinaryOpenParabolicImageFilter.wrap new file mode 100644 index 00000000000..65a347e2b08 --- /dev/null +++ b/Modules/Filtering/ParabolicMorphology/wrapping/itkBinaryOpenParabolicImageFilter.wrap @@ -0,0 +1,3 @@ +itk_wrap_class("itk::BinaryOpenParabolicImageFilter" POINTER_WITH_SUPERCLASS) +itk_wrap_image_filter("${WRAP_ITK_SCALAR}" 2 2+) +itk_end_wrap_class() diff --git a/Modules/Filtering/ParabolicMorphology/wrapping/itkMorphologicalDistanceTransformImageFilter.wrap b/Modules/Filtering/ParabolicMorphology/wrapping/itkMorphologicalDistanceTransformImageFilter.wrap new file mode 100644 index 00000000000..e8a46c09f46 --- /dev/null +++ b/Modules/Filtering/ParabolicMorphology/wrapping/itkMorphologicalDistanceTransformImageFilter.wrap @@ -0,0 +1,3 @@ +itk_wrap_class("itk::MorphologicalDistanceTransformImageFilter" POINTER) +itk_wrap_image_filter_combinations("${WRAP_ITK_USIGN_INT}" "${WRAP_ITK_REAL}") +itk_end_wrap_class() diff --git a/Modules/Filtering/ParabolicMorphology/wrapping/itkMorphologicalSharpeningImageFilter.wrap b/Modules/Filtering/ParabolicMorphology/wrapping/itkMorphologicalSharpeningImageFilter.wrap new file mode 100644 index 00000000000..e3f08df8218 --- /dev/null +++ b/Modules/Filtering/ParabolicMorphology/wrapping/itkMorphologicalSharpeningImageFilter.wrap @@ -0,0 +1,3 @@ +itk_wrap_class("itk::MorphologicalSharpeningImageFilter" POINTER) +itk_wrap_image_filter("${WRAP_ITK_SCALAR}" 2 2+) +itk_end_wrap_class() diff --git a/Modules/Filtering/ParabolicMorphology/wrapping/itkMorphologicalSignedDistanceTransformImageFilter.wrap b/Modules/Filtering/ParabolicMorphology/wrapping/itkMorphologicalSignedDistanceTransformImageFilter.wrap new file mode 100644 index 00000000000..9d3c1363e96 --- /dev/null +++ b/Modules/Filtering/ParabolicMorphology/wrapping/itkMorphologicalSignedDistanceTransformImageFilter.wrap @@ -0,0 +1,3 @@ +itk_wrap_class("itk::MorphologicalSignedDistanceTransformImageFilter" POINTER) +itk_wrap_image_filter_combinations("${WRAP_ITK_USIGN_INT}" "${WRAP_ITK_REAL}") +itk_end_wrap_class() diff --git a/Modules/Filtering/ParabolicMorphology/wrapping/itkParabolicCloseImageFilter.wrap b/Modules/Filtering/ParabolicMorphology/wrapping/itkParabolicCloseImageFilter.wrap new file mode 100644 index 00000000000..04583e77338 --- /dev/null +++ b/Modules/Filtering/ParabolicMorphology/wrapping/itkParabolicCloseImageFilter.wrap @@ -0,0 +1,3 @@ +itk_wrap_class("itk::ParabolicCloseImageFilter" POINTER_WITH_SUPERCLASS) +itk_wrap_image_filter("${WRAP_ITK_SCALAR}" 2 2+) +itk_end_wrap_class() diff --git a/Modules/Filtering/ParabolicMorphology/wrapping/itkParabolicDilateImageFilter.wrap b/Modules/Filtering/ParabolicMorphology/wrapping/itkParabolicDilateImageFilter.wrap new file mode 100644 index 00000000000..0dd82281082 --- /dev/null +++ b/Modules/Filtering/ParabolicMorphology/wrapping/itkParabolicDilateImageFilter.wrap @@ -0,0 +1,3 @@ +itk_wrap_class("itk::ParabolicDilateImageFilter" POINTER_WITH_SUPERCLASS) +itk_wrap_image_filter("${WRAP_ITK_SCALAR}" 2 2+) +itk_end_wrap_class() diff --git a/Modules/Filtering/ParabolicMorphology/wrapping/itkParabolicErodeImageFilter.wrap b/Modules/Filtering/ParabolicMorphology/wrapping/itkParabolicErodeImageFilter.wrap new file mode 100644 index 00000000000..41b67f4eda4 --- /dev/null +++ b/Modules/Filtering/ParabolicMorphology/wrapping/itkParabolicErodeImageFilter.wrap @@ -0,0 +1,3 @@ +itk_wrap_class("itk::ParabolicErodeImageFilter" POINTER_WITH_SUPERCLASS) +itk_wrap_image_filter("${WRAP_ITK_SCALAR}" 2 2+) +itk_end_wrap_class() diff --git a/Modules/Filtering/ParabolicMorphology/wrapping/itkParabolicOpenImageFilter.wrap b/Modules/Filtering/ParabolicMorphology/wrapping/itkParabolicOpenImageFilter.wrap new file mode 100644 index 00000000000..c8f763a1201 --- /dev/null +++ b/Modules/Filtering/ParabolicMorphology/wrapping/itkParabolicOpenImageFilter.wrap @@ -0,0 +1,3 @@ +itk_wrap_class("itk::ParabolicOpenImageFilter" POINTER_WITH_SUPERCLASS) +itk_wrap_image_filter("${WRAP_ITK_SCALAR}" 2 2+) +itk_end_wrap_class() diff --git a/Modules/Remote/ParabolicMorphology.remote.cmake b/Modules/Remote/ParabolicMorphology.remote.cmake deleted file mode 100644 index 5e92848c4e0..00000000000 --- a/Modules/Remote/ParabolicMorphology.remote.cmake +++ /dev/null @@ -1,53 +0,0 @@ -#-- # Grading Level Criteria Report -#-- EVALUATION DATE: 2020-03-01 -#-- EVALUATORS: [<>,<>] -#-- -#-- ## Compliance level 5 star (AKA ITK main modules, or remote modules that could become core modules) -#-- - [ ] Widespread community dependance -#-- - [ ] Above 90% code coverage -#-- - [ ] CI dashboards and testing monitored rigorously -#-- - [ ] Key API features are exposed in wrapping interface -#-- - [ ] All requirements of Levels 4,3,2,1 -#-- -#-- ## Compliance Level 4 star (Very high-quality code, perhaps small community dependance) -#-- - [ ] Meets all ITK code style standards -#-- - [ ] No external requirements beyond those needed by ITK proper -#-- - [ ] Builds and passes tests on all supported platforms within 1 month of each core tagged release -#-- - [ ] Windows Shared Library Build with Visual Studio -#-- - [ ] Mac with clang compiller -#-- - [ ] Linux with gcc compiler -#-- - [ ] Active developer community dedicated to maintaining code-base -#-- - [ ] 75% code coverage demonstrated for testing suite -#-- - [ ] Continuous integration testing performed -#-- - [ ] All requirements of Levels 3,2,1 -#-- -#-- ## Compliance Level 3 star (Quality beta code) -#-- - [ ] API | executable interface is considered mostly stable and feature complete -#-- - [ ] 10% C0-code coverage demonstrated for testing suite -#-- - [ ] Some tests exist and pass on at least some platform -#-- - [X] All requirements of Levels 2,1 -#-- -#-- ## Compliance Level 2 star (Alpha code feature API development or niche community/execution environment dependance ) -#-- - [X] Compiles for at least 1 niche set of execution envirionments, and perhaps others -#-- (may depend on specific external tools like a java environment, or specific external libraries to work ) -#-- - [X] All requirements of Levels 1 -#-- -#-- ## Compliance Level 1 star (Pre-alpha features under development and code of unknown quality) -#-- - [X] Code complies on at least 1 platform -#-- -#-- ## Compliance Level 0 star ( Code/Feature of known poor-quality or deprecated status ) -#-- - [ ] Code reviewed and explicitly identified as not recommended for use -#-- -#-- ### Please document here any justification for the criteria above -# Code style enforced by clang-format on 2020-02-19, and clang-tidy modernizations completed - -itk_fetch_module( - ParabolicMorphology - "Classes performing morphology using parabolic functions. - Fast distance transforms and binary erosions/dilations/openings/closings - by spheres, sharpenings and grayscale operations. - https://doi.org/10.54294/aq68pt" - MODULE_COMPLIANCE_LEVEL 2 - GIT_REPOSITORY https://github.com/InsightSoftwareConsortium/ITKParabolicMorphology.git - GIT_TAG 68286edc7459072eb33422a7d53c855318496922 - ) diff --git a/pyproject.toml b/pyproject.toml index e0770d91a15..a6f86c46f2f 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -69,6 +69,7 @@ cmd = '''cmake \ -DModule_IOMeshMZ3:BOOL=ON \ -DModule_IOFDF:BOOL=ON \ -DModule_MorphologicalContourInterpolation:BOOL=ON \ + -DModule_ParabolicMorphology:BOOL=ON \ -DModule_MultipleImageIterator:BOOL=ON \ -DModule_RLEImage:BOOL=ON \ -DModule_SubdivisionQuadEdgeMeshFilter:BOOL=ON \