diff --git a/Modules/Filtering/AnisotropicDiffusionLBR/CMakeLists.txt b/Modules/Filtering/AnisotropicDiffusionLBR/CMakeLists.txt new file mode 100644 index 00000000000..b746377e0ab --- /dev/null +++ b/Modules/Filtering/AnisotropicDiffusionLBR/CMakeLists.txt @@ -0,0 +1,26 @@ +cmake_minimum_required(VERSION 3.16.3) + +project( + AnisotropicDiffusionLBR + VERSION + 5.1.0 #Version should track with ITK + LANGUAGES + CXX +) + +# Header only module, no libraries set(AnisotropicDiffusionLBR_LIBRARIES AnisotropicDiffusionLBR) + +if(NOT ITK_SOURCE_DIR) + find_package(ITK 5.1 REQUIRED) + list(APPEND CMAKE_MODULE_PATH ${ITK_CMAKE_DIR}) + include(ITKModuleExternal) +else() + set(ITK_DIR ${CMAKE_BINARY_DIR}) + itk_module_impl() +endif() + +# itk_module_examples() is intentionally omitted: the upstream +# examples/ directory is not part of the whitelist-based ingest +# (see Modules/Filtering/AnisotropicDiffusionLBR/README.md), so +# expanding this call under BUILD_EXAMPLES=ON would fail at CMake +# configure time looking for a non-existent subdirectory. diff --git a/Modules/Filtering/AnisotropicDiffusionLBR/README.md b/Modules/Filtering/AnisotropicDiffusionLBR/README.md new file mode 100644 index 00000000000..d0710dbd3df --- /dev/null +++ b/Modules/Filtering/AnisotropicDiffusionLBR/README.md @@ -0,0 +1,60 @@ +# AnisotropicDiffusionLBR + +In-tree ITK module providing coherence-enhancing (CED) and +edge-enhancing (EED) anisotropic diffusion filters built on the +Lattice-Basis-Reduction (LBR) stencil scheme by Jean-Marie Mirebeau. + +## Origin + +Ingested from the standalone remote module +[**InsightSoftwareConsortium/ITKAnisotropicDiffusionLBR**](https://github.com/InsightSoftwareConsortium/ITKAnisotropicDiffusionLBR) +on 2026-04-22, at upstream commit +[`203260b9`](https://github.com/InsightSoftwareConsortium/ITKAnisotropicDiffusionLBR/commit/203260b929acda68a1f64b1267b0d89e825904ec). +The upstream repository will be archived read-only after this PR +merges; it remains reachable at the URL above. + +## What lives here + +Per the v3 ingestion strategy (see +`Utilities/Maintenance/RemoteModuleIngest/INGESTION_STRATEGY.md`), only +paths matching the narrow whitelist (code, headers, tests, wrapping, +module CMake) crossed the merge boundary: + +- `include/` — public C++ headers. +- `test/` — Google-test / CTest drivers and content-link stubs. +- `wrapping/` — Python wrapping descriptors. +- `CMakeLists.txt`, `itk-module.cmake` — build + module descriptors. + +Every surviving commit preserves original authorship; `git blame` +walks across the merge boundary to upstream authors from 2014. + +## What was intentionally left upstream + +Everything outside the whitelist stays in the archived upstream repo. +If you need any of it, clone +. + +| Content in upstream | Why it did not ingest | +|---|---| +| `Old/` | Pre-refactor legacy trees; no current test depends on them. | +| `examples/` + `examples/Data/` | Per-module `examples/` are routed to top-level [`Examples/`](https://github.com/InsightSoftwareConsortium/ITK/tree/main/Examples) via a separate follow-up PR. | +| `README.rst` | Algorithm citations + background are folded into the Doxygen on the filter headers. | +| `test/azure-pipelines.yml`, `test/Docker/` | Standalone-build CI scaffolding; not useful in-tree. | +| `CTestConfig.cmake`, `pyproject.toml`, `.github/`, `.clang-format`, `LICENSE` | Packaging / CI / style scaffolding superseded by ITK root. | + +## Long-form documentation + +- **Algorithm description** — see the Doxygen on + `itkAnisotropicDiffusionLBRImageFilter` and + `itkCoherenceEnhancingDiffusionImageFilter`, plus embedded citations + (Mirebeau et al.). +- **Standalone build + usage examples** — see the archived upstream at + . + +## Content-link status + +The 24 baseline and input content-links under `test/Baseline/` and +`test/Input/` are currently `.md5` stubs. Per the v3 strategy they +should be normalized to `.cid` (IPFS Content Identifier) before this +PR merges; see +`Utilities/Maintenance/RemoteModuleIngest/cid-normalize.sh`. diff --git a/Modules/Filtering/AnisotropicDiffusionLBR/include/CoherenceEnhancingDiffusionCommandLine.h b/Modules/Filtering/AnisotropicDiffusionLBR/include/CoherenceEnhancingDiffusionCommandLine.h new file mode 100644 index 00000000000..0f44efdb2ff --- /dev/null +++ b/Modules/Filtering/AnisotropicDiffusionLBR/include/CoherenceEnhancingDiffusionCommandLine.h @@ -0,0 +1,262 @@ +/*========================================================================= + * + * 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. + * + *=========================================================================*/ +// +// +// Created by Jean-Marie Mirebeau on 20/11/2014. +// +// + +#ifndef itkDiffusion_CoherenceEnhancingDiffusionCommandLine_h +#define itkDiffusion_CoherenceEnhancingDiffusionCommandLine_h + +#include "itkImageFileReader.h" +#include "itkImageFileWriter.h" +#include "itkCoherenceEnhancingDiffusionImageFilter.h" +#include "LinearAnisotropicDiffusionCommandLine.h" +#include "itkTimeProbe.h" + +namespace CoherenceEnhancingDiffusionCommandLine +{ + +inline void +Usage() +{ + std::cerr << "Input image filename. 2D and 3D images supported. Required.\n" + << "Output image filename. Required.\n" + << "Diffusion time. Suggested range: 0.5-5, up to 50 for strong noise. Default: 2.\n" + << "Lambda. Small values detect more edges. Suggested range: 0.05, 0.0001. Default: 0.01\n" + << "Weickert diffusion type. Accepted values: CED, cCED, EED, cEED, Isotropic. Default: CED.\n" + << "Noise scale. Suggested range: 0.5 - 4. Default 1.\n" + << "Feature scale. Suggested range: 2-6. Default 2.\n" + << "\n"; +} + +using namespace itk; + +int +Execute(int argc, char * argv[]); + +template +int +Execute(int argc, char * argv[], itk::ImageIOBase::IOComponentEnum, int nComponents); + +template +int +Execute(int argc, char * argv[], int nComponents); + +template +int +Execute(int argc, char * argv[]); + + +using ReportProgressToCOutType = LinearAnisotropicDiffusionCommandLine::ReportProgressToCOutType; + +inline int +Execute(int argc, char * argv[]) +{ + using std::cerr; + using std::endl; + using namespace itk; + + if (argc < 2 + 1) + { + Usage(); + return EXIT_SUCCESS; + } + + const char * imageFileName = argv[1]; + + itk::ImageIOBase::Pointer imageIO = itk::ImageIOFactory::CreateImageIO(imageFileName, itk::IOFileModeEnum::ReadMode); + if (imageIO.IsNull()) + { + std::cerr << "Could not create ImageIO" << std::endl; + return EXIT_FAILURE; + } + imageIO->SetFileName(imageFileName); + imageIO->ReadImageInformation(); + + const unsigned int imageDimension = imageIO->GetNumberOfDimensions(); + const itk::ImageIOBase::IOComponentEnum componentType = imageIO->GetComponentType(); + const unsigned int nComponents = imageIO->GetNumberOfComponents(); + + switch (imageDimension) + { + case 2: + return Execute<2>(argc, argv, componentType, nComponents); + case 3: + return Execute<3>(argc, argv, componentType, nComponents); + default: + itkGenericExceptionMacro("Sorry, unsupported image dimension."); + } +} + +template +int +Execute(int argc, char * argv[], itk::IOComponentEnum componentType, int nComponents) +{ + switch (componentType) + { + case itk::IOComponentEnum::UCHAR: + return Execute(argc, argv, nComponents); + case itk::IOComponentEnum::FLOAT: + return Execute(argc, argv, nComponents); + case itk::IOComponentEnum::DOUBLE: + return Execute(argc, argv, nComponents); + default: + itkGenericExceptionMacro("Sorry, unsupported component type"); + } +} + +template +int +Execute(int argc, char * argv[], int nComponents) +{ + switch (nComponents) + { + case 1: + return Execute(argc, argv); + case 2: + return Execute, Vector>(argc, argv); + case 3: + return Execute, Vector>(argc, argv); + default: + itkGenericExceptionMacro("Sorry, unsupported number of components"); + } +} + +template +int +Execute(int argc, char * argv[]) +{ + using ImageType = Image; + + using ReaderType = ImageFileReader; + typename ReaderType::Pointer reader = ReaderType::New(); + + const char * imageFileName = argv[1]; + const char * outputFileName = argv[2]; + reader->SetFileName(imageFileName); + + using DiffusionFilterType = CoherenceEnhancingDiffusionImageFilter; + typename DiffusionFilterType::Pointer diffusionFilter = DiffusionFilterType::New(); + diffusionFilter->SetInput(reader->GetOutput()); + + ReportProgressToCOutType::Pointer reportDiffusionProgress = ReportProgressToCOutType::New(); + diffusionFilter->AddObserver(ProgressEvent(), reportDiffusionProgress); + + int argIndex = 3; + if (argIndex < argc) + { + const double diffusionTime = std::stod(argv[argIndex++]); + if (diffusionTime == 0) + itkGenericExceptionMacro("Error: Unrecognized diffusion time (third argument).\n"); + diffusionFilter->SetDiffusionTime(diffusionTime); + } + + if (argIndex < argc) + { + const double lambda = std::stod(argv[argIndex++]); + if (lambda == 0.) + itkGenericExceptionMacro("Error: Unrecognized lambda (fourth argument).\n"); + diffusionFilter->SetLambda(lambda); + } + + if (argIndex < argc) + { + const char * enhancement = argv[argIndex++]; + if (!strcmp(enhancement, "EED")) + diffusionFilter->SetEnhancement(DiffusionFilterType::EED); // Weickert's exponent : 4. + else if (!strcmp(enhancement, "cEED")) + diffusionFilter->SetEnhancement(DiffusionFilterType::cEED); // Weickert's exponent : 4. + else if (!strcmp(enhancement, "CED")) + diffusionFilter->SetEnhancement(DiffusionFilterType::CED); // Weickert's exponent : 2. + else if (!strcmp(enhancement, "cCED")) + diffusionFilter->SetEnhancement(DiffusionFilterType::cCED); // Weickert's exponent : 2. + else if (!strcmp(enhancement, "Isotropic")) + diffusionFilter->SetEnhancement(DiffusionFilterType::Isotropic); // Perona-Mali's exponent: 2. + else + itkGenericExceptionMacro("Error: Unrecognized enhancement (fifth argument).\n"); + } + + if (argIndex < argc) + { + const double noiseScale = std::stod(argv[argIndex++]); + if (noiseScale == 0.) + itkGenericExceptionMacro("Error: Unrecognized noiseScale (sixth argument).\n"); + diffusionFilter->SetNoiseScale(noiseScale); + } + + if (argIndex < argc) + { + const double featureScale = std::stod(argv[argIndex++]); + if (featureScale == 0.) + itkGenericExceptionMacro("Error: Unrecognized featureScale (seventh argument).\n"); + diffusionFilter->SetFeatureScale(featureScale); + } + + if (argIndex < argc) + { + const double exponent = std::stod(argv[argIndex++]); + if (exponent == 0.) + itkGenericExceptionMacro("Error: Unrecognized exponent (eighth argument).\n"); + diffusionFilter->SetExponent(exponent); + } + + if (argIndex < argc) + { + itkGenericExceptionMacro("Error: excessive number of arguments"); + } + + /*{ + std::cerr << + "T: " << diffusionFilter->GetDiffusionTime() << "\n" << + "Lambda: " << diffusionFilter->GetLambda() << "\n" << + "argc: " << argc << "\n"; + + diffusionFilter->Update(); + auto image = diffusionFilter->GetOutput(); + std::cerr << + image->GetBufferedRegion() << "\n\n" << + "pixel export size: " << sizeof(ExportPixelType) << "\n"; +; + + }*/ + + using ExportImageType = Image; + using CasterType = CastImageFilter; + typename CasterType::Pointer caster = CasterType::New(); + caster->SetInput(diffusionFilter->GetOutput()); + + // using ScalarImageType = typename DiffusionFilterType::ScalarImageType; + using WriterType = ImageFileWriter; + typename WriterType::Pointer writer = WriterType::New(); + writer->SetInput(caster->GetOutput()); + writer->SetFileName(outputFileName); + + itk::TimeProbe clock; + clock.Start(); + writer->Update(); + clock.Stop(); + std::cout << "Took: " << clock.GetMean() << " seconds\n"; + + return EXIT_SUCCESS; +} + +} // end namespace CoherenceEnhancingDiffusionCommandLine + +#endif diff --git a/Modules/Filtering/AnisotropicDiffusionLBR/include/LinearAnisotropicDiffusionCommandLine.h b/Modules/Filtering/AnisotropicDiffusionLBR/include/LinearAnisotropicDiffusionCommandLine.h new file mode 100644 index 00000000000..5cfed7d8b6a --- /dev/null +++ b/Modules/Filtering/AnisotropicDiffusionLBR/include/LinearAnisotropicDiffusionCommandLine.h @@ -0,0 +1,252 @@ +/*========================================================================= + * + * 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. + * + *=========================================================================*/ + +// +// Created by Jean-Marie Mirebeau on 01/12/2014. +// +// + +#ifndef itkLinearAnisotropicDiffusionCommandLine_h +#define itkLinearAnisotropicDiffusionCommandLine_h + +#include "itkImageFileReader.h" +#include "itkImageFileWriter.h" +#include "itkCoherenceEnhancingDiffusionImageFilter.h" + +namespace LinearAnisotropicDiffusionCommandLine +{ + +inline void +Usage() +{ + std::cerr << "Input image filename. 2D and 3D images supported. Required.\n" + << "Output tensor field filename. Required.\n" + << "Diffusion time. Required.\n" + << "Output image filename. Required.\n" + << "RatioToMaxStableStep. Range: ]0,1]. Default: 0.7. Optionnal.\n" + << "MaxNumberOfIterations. Range: 1...Infinity. Default: 200. Optionnal.\n" + << "\n"; +} + +using namespace itk; + +int +Execute(int argc, char * argv[]); + +template +int +Execute(int argc, char * argv[], itk::ImageIOBase::IOComponentEnum, int nComponents); + +template +int +Execute(int argc, char * argv[], int nComponents); + +template +int +Execute(int argc, char * argv[]); + + +struct ReportProgressToCOutType : public itk::Command +{ + itkNewMacro(ReportProgressToCOutType); + void + Execute(itk::Object * caller, const itk::EventObject & event) override + { + Execute((const itk::Object *)caller, event); + } + + void + Execute(const itk::Object * object, const itk::EventObject &) override + { + std::cout << object->GetNameOfClass() + << " has completed: " << int(100 * dynamic_cast(object)->GetProgress()) << "%" + << std::endl; + } +}; + + +inline int +Execute(int argc, char * argv[]) +{ + using std::cerr; + using std::endl; + using namespace itk; + + if (argc < 4 + 1) + { + Usage(); + return EXIT_SUCCESS; + } + + const char * imageFileName = argv[0 + 1]; + using ReaderType = ImageFileReader>; + ReaderType::Pointer reader = ReaderType::New(); + reader->SetFileName(imageFileName); + + reader->UpdateOutputInformation(); + + const ImageIOBase * io = reader->GetImageIO(); + const int ImageDimension = io->GetNumberOfDimensions(); + const itk::ImageIOBase::IOComponentEnum componentType = io->GetComponentType(); + const int nComponents = io->GetNumberOfComponents(); + + { + const char * tensorFileName = argv[1 + 1]; + ReaderType::Pointer reader2 = ReaderType::New(); + reader2->SetFileName(tensorFileName); + reader2->UpdateOutputInformation(); + const ImageIOBase * io2 = reader2->GetImageIO(); + if (io2->GetComponentType() != componentType) + std::cerr << "Warning: image and tensors have distinct component types.\n"; + if (ImageDimension != (int)io2->GetNumberOfDimensions()) + itkGenericExceptionMacro("Error: image and tensor image have distinct dimension"); + const int TensorDimension = (ImageDimension * (ImageDimension + 1)) / 2; + if (TensorDimension != (int)io2->GetNumberOfComponents()) + itkGenericExceptionMacro("Error: wrong tensor dimension, should be n*(n+1)/2 where n=ImageDimension."); + if (io2->GetPixelType() == itk::IOPixelEnum::SYMMETRICSECONDRANKTENSOR) + std::cerr << "Warning: tensor image pixel type not marked as Symmetric Second Rank Tensor.\n"; + } + + switch (ImageDimension) + { + case 2: + return Execute<2>(argc, argv, componentType, nComponents); + // case 3: return Execute<3>(argc,argv,componentType,nComponents); + default: + itkGenericExceptionMacro("Sorry, unsupported image dimension."); + } +} + +template +int +Execute(int argc, char * argv[], itk::ImageIOBase::IOComponentEnum componentType, int nComponents) +{ + switch (componentType) + { + case itk::IOComponentEnum::UCHAR: + return Execute(argc, argv, nComponents); + case itk::IOComponentEnum::FLOAT: + return Execute(argc, argv, nComponents); + case itk::IOComponentEnum::DOUBLE: + return Execute(argc, argv, nComponents); + default: + itkGenericExceptionMacro("Sorry, unsupported component type"); + } +} + +template +int +Execute(int argc, char * argv[], int nComponents) +{ + switch (nComponents) + { + case 1: + return Execute(argc, argv); + case 2: + return Execute, Vector>(argc, argv); + case 3: + return Execute, Vector>(argc, argv); + default: + itkGenericExceptionMacro("Sorry, unsupported number of components"); + } +} + +template +int +Execute(int argc, char * argv[]) +{ + + // Import image + using ImageType = Image; + using ReaderType = ImageFileReader; + typename ReaderType::Pointer reader = ReaderType::New(); + const char * imageFileName = argv[0 + 1]; + reader->SetFileName(imageFileName); + + // Import tensors + using TensorType = SymmetricSecondRankTensor; + using TensorImageType = Image; + using TensorReaderType = ImageFileReader; + typename TensorReaderType::Pointer tensorReader = TensorReaderType::New(); + const char * tensorImageFileName = argv[1 + 1]; + tensorReader->SetFileName(tensorImageFileName); + + // Import diffusion time + const double diffusionTime = std::stod(argv[2 + 1]); + if (diffusionTime == 0) + itkGenericExceptionMacro("Error: Unrecognized diffusion time (third argument).\n"); + + // Import output image filename + const char * outputFileName = argv[3 + 1]; + + // Setup diffusion filter + using DiffusionFilterType = LinearAnisotropicDiffusionLBRImageFilter; + typename DiffusionFilterType::Pointer diffusionFilter = DiffusionFilterType::New(); + diffusionFilter->SetInputImage(reader->GetOutput()); + diffusionFilter->SetInputTensor(tensorReader->GetOutput()); + diffusionFilter->SetMaxDiffusionTime(diffusionTime); + + int argIndex = 4 + 1; + if (argIndex < argc) + { + const double ratioToMaxStableTimeStep = std::stod(argv[argIndex++]); + if (ratioToMaxStableTimeStep == 0) + itkGenericExceptionMacro("Error: Unrecognized RatioToMaxStableTimeStep (fourth argument).\n"); + diffusionFilter->SetRatioToMaxStableTimeStep(ratioToMaxStableTimeStep); + } + + if (argIndex < argc) + { + const int maxNumberOfTimeSteps = std::stoi(argv[argIndex++]); + if (maxNumberOfTimeSteps == 0) + itkGenericExceptionMacro("Error: Unrecognized maxNumberOfTimeSteps (fifth argument).\n"); + diffusionFilter->SetMaxNumberOfTimeSteps(maxNumberOfTimeSteps); + } + else + diffusionFilter->SetMaxNumberOfTimeSteps(200); + + ReportProgressToCOutType::Pointer reportDiffusionProgress = ReportProgressToCOutType::New(); + diffusionFilter->AddObserver(ProgressEvent(), reportDiffusionProgress); + + using ExportImageType = Image; + using CasterType = CastImageFilter; + typename CasterType::Pointer caster = CasterType::New(); + caster->SetInput(diffusionFilter->GetOutput()); + + // using ScalarImageType = typename DiffusionFilterType::ScalarImageType; + using WriterType = ImageFileWriter; + typename WriterType::Pointer writer = WriterType::New(); + writer->SetInput(caster->GetOutput()); + writer->SetFileName(outputFileName); + writer->Update(); + + const ScalarType effectiveDiffusionTime = diffusionFilter->GetEffectiveDiffusionTime(); + if (effectiveDiffusionTime < 0.99 * diffusionTime) + { + std::cerr << "Warning: early abort at effective diffusion time: " << effectiveDiffusionTime + << ", you may want to increase the max number of time steps: " + << diffusionFilter->GetMaxNumberOfTimeSteps() << "\n"; + Usage(); + } + + return EXIT_SUCCESS; +} + +} // end namespace LinearAnisotropicDiffusionCommandLine + +#endif diff --git a/Modules/Filtering/AnisotropicDiffusionLBR/include/itkAnisotropicDiffusionLBRImageFilter.h b/Modules/Filtering/AnisotropicDiffusionLBR/include/itkAnisotropicDiffusionLBRImageFilter.h new file mode 100644 index 00000000000..8c9030f1de6 --- /dev/null +++ b/Modules/Filtering/AnisotropicDiffusionLBR/include/itkAnisotropicDiffusionLBRImageFilter.h @@ -0,0 +1,146 @@ +/*========================================================================= + * + * 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. + * + *=========================================================================*/ + + +// +// Created by Jean-Marie Mirebeau on 28/02/2014. +// +// + +#ifndef itkAnisotropicDiffusionLBRImageFilter_h +#define itkAnisotropicDiffusionLBRImageFilter_h + +#include "itkLinearAnisotropicDiffusionLBRImageFilter.h" +#include "itkStructureTensorImageFilter.h" + +namespace itk +{ + +/** \class AnisotropicDiffusionLBRImageFilter + * + * \brief Non-linear anisotropic diffusion using lattice basis reduction. + * + * This class repeatedly calls the LinearAnisotropicDiffusionLBRImageFilter, + * with non-linear diffusion tensors built on the fly. These tensors are + * obtained by computing the image structure tensors, and appropriately + * modifying their eigenvalues with the method EigenValuesTransform. The + * latter method is not implemented, and needs to be provided in a subclass, + * such as CoherenceEnhancingDiffusionImageFilter. + * + * \ingroup AnisotropicDiffusionLBR + */ +template ::RealType> +class AnisotropicDiffusionLBRImageFilter : public ImageToImageFilter +{ +public: + ITK_DISALLOW_COPY_AND_MOVE(AnisotropicDiffusionLBRImageFilter); + + using Self = AnisotropicDiffusionLBRImageFilter; + using Superclass = ImageToImageFilter; + using Pointer = SmartPointer; + using ConstPointer = SmartPointer; + + /// Method for creation through the object factory. + itkNewMacro(Self); + /// Run-time type information (and related methods). + itkOverrideGetNameOfClassMacro(AnisotropicDiffusionLBRImageFilter); + + using ImageType = TImage; + using PixelType = typename ImageType::PixelType; + using ScalarType = TScalar; + + using ImageDimensionType = typename ImageType::ImageDimensionType; + static constexpr ImageDimensionType ImageDimension = ImageType::ImageDimension; + + using TensorType = SymmetricSecondRankTensor; + using TensorImageType = Image; + + using StructureTensorFilterType = StructureTensorImageFilter; + using LinearDiffusionFilterType = LinearAnisotropicDiffusionLBRImageFilter; + + /** Passed to a StructureTensorImageFilter. */ + itkSetMacro(NoiseScale, ScalarType); + itkGetConstMacro(NoiseScale, ScalarType); + itkSetMacro(FeatureScale, ScalarType); + itkGetConstMacro(FeatureScale, ScalarType); + + /** Passed to a LinearAnisotropicDiffusion Filter. */ + itkSetMacro(RatioToMaxStableTimeStep, ScalarType); + itkGetConstMacro(RatioToMaxStableTimeStep, ScalarType); + itkSetMacro(MaxTimeStepsBetweenTensorUpdates, int); + itkGetConstMacro(MaxTimeStepsBetweenTensorUpdates, int); + + itkSetMacro(DiffusionTime, ScalarType); + itkGetConstMacro(DiffusionTime, ScalarType); + + /** If true, uses unit pixel spacing, and rescales structure + * tensors for uni maximum trace. */ + itkSetMacro(Adimensionize, bool); + itkGetConstMacro(Adimensionize, bool); + + using EigenValuesArrayType = typename TensorType::EigenValuesArrayType; + /** Transformation of the Structure tensor eigenvalues into the diffusion + * tensor eigenvalues. Needs to be overloaded in a subclass. + * (Structure tensor eigenvalues are sorted by increasing order for convenience). */ + virtual EigenValuesArrayType + EigenValuesTransform(const EigenValuesArrayType &) const + { + itkExceptionMacro("Undefined tensor eigenvalues transform"); + } + + virtual typename TensorImageType::Pointer + GetLastTensorImage() + { + return m_TensorImage; + } + using EffectiveTimesAndIterationsType = std::vector>; + itkGetConstReferenceMacro(LinearFilterEffectiveTimesAndIterations, EffectiveTimesAndIterationsType); + +protected: + ScalarType m_NoiseScale; + ScalarType m_FeatureScale; + + ScalarType m_RatioToMaxStableTimeStep; + int m_MaxTimeStepsBetweenTensorUpdates{ 5 }; + + AnisotropicDiffusionLBRImageFilter(); + ~AnisotropicDiffusionLBRImageFilter() override = default; + + typename TensorImageType::Pointer m_TensorImage; + + virtual void + ComputeDiffusionTensors(ImageType *); + + ScalarType m_DiffusionTime; + bool m_Adimensionize{ true }; + + void + GenerateData() override; + + EffectiveTimesAndIterationsType m_LinearFilterEffectiveTimesAndIterations; + + struct DiffusionTensorFunctor; +}; + +} // end namespace itk + +#ifndef ITK_MANUAL_INSTANTIATION +# include "itkAnisotropicDiffusionLBRImageFilter.hxx" +#endif + +#endif diff --git a/Modules/Filtering/AnisotropicDiffusionLBR/include/itkAnisotropicDiffusionLBRImageFilter.hxx b/Modules/Filtering/AnisotropicDiffusionLBR/include/itkAnisotropicDiffusionLBRImageFilter.hxx new file mode 100644 index 00000000000..dfa544ff390 --- /dev/null +++ b/Modules/Filtering/AnisotropicDiffusionLBR/include/itkAnisotropicDiffusionLBRImageFilter.hxx @@ -0,0 +1,169 @@ +/*========================================================================= + * + * 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. + * + *=========================================================================*/ + +// +// Created by Jean-Marie Mirebeau on 28/02/2014. +// +// + +#ifndef itkAnisotropicDiffusionLBRImageFilter_hxx +#define itkAnisotropicDiffusionLBRImageFilter_hxx + + +namespace itk +{ + +template +AnisotropicDiffusionLBRImageFilter::AnisotropicDiffusionLBRImageFilter() + : m_NoiseScale(0.5) + , m_FeatureScale(2) + , m_RatioToMaxStableTimeStep(0.7) + , + + m_DiffusionTime(1) + +{} + + +template +void +AnisotropicDiffusionLBRImageFilter::GenerateData() +{ + typename ImageType::Pointer inputImage = const_cast(this->GetInput()); + typename ImageType::Pointer image = inputImage; + + using SpacingType = typename ImageType::SpacingType; + const SpacingType referenceSpacing = inputImage->GetSpacing(); + + // const SpacingType unitSpacing(1); // Better below for non-uniform spacing. + double minSpacing = referenceSpacing[0]; + for (ImageDimensionType i = 1; i < ImageDimension; ++i) + { + minSpacing = std::min(minSpacing, referenceSpacing[i]); + } + const SpacingType unitSpacing = referenceSpacing / minSpacing; + + if (m_Adimensionize) + { + inputImage->SetSpacing(unitSpacing); + } + + ScalarType remainingTime = m_DiffusionTime; + + while (remainingTime > 0) + { + ComputeDiffusionTensors(image); + typename LinearDiffusionFilterType::Pointer linearDiffusionFilter = LinearDiffusionFilterType::New(); + linearDiffusionFilter->SetMaxNumberOfTimeSteps(m_MaxTimeStepsBetweenTensorUpdates); + linearDiffusionFilter->SetRatioToMaxStableTimeStep(m_RatioToMaxStableTimeStep); + + linearDiffusionFilter->SetInputImage(image); + linearDiffusionFilter->SetInputTensor(m_TensorImage); + linearDiffusionFilter->SetMaxDiffusionTime(remainingTime); + linearDiffusionFilter->Update(); + image = linearDiffusionFilter->GetOutput(); + remainingTime -= linearDiffusionFilter->GetEffectiveDiffusionTime(); + + m_LinearFilterEffectiveTimesAndIterations.push_back(std::pair( + linearDiffusionFilter->GetEffectiveDiffusionTime(), linearDiffusionFilter->GetEffectiveNumberOfTimeSteps())); + + this->UpdateProgress(1. - remainingTime / m_DiffusionTime); + } + + if (m_Adimensionize) + { + inputImage->SetSpacing(referenceSpacing); + image->SetSpacing(referenceSpacing); + } + this->GraftOutput(image); +} + + +template +struct AnisotropicDiffusionLBRImageFilter::DiffusionTensorFunctor +{ + Self * eigenValuesFunctor; + struct OrderingType; + TensorType + operator()(const TensorType & S) + { + EigenValuesArrayType eigenValues; + typename TensorType::EigenVectorsMatrixType eigenVectors; + S.ComputeEigenAnalysis(eigenValues, eigenVectors); + + // For convenience, eigenvalues are sorted by increasing order + Vector order; + for (ImageDimensionType i = 0; i < ImageDimension; ++i) + order[i] = i; + + OrderingType ordering(eigenValues); + + std::sort(order.Begin(), order.End(), ordering); + + std::sort(eigenValues.Begin(), eigenValues.End()); + EigenValuesArrayType ev = this->eigenValuesFunctor->EigenValuesTransform(eigenValues); + + TensorType DiffusionTensor; + for (ImageDimensionType i = 0; i < ImageDimension; ++i) + { + DiffusionTensor(order[i], order[i]) = ev[i]; + for (ImageDimensionType j = 0; j < i; ++j) + DiffusionTensor(i, j) = 0.; + } + return DiffusionTensor.Rotate(eigenVectors.GetTranspose()); + } +}; + + +// c++ 11 would be : [& eigenValues](int i, int j)->bool {return eigenValues[i] +struct AnisotropicDiffusionLBRImageFilter::DiffusionTensorFunctor ::OrderingType +{ + bool + operator()(int i, int j) const + { + return this->e[i] < this->e[j]; + } + const EigenValuesArrayType & e; + OrderingType(const EigenValuesArrayType & e_) + : e(e_) {}; +}; + +template +void +AnisotropicDiffusionLBRImageFilter::ComputeDiffusionTensors(ImageType * image) +{ + typename StructureTensorFilterType::Pointer structureTensorFilter = StructureTensorFilterType::New(); + + structureTensorFilter->SetNoiseScale(m_NoiseScale); + structureTensorFilter->SetFeatureScale(m_FeatureScale); + structureTensorFilter->SetRescaleForUnitMaximumTrace(m_Adimensionize); + structureTensorFilter->SetInput(image); + + using ImageFunctorType = UnaryFunctorImageFilter; + typename ImageFunctorType::Pointer imageFunctor = ImageFunctorType::New(); + imageFunctor->GetFunctor().eigenValuesFunctor = this; + imageFunctor->SetInput(structureTensorFilter->GetOutput()); + + imageFunctor->Update(); + m_TensorImage = imageFunctor->GetOutput(); +} + +} // end namespace itk + +#endif diff --git a/Modules/Filtering/AnisotropicDiffusionLBR/include/itkAnisotropicDiffusionLBRMacro.h b/Modules/Filtering/AnisotropicDiffusionLBR/include/itkAnisotropicDiffusionLBRMacro.h new file mode 100644 index 00000000000..b9cf06c274a --- /dev/null +++ b/Modules/Filtering/AnisotropicDiffusionLBR/include/itkAnisotropicDiffusionLBRMacro.h @@ -0,0 +1,40 @@ +/*========================================================================= + * + * 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. + * + *=========================================================================*/ + +// +// Created by Jean-Marie Mirebeau on 07/03/2014. +// +// + +#ifndef itkAnisotropicDiffusionLBRMacro_h +#define itkAnisotropicDiffusionLBRMacro_h + +/** + Getters and setters for functor types, inspired by UnaryFunctorImageFilter. + No equality test performed for the setter. + */ +#define GetSetFunctorMacro(name, type) \ + virtual type & Get##name() { return this->m_##name; } \ + virtual const type & Get##name() const { return this->m_##name; } \ + virtual void Set##name(const type & _arg) \ + { \ + m_##name = _arg; \ + this->Modified(); \ + } + +#endif diff --git a/Modules/Filtering/AnisotropicDiffusionLBR/include/itkCoherenceEnhancingDiffusionImageFilter.h b/Modules/Filtering/AnisotropicDiffusionLBR/include/itkCoherenceEnhancingDiffusionImageFilter.h new file mode 100644 index 00000000000..9c3a498ff6b --- /dev/null +++ b/Modules/Filtering/AnisotropicDiffusionLBR/include/itkCoherenceEnhancingDiffusionImageFilter.h @@ -0,0 +1,137 @@ +/*========================================================================= + * + * 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. + * + *=========================================================================*/ + +// +// Created by Jean-Marie Mirebeau on 06/03/2014. +// +// + +#ifndef itkCoherenceEnhancingDiffusionImageFilter_h +#define itkCoherenceEnhancingDiffusionImageFilter_h + +#include "itkAnisotropicDiffusionLBRImageFilter.h" + + +namespace itk +{ +/** + * \class CoherenceEnhancingDiffusionImageFilter + * + * \brief Coherence enhancing diffusion and edge enhancing diffusion. + * + * Implementation of Coherence Enhancing Diffusion (CED), and + * Edge Enhancing Diffusion (EED), as described by Weickert. + * + * CED heuristically smoothes everywhere except accross image contours, + * while EED smoothes nowhere but tangentially to image contours. + * + * The non-linear diffusion tensor is defined in terms of the structure tensor. + * + * Denote by \f$\mu_i\f$ the structure tensor eigenvalues, at a given point \f$x\f$, + * with \f$0\leq i < d\f$. Let also \f$\mu_{\rm min}\f$ and \f$\mu_{\rm max}\f$, + * be the smallest and largest eigenvalues respectively. The diffusion tensor is + * defined by the same eigenvectors, but with modified with eigenvalues \f$\lambda_i\f$. + * + * Edge Enhancing Diffusion: + * + * \f$\lambda_i := g(\mu_i - \mu_{\rm min})\f$, where \f$g(s) = 1 - (1-\alpha)*exp(-(\lambda/s)^m)\f$ + * + * Note the limit values \f$g(0) = 1\f$, \f$g(\infty) = \alpha\f$. + * + * Coherence Enhancing Diffusion: + * + * \f$\lambda_i := g(\mu_{\rm max} - \mu_i)\f$, where \f$g(s) = \alpha + (1-\alpha)*exp(-(\lambda/s)^m)\f$ + * + * Note the limit values \f$g(0) = \alpha\f$, \f$g(\infty) = 1\f$. + * + * \ingroup AnisotropicDiffusionLBR + */ +template ::RealType> +class CoherenceEnhancingDiffusionImageFilter : public AnisotropicDiffusionLBRImageFilter +{ +public: + ITK_DISALLOW_COPY_AND_MOVE(CoherenceEnhancingDiffusionImageFilter); + + using Self = CoherenceEnhancingDiffusionImageFilter; + using Superclass = AnisotropicDiffusionLBRImageFilter; + using Pointer = SmartPointer; + using ConstPointer = SmartPointer; + + /** Method for creation through the object factory. */ + itkNewMacro(Self); + + /** Run-time type information (and related methods). */ + itkOverrideGetNameOfClassMacro(CoherenceEnhancingDiffusionImageFilter); + + using InputImageDimensionType = typename Superclass::InputImageType::ImageDimensionType; + static constexpr InputImageDimensionType InputImageDimension = Superclass::InputImageType::ImageDimension; + + using EigenValuesArrayType = typename Superclass::EigenValuesArrayType; + EigenValuesArrayType + EigenValuesTransform(const EigenValuesArrayType &) const override; + + using ScalarType = typename Superclass::ScalarType; + /** Exponent m involved in the function g defining eigenvalues. */ + itkSetMacro(Exponent, ScalarType); + itkSetMacro(Lambda, ScalarType); + itkSetMacro(Alpha, ScalarType); + + itkGetMacro(Exponent, ScalarType); + itkGetMacro(Lambda, ScalarType); + itkGetMacro(Alpha, ScalarType); + + enum EnhancementType + { + CED, + cCED, + EED, + cEED, + Isotropic + }; + /// Switch between CED, EED, and variants. + itkSetEnumMacro(Enhancement, EnhancementType); + itkGetEnumMacro(Enhancement, EnhancementType); + +protected: + ScalarType m_Lambda; + ScalarType m_Exponent; + ScalarType m_Alpha; + EnhancementType m_Enhancement; + + ScalarType + g_CED(ScalarType s) const + { + return s <= 0 ? m_Alpha : m_Alpha + (1 - m_Alpha) * exp(-pow(m_Lambda / s, m_Exponent)); + } + ScalarType + g_EED(ScalarType s) const + { + return s <= 0 ? 1 : 1 - (1 - m_Alpha) * exp(-pow(m_Lambda / s, m_Exponent)); + } + + CoherenceEnhancingDiffusionImageFilter(); + ~CoherenceEnhancingDiffusionImageFilter() override = default; +}; + +} // end namespace itk + +#ifndef ITK_MANUAL_INSTANTIATION +# include "itkCoherenceEnhancingDiffusionImageFilter.hxx" +#endif + +#endif diff --git a/Modules/Filtering/AnisotropicDiffusionLBR/include/itkCoherenceEnhancingDiffusionImageFilter.hxx b/Modules/Filtering/AnisotropicDiffusionLBR/include/itkCoherenceEnhancingDiffusionImageFilter.hxx new file mode 100644 index 00000000000..1e00bb70bd1 --- /dev/null +++ b/Modules/Filtering/AnisotropicDiffusionLBR/include/itkCoherenceEnhancingDiffusionImageFilter.hxx @@ -0,0 +1,93 @@ +/*========================================================================= + * + * 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 itkCoherenceEnhancingDiffusionImageFilter_hxx +#define itkCoherenceEnhancingDiffusionImageFilter_hxx + + +namespace itk +{ + +template +CoherenceEnhancingDiffusionImageFilter::CoherenceEnhancingDiffusionImageFilter() + : m_Lambda(0.05) + , m_Exponent(2) + , m_Alpha(0.01) + , m_Enhancement(CED) +{} + + +template +typename CoherenceEnhancingDiffusionImageFilter::EigenValuesArrayType +CoherenceEnhancingDiffusionImageFilter::EigenValuesTransform(const EigenValuesArrayType & ev0) const +{ + const ScalarType evMin = ev0[0]; + const ScalarType evMax = ev0[InputImageDimension - 1]; + + EigenValuesArrayType ev; + switch (m_Enhancement) + { + // Weickert's filter. + case CED: + for (InputImageDimensionType i = 0; i < InputImageDimension; ++i) + { + ev[i] = g_CED(evMax - ev0[i]); + } + break; + + // A variance, requiring stronger coherence. + case cCED: + for (InputImageDimensionType i = 0; i < InputImageDimension; ++i) + { + ev[i] = g_CED((evMax - ev0[i]) / (1. + ev0[i] / m_Lambda)); + } + break; + + // Weickert's filter. + case EED: + for (InputImageDimensionType i = 0; i < InputImageDimension; ++i) + { + ev[i] = g_EED(ev0[i] - evMin); + } + break; + + // A variant, promoting diffusion in at least one direction at each point. + case cEED: + for (InputImageDimensionType i = 0; i < InputImageDimension; ++i) + { + ev[i] = g_EED(ev0[i]); + } + break; + + // Isotropic tensors, closely related to Perona-Malik's approach. + case Isotropic: + for (InputImageDimensionType i = 0; i < InputImageDimension; ++i) + { + ev[i] = g_EED(evMax); + } + break; + + default: + itkExceptionMacro("Unsupported diffusion type"); + } + return ev; +} + +} // end namespace itk + +#endif diff --git a/Modules/Filtering/AnisotropicDiffusionLBR/include/itkLinearAnisotropicDiffusionLBRImageFilter.h b/Modules/Filtering/AnisotropicDiffusionLBR/include/itkLinearAnisotropicDiffusionLBRImageFilter.h new file mode 100644 index 00000000000..92f3c647afe --- /dev/null +++ b/Modules/Filtering/AnisotropicDiffusionLBR/include/itkLinearAnisotropicDiffusionLBRImageFilter.h @@ -0,0 +1,175 @@ +/*========================================================================= + * + * 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. + * + *=========================================================================*/ + +// +// Created by Jean-Marie Mirebeau on 28/02/2014. +// +// + +#ifndef itkLinearAnisotropicDiffusionLBRImageFilter_h +#define itkLinearAnisotropicDiffusionLBRImageFilter_h + + +#include "itkImageToImageFilter.h" +#include "itkSymmetricSecondRankTensor.h" + +namespace itk +{ +/** + * \class LinearAnisotropicDiffusionLBRImageFilter + * + * \brief Anisotropic diffusion using lattice basis reduction. + * + * \f[\partial_t u = {\rm div} (D \nabla u),\f] + * + * with Neumann boundary conditions. The numerical scheme is stable and + * satisfies the maximum principle, even for strongly anisotropic tensors, + * thanks to an adaptive discretization using arithmetic techniques + * (Lattice Basis Reduction, LBR). + * + * \ingroup AnisotropicDiffusionLBR + */ +template ::RealType> +class LinearAnisotropicDiffusionLBRImageFilter : public ImageToImageFilter +{ +public: + ITK_DISALLOW_COPY_AND_MOVE(LinearAnisotropicDiffusionLBRImageFilter); + + /** Standard class type alias. */ + using Self = LinearAnisotropicDiffusionLBRImageFilter; + using Superclass = ImageToImageFilter; + using Pointer = SmartPointer; + using ConstPointer = SmartPointer; + + /** Method for creation through the object factory. */ + itkNewMacro(Self); + + /** Run-time type information (and related methods). */ + itkOverrideGetNameOfClassMacro(LinearAnisotropicDiffusionLBRImageFilter); + + using ImageType = TImage; + using PixelType = typename ImageType::PixelType; + + using ImageDimensionType = typename ImageType::ImageDimensionType; + static constexpr ImageDimensionType ImageDimension = ImageType::ImageDimension; + + using ScalarType = TScalar; + using TensorType = SymmetricSecondRankTensor; + using TensorImageType = Image; + using RegionType = ImageRegion; + + void + SetInputImage(const ImageType * image); + void + SetInputTensor(const TensorImageType * tensorImage); + + void + SetMaxDiffusionTime(ScalarType time); + itkGetConstMacro(DiffusionTime, ScalarType); + + void + SetMaxNumberOfTimeSteps(int n); + itkGetConstMacro(MaxNumberOfTimeSteps, int); + + void + SetRatioToMaxStableTimeStep(ScalarType ratio); + itkGetConstMacro(RatioToMaxStableTimeStep, ScalarType); + + itkGetConstMacro(EffectiveDiffusionTime, ScalarType); + itkGetConstMacro(EffectiveNumberOfTimeSteps, int); + +protected: + LinearAnisotropicDiffusionLBRImageFilter(); + ~LinearAnisotropicDiffusionLBRImageFilter() override = default; + + typename ImageType::ConstPointer + GetInputImage(); + typename TensorImageType::ConstPointer + GetInputTensor(); + + using IndexType = Index; + + // ******* Containers for the stencils used in the discretization + static const unsigned int HalfStencilSize = (ImageDimension == 2) ? 3 : 6; + static const unsigned int StencilSize = 2 * HalfStencilSize; + + using StencilCoefficientsType = Vector; + using OffsetType = Offset; + using StencilOffsetsType = Vector; + + using InternalSizeT = int; + using StencilBufferIndicesType = Vector; + + + // *************** Computation ***************** + void + GenerateData() override; + virtual void + GenerateStencils(); /// Automatically called by GenerateData + virtual void + ImageUpdateLoop(); /// Automatically called by GenerateData + + using StencilType = std::pair; + using StencilImageType = Image; + typename StencilImageType::Pointer m_StencilImage; + + using ScalarImageType = Image; + typename ScalarImageType::Pointer m_DiagonalCoefficients; + + virtual ScalarType + MaxStableTimeStep(); + + ScalarType m_DiffusionTime; + ScalarType m_RatioToMaxStableTimeStep; + int m_MaxNumberOfTimeSteps{ 10 }; + + ScalarType m_EffectiveDiffusionTime; + int m_EffectiveNumberOfTimeSteps{ 0 }; + + virtual void + ImageUpdate(ScalarType delta); + typename ImageType::Pointer m_PreviousImage; + typename ImageType::Pointer m_NextImage; + + virtual RegionType + GetRequestedRegion() + { + return GetInputImage()->GetRequestedRegion(); + } + + InternalSizeT + OutsideBufferIndex() const + { + return NumericTraits::max(); + } + + struct StencilFunctor; + struct FunctorType; + + using VectorType = Vector; + static ScalarType + ScalarProduct(const TensorType &, const VectorType &, const VectorType &); +}; +} // end namespace itk + + +#ifndef ITK_MANUAL_INSTANTIATION +# include "itkLinearAnisotropicDiffusionLBRImageFilter.hxx" +#endif + +#endif diff --git a/Modules/Filtering/AnisotropicDiffusionLBR/include/itkLinearAnisotropicDiffusionLBRImageFilter.hxx b/Modules/Filtering/AnisotropicDiffusionLBR/include/itkLinearAnisotropicDiffusionLBRImageFilter.hxx new file mode 100644 index 00000000000..e6027776aef --- /dev/null +++ b/Modules/Filtering/AnisotropicDiffusionLBR/include/itkLinearAnisotropicDiffusionLBRImageFilter.hxx @@ -0,0 +1,544 @@ +/*========================================================================= + * + * 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. + * + *=========================================================================*/ + +// +// Created by Jean-Marie Mirebeau on 28/02/2014. +// +// + +#ifndef itkLinearAnisotropicDiffusionLBRImageFilter_hxx +#define itkLinearAnisotropicDiffusionLBRImageFilter_hxx + +#include "itkUnaryFunctorImageFilter.h" +#include "itkImageRegionIterator.h" +#include "itkMinimumMaximumImageCalculator.h" +#include "itkCastImageFilter.h" +#include "itkExtractImageFilter.h" +#include "itkUnaryFunctorWithIndexImageFilter.h" +#include "itkTernaryFunctorImageFilter.h" + +namespace itk +{ + +template +LinearAnisotropicDiffusionLBRImageFilter::LinearAnisotropicDiffusionLBRImageFilter() + : m_DiffusionTime(1) + , m_RatioToMaxStableTimeStep(0.7) + , + + m_EffectiveDiffusionTime(0) + +{ + this->SetNumberOfRequiredInputs(2); +} + + +template +void +LinearAnisotropicDiffusionLBRImageFilter::SetInputImage(const ImageType * image) +{ + this->SetNthInput(0, const_cast(image)); +} + + +template +void +LinearAnisotropicDiffusionLBRImageFilter::SetInputTensor(const TensorImageType * tensorImage) +{ + this->SetNthInput(1, const_cast(tensorImage)); +} + + +template +typename TImage::ConstPointer +LinearAnisotropicDiffusionLBRImageFilter::GetInputImage() +{ + return static_cast(this->ProcessObject::GetInput(0)); +} + + +template +typename LinearAnisotropicDiffusionLBRImageFilter::TensorImageType::ConstPointer +LinearAnisotropicDiffusionLBRImageFilter::GetInputTensor() +{ + return static_cast(this->ProcessObject::GetInput(1)); +} + + +template +void +LinearAnisotropicDiffusionLBRImageFilter::GenerateData() +{ + GenerateStencils(); + this->UpdateProgress(0.5); + + this->ImageUpdateLoop(); +} + +// **************************** Computation *********************** +template +struct LinearAnisotropicDiffusionLBRImageFilter::StencilFunctor +{ +public: + using SpacingType = typename TensorImageType::SpacingType; + void + Initialize(RegionType region_, SpacingType spacing) + { + region = region_; + prod[0] = 1; + for (ImageDimensionType i = 1; i < ImageDimension; ++i) + { + prod[i] = prod[i - 1] * region.GetSize()[i - 1]; + } + for (ImageDimensionType i = 0; i < ImageDimension; ++i) + { + invSpacing[i] = ScalarType(1) / spacing[i]; + } + } + + InternalSizeT + BufferIndex(const IndexType & x) const + { + IndexValueType ans = 0; + for (ImageDimensionType i = 0; i < ImageDimension; ++i) + { + ans += this->prod[i] * (x[i] - this->region.GetIndex()[i]); + } + return ans; + } + + StencilType + operator()(const TensorType & tensor, const IndexType & x) const + { + StencilType stencil; + StencilOffsetsType offsets; + + // Diffusion tensors are homogeneous to the inverse of norms, and are thus rescaled with an inverse spacing. + + TensorType D; + for (ImageDimensionType i = 0; i < ImageDimension; ++i) + for (ImageDimensionType j = i; j < ImageDimension; ++j) + D(i, j) = tensor(i, j) * this->invSpacing[i] * this->invSpacing[j]; + this->Stencil(Dispatch(), D, offsets, stencil.second); + + InternalSizeT * yIndex = &stencil.first[0]; + + // Compute buffer offsets from geometrical offsets + for (unsigned int i = 0; i < HalfStencilSize; ++i) + { + for (auto orientation = 0; orientation < 2; ++orientation, ++yIndex) + { + const IndexType y = orientation ? x - offsets[i] : x + offsets[i]; + if (this->region.IsInside(y)) + { + *yIndex = this->BufferIndex(y); + } + else + { + // Neumann boundary conditions. + *yIndex = this->OutsideBufferIndex(); + } // if y + } // for eps + } // for i + return stencil; + } + +private: + struct DispatchBase + {}; + template + struct Dispatch : public DispatchBase + {}; + + static void + Stencil(const Dispatch<2> &, + const TensorType & D, + StencilOffsetsType & offsets, + StencilCoefficientsType & coefficients) + { + // Construct a superbase, and make it obtuse with Selling's algorithm + VectorType sb[ImageDimension + 1]; // SuperBase + for (ImageDimensionType i = 0; i < ImageDimension; ++i) + { + for (ImageDimensionType j = 0; j < ImageDimension; ++j) + { + sb[i][j] = (i == j); + } + } + + sb[ImageDimension] = -(sb[0] + sb[1]); + constexpr int maxIter = 200; + int iter = 0; + for (; iter < maxIter; ++iter) + { + bool same = true; + for (ImageDimensionType i = 1; i <= ImageDimension && same; ++i) + { + for (ImageDimensionType j = 0; j < i && same; ++j) + { + if (ScalarProduct(D, sb[i], sb[j]) > 0) + { + const VectorType u = sb[i], v = sb[j]; + sb[0] = v - u; + sb[1] = u; + sb[2] = -v; + same = false; + } + } + } + if (same) + { + break; + } + } + if (iter == maxIter) + { + std::cerr << "Warning: Selling's algorithm not stabilized." << std::endl; + } + + for (ImageDimensionType i = 0; i < 3; ++i) + { + coefficients[i] = (-0.5) * ScalarProduct(D, sb[(i + 1) % 3], sb[(i + 2) % 3]); + assert(coefficients[i] >= 0); + + offsets[i][0] = static_cast(-sb[i][1]); + offsets[i][1] = static_cast(sb[i][0]); + } + } + static void + Stencil(const Dispatch<3> &, + const TensorType & D, + StencilOffsetsType & offsets, + StencilCoefficientsType & coefficients) + { + // Construct a superbase, and make it obtuse with Selling's algorithm + VectorType sb[ImageDimension + 1]; + for (ImageDimensionType i = 0; i < ImageDimension; ++i) + { + for (ImageDimensionType j = 0; j < ImageDimension; ++j) + { + sb[i][j] = (i == j); + } + } + sb[ImageDimension] = -(sb[0] + sb[1] + sb[2]); + + constexpr int maxIter = 200; + int iter = 0; + for (; iter < maxIter; ++iter) + { + bool same = true; + for (ImageDimensionType i = 1; i <= ImageDimension && same; ++i) + for (ImageDimensionType j = 0; j < i && same; ++j) + if (ScalarProduct(D, sb[i], sb[j]) > 0) + { + const VectorType u = sb[i], v = sb[j]; + for (ImageDimensionType k = 0, l = 0; k <= ImageDimension; ++k) + if (k != i && k != j) + sb[l++] = sb[k] + u; + sb[2] = -u; + sb[3] = v; + same = false; + } + if (same) + break; + } + if (iter == maxIter) + { + std::cerr << "Warning: Selling's algorithm not stabilized." << std::endl; + } + + // Computation of the weights + SymmetricSecondRankTensor Weights; + for (ImageDimensionType i = 1; i < ImageDimension + 1; ++i) + { + for (ImageDimensionType j = 0; j < i; ++j) + { + Weights(i, j) = (-0.5) * ScalarProduct(D, sb[i], sb[j]); + } + } + + // Now that the obtuse superbasis has been created, generate the stencil. + // First get the dual basis. Obtained by computing the comatrix of Basis[1..ImageDimension]. + + for (ImageDimensionType i = 0; i < ImageDimension; ++i) + { + for (ImageDimensionType j = 0; j < ImageDimension; ++j) + { + offsets[i][j] = sb[(i + 1) % ImageDimension][(j + 1) % ImageDimension] * + sb[(i + 2) % ImageDimension][(j + 2) % ImageDimension] - + sb[(i + 2) % ImageDimension][(j + 1) % ImageDimension] * + sb[(i + 1) % ImageDimension][(j + 2) % ImageDimension]; + } + } + + offsets[ImageDimension] = offsets[0] - offsets[1]; + offsets[ImageDimension + 1] = offsets[0] - offsets[2]; + offsets[ImageDimension + 2] = offsets[1] - offsets[2]; + + // The corresponding coefficients are given by the scalar products. + for (ImageDimensionType i = 0; i < ImageDimension; ++i) + { + coefficients[i] = Weights(i, 3); + } + + coefficients[ImageDimension] = Weights(0, 1); + coefficients[ImageDimension + 1] = Weights(0, 2); + coefficients[ImageDimension + 2] = Weights(1, 2); + } + + RegionType region; + IndexType prod; + SpacingType invSpacing; + InternalSizeT + OutsideBufferIndex() const + { + return NumericTraits::max(); + } +}; + + +template +void +LinearAnisotropicDiffusionLBRImageFilter::GenerateStencils() +{ + // Stencil type is a pair type because itk::UnaryFunctorImage filter + // only produces one output + // using SpacingType = typename TensorImageType::SpacingType; + const RegionType region = GetRequestedRegion(); + + using FunctorFilterType = UnaryFunctorWithIndexImageFilter; + typename FunctorFilterType::Pointer filter = FunctorFilterType::New(); + filter->SetInput(GetInputTensor()); + filter->GetFunctor().Initialize(region, GetInputTensor()->GetSpacing()); + filter->Update(); + m_StencilImage = filter->GetOutput(); + + + // setup diagonal coefficients. Cannot be parallelized due to non-local modifications of diagBuffer. + + m_DiagonalCoefficients = ScalarImageType::New(); + m_DiagonalCoefficients->CopyInformation(GetInputTensor()); + m_DiagonalCoefficients->SetRegions(GetRequestedRegion()); + m_DiagonalCoefficients->Allocate(); + m_DiagonalCoefficients->FillBuffer(ScalarType(0)); + + ImageRegionConstIterator stencilIt(m_StencilImage, region); + ImageRegionIterator diagIt(m_DiagonalCoefficients, region); + ScalarType * diagBuffer = m_DiagonalCoefficients->GetBufferPointer(); + + for (stencilIt.GoToBegin(), diagIt.GoToBegin(); !stencilIt.IsAtEnd(); ++stencilIt, ++diagIt) + { + for (unsigned int i = 0; i < StencilSize; ++i) + { + const InternalSizeT yIndex = stencilIt.Value().first[i]; + if (yIndex != OutsideBufferIndex()) + { + const ScalarType coefficient = stencilIt.Value().second[i / 2]; + diagIt.Value() += coefficient; + diagBuffer[yIndex] += coefficient; + } // if y + } // for i + } // for stencilIt, diagIt +} + + +template +typename LinearAnisotropicDiffusionLBRImageFilter::ScalarType +LinearAnisotropicDiffusionLBRImageFilter::MaxStableTimeStep() +{ + using MaxCalculatorType = MinimumMaximumImageCalculator; + typename MaxCalculatorType::Pointer maximumCalculator = MaxCalculatorType::New(); + maximumCalculator->SetImage(m_DiagonalCoefficients); + maximumCalculator->SetRegion(GetRequestedRegion()); + maximumCalculator->ComputeMaximum(); + return 1. / maximumCalculator->GetMaximum(); +} + + +template +void +LinearAnisotropicDiffusionLBRImageFilter::SetMaxDiffusionTime(ScalarType time) +{ + if (time < 0) + { + itkExceptionMacro("diffusion time must be finite and positive"); + } + m_DiffusionTime = time; +} + + +template +void +LinearAnisotropicDiffusionLBRImageFilter::SetRatioToMaxStableTimeStep(ScalarType ratio) +{ + if (ratio <= 0 || ratio > 1) + { + itkExceptionMacro("Ratio to max time step " << ratio << "should be within ]0,1]"); + } + m_RatioToMaxStableTimeStep = ratio; +} + + +template +void +LinearAnisotropicDiffusionLBRImageFilter::SetMaxNumberOfTimeSteps(int n) +{ + if (n <= 0) + { + itkExceptionMacro("Max number of time steps must be positive"); + } + m_MaxNumberOfTimeSteps = n; +} + + +template +void +LinearAnisotropicDiffusionLBRImageFilter::ImageUpdateLoop() +{ + ScalarType delta = MaxStableTimeStep() * m_RatioToMaxStableTimeStep; + int n = ceil(m_DiffusionTime / delta); + if (n > m_MaxNumberOfTimeSteps) + { + n = m_MaxNumberOfTimeSteps; + m_EffectiveDiffusionTime = n * delta; + } + else + { + delta = m_DiffusionTime / n; + m_EffectiveDiffusionTime = m_DiffusionTime; + } + m_EffectiveNumberOfTimeSteps = n; + + // Extraction of the region of interest is required for image buffer access. + using InputCasterType = ExtractImageFilter; + typename InputCasterType::Pointer inputCaster = InputCasterType::New(); + inputCaster->SetInput(GetInputImage()); + inputCaster->SetExtractionRegion(GetRequestedRegion()); + inputCaster->SetDirectionCollapseToIdentity(); + inputCaster->Update(); + m_PreviousImage = inputCaster->GetOutput(); + + m_NextImage = ImageType::New(); + m_NextImage->CopyInformation(m_PreviousImage); + m_NextImage->SetRegions(m_PreviousImage->GetBufferedRegion()); + m_NextImage->Allocate(); + + for (auto k = 0; k < n; ++k) + { + ImageUpdate(delta); + std::swap(m_PreviousImage, m_NextImage); + this->UpdateProgress(0.5 + 0.5 * k / float(n)); + } + this->GraftOutput(m_PreviousImage); +} + + +template +struct LinearAnisotropicDiffusionLBRImageFilter::FunctorType +{ + ScalarType delta; + PixelType + operator()(PixelType output, PixelType input, ScalarType diag) + { + return output * this->delta + input * (ScalarType(1) - this->delta * diag); + } +}; + + +template +void +LinearAnisotropicDiffusionLBRImageFilter::ImageUpdate(ScalarType delta) +{ + // Setting up iterators + ImageRegion region = GetRequestedRegion(); + + ImageRegionConstIterator inputIt(m_PreviousImage, region); + ImageRegionIterator outputIt(m_NextImage, region); + + const PixelType * inputBuffer = m_PreviousImage->GetBufferPointer(); + PixelType * outputBuffer = m_NextImage->GetBufferPointer(); + + ImageRegionConstIterator diagIt(m_DiagonalCoefficients, region); + ImageRegionConstIterator stencilIt(m_StencilImage, region); + + // Rest of function is a hand-made (sparse matrix)*vector product. + m_NextImage->FillBuffer({}); + + // Taking care of Off-Diagonal matrix elements. Cannot be parallelized due to non-local modifications of outputBuffer + for (inputIt.GoToBegin(), outputIt.GoToBegin(), stencilIt.GoToBegin(); !inputIt.IsAtEnd(); + ++inputIt, ++outputIt, ++stencilIt) + { + for (unsigned int i = 0; i < StencilSize; ++i) + { + const InternalSizeT yIndex = stencilIt.Value().first[i]; + if (yIndex != OutsideBufferIndex()) + { + const ScalarType coefficient = stencilIt.Value().second[i / 2]; + outputIt.Value() += coefficient * (inputBuffer[yIndex]); + outputBuffer[yIndex] += coefficient * inputIt.Value(); + } + } + } + + using ImageFunctorType = TernaryFunctorImageFilter; + typename ImageFunctorType::Pointer imageFunctor = ImageFunctorType::New(); + imageFunctor->SetInput1(m_NextImage); + imageFunctor->SetInput2(m_PreviousImage); + imageFunctor->SetInput3(m_DiagonalCoefficients); + imageFunctor->GetFunctor().delta = delta; + + assert(imageFunctor->CanRunInPlace()); + imageFunctor->InPlaceOn(); + imageFunctor->Update(); + m_NextImage = imageFunctor->GetOutput(); + + /* + // Old Serial version for diagonal elements + for(inputIt.GoToBegin(), outputIt.GoToBegin(), diagIt.GoToBegin(); + !inputIt.IsAtEnd(); + ++inputIt, ++outputIt, ++diagIt) + outputIt.Value() = delta*outputIt.Value() + (1-delta*diagIt.Value())*inputIt.Value(); + */ +} + +// **************************** subclass SSRT_Traits call method ************************** + +template +TScalar +LinearAnisotropicDiffusionLBRImageFilter::ScalarProduct(const TensorType & m, + const VectorType & u, + const VectorType & v) +{ + ScalarType result(0); + for (ImageDimensionType i = 0; i < ImageDimension; ++i) + { + result += m(i, i) * u[i] * v[i]; + } + for (ImageDimensionType i = 0; i < ImageDimension; ++i) + { + for (ImageDimensionType j = i + 1; j < ImageDimension; ++j) + { + result += m(i, j) * (u[i] * v[j] + u[j] * v[i]); + } + } + return result; +} + +} // end namespace itk + +#endif diff --git a/Modules/Filtering/AnisotropicDiffusionLBR/include/itkStructureTensorImageFilter.h b/Modules/Filtering/AnisotropicDiffusionLBR/include/itkStructureTensorImageFilter.h new file mode 100644 index 00000000000..103a574eac9 --- /dev/null +++ b/Modules/Filtering/AnisotropicDiffusionLBR/include/itkStructureTensorImageFilter.h @@ -0,0 +1,158 @@ +/*========================================================================= + * + * 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. + * + *=========================================================================*/ + +// +// Created by Jean-Marie Mirebeau on 05/03/2014. +// +// + +#ifndef itkStructureTensorImageFilter_h +#define itkStructureTensorImageFilter_h + +#include "itkCastImageFilter.h" +#include "itkGradientRecursiveGaussianImageFilter.h" +#include "itkAddImageFilter.h" +#include "itkVectorIndexSelectionCastImageFilter.h" +#include "itkGradientImageFilter.h" +#include "itkSymmetricSecondRankTensor.h" + +namespace itk +{ +/** + * \class StructureTensorImageFilter + * + * \brief Computes the structure tensor. + * + * Implementation of the structure tensor, defined by + * + * \f[K_\rho (\nabla u_\sigma \otimes \nabla u_\sigma),\f] + * + * where \f$K_\rho\f$ denotes the gaussian kernel of standard deviation \f$\rho\f$, + * and \f$u_\sigma := K_\sigma * u\f$. + * + * \ingroup AnisotropicDiffusionLBR + */ +template , + TImage::ImageDimension>> +class StructureTensorImageFilter : public ImageToImageFilter +{ +public: + ITK_DISALLOW_COPY_AND_MOVE(StructureTensorImageFilter); + + using Self = StructureTensorImageFilter; + using Superclass = ImageToImageFilter; + using Pointer = SmartPointer; + using ConstPointer = SmartPointer; + + /// Method for creation through the object factory. + itkNewMacro(Self); + /// Run-time type information (and related methods). + itkOverrideGetNameOfClassMacro(StructureTensorImageFilter); + + using InputImageDimensionType = typename Superclass::InputImageType::ImageDimensionType; + static constexpr InputImageDimensionType InputImageDimension = Superclass::InputImageType::ImageDimension; + + using ImageType = TImage; + using PixelType = typename ImageType::PixelType; + using TensorImageType = TTensorImage; + using TensorType = typename TensorImageType::PixelType; + using ScalarType = typename TensorType::ComponentType; + using ScalarImageType = Image; + + /// Parameter \f$\sigma\f$ of the structure tensor definition. + itkSetMacro(NoiseScale, ScalarType); + /// Parameter \f$\rho\f$ of the structure tensor definition. + itkSetMacro(FeatureScale, ScalarType); + /// Rescales all structure tensors by a common factor, so that the maximum trace is 1. + itkSetMacro(RescaleForUnitMaximumTrace, bool); + + itkGetConstMacro(NoiseScale, ScalarType); + itkGetConstMacro(FeatureScale, ScalarType); + itkGetConstMacro(RescaleForUnitMaximumTrace, bool); + itkGetConstMacro(PostRescaling, ScalarType); /// Global rescaling constant used. + +protected: + void + GenerateData() override; + + ScalarType m_FeatureScale; + ScalarType m_NoiseScale; + bool m_RescaleForUnitMaximumTrace{ false }; + ScalarType m_PostRescaling; + bool m_UseGradientRecursiveGaussianImageFilter{ true }; + + struct DispatchBase + {}; + template + struct Dispatch : public DispatchBase + {}; + + void + IntermediateFilter(const Dispatch &); + void + IntermediateFilter(const Dispatch &); + typename TensorImageType::Pointer m_IntermediateResult; + + using CovariantVectorType = CovariantVector; + using CovariantImageType = Image; + + struct OuterFunctor + { + TensorType + operator()(const CovariantVectorType & u) const + { + TensorType m; + for (InputImageDimensionType i = 0; i < InputImageDimension; ++i) + { + for (InputImageDimensionType j = i; j < InputImageDimension; ++j) + { + m(i, j) = u[i] * u[j]; + } + } + return m; + } + }; + struct TraceFunctor + { + ScalarType + operator()(const TensorType & t) const + { + return t.GetTrace(); + } + }; + struct ScaleFunctor + { + ScalarType scaling; + TensorType + operator()(const TensorType & t) const + { + return t * scaling; + } + }; + + StructureTensorImageFilter(); +}; + +} // end namespace itk + +#ifndef ITK_MANUAL_INSTANTIATION +# include "itkStructureTensorImageFilter.hxx" +#endif + +#endif diff --git a/Modules/Filtering/AnisotropicDiffusionLBR/include/itkStructureTensorImageFilter.hxx b/Modules/Filtering/AnisotropicDiffusionLBR/include/itkStructureTensorImageFilter.hxx new file mode 100644 index 00000000000..82e8d640edb --- /dev/null +++ b/Modules/Filtering/AnisotropicDiffusionLBR/include/itkStructureTensorImageFilter.hxx @@ -0,0 +1,162 @@ +/*========================================================================= + * + * 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. + * + *=========================================================================*/ + +// +// Created by Jean-Marie Mirebeau on 21/11/2014. +// +// + +#ifndef itkStructureTensorImageFilter_hxx +#define itkStructureTensorImageFilter_hxx + +#include "itkMinimumMaximumImageCalculator.h" + +namespace itk +{ + +template +StructureTensorImageFilter::StructureTensorImageFilter() + : m_FeatureScale(2) + , m_NoiseScale(1) + +{} + + +template +void +StructureTensorImageFilter::IntermediateFilter(const Dispatch &) +{ + using GradientFilterType = GradientRecursiveGaussianImageFilter; + typename GradientFilterType::Pointer gradientFilter = GradientFilterType::New(); + gradientFilter->SetInput(this->GetInput()); + gradientFilter->SetSigma(this->m_NoiseScale); + + using OuterFilterType = UnaryFunctorImageFilter; + typename OuterFilterType::Pointer outerFilter = OuterFilterType::New(); + outerFilter->SetInput(gradientFilter->GetOutput()); + + outerFilter->Update(); + this->m_IntermediateResult = outerFilter->GetOutput(); +} + + +template +void +StructureTensorImageFilter::IntermediateFilter(const Dispatch &) +{ + typename Self::ImageType::ConstPointer input = this->GetInput(); + typename Self::TensorImageType::Pointer output = Self::TensorImageType::New(); + output->CopyInformation(input); + output->SetRegions(input->GetRequestedRegion()); + output->Allocate(); + output->FillBuffer(Self::TensorType(0.)); + + for (unsigned int index = 0; index < Self::PixelType::Dimension; ++index) + { + using SelectionFilterType = VectorIndexSelectionCastImageFilter; + typename SelectionFilterType::Pointer selectionFilter = SelectionFilterType::New(); + selectionFilter->SetIndex(index); + selectionFilter->SetInput(input); + + using GaussianFilterType = RecursiveGaussianImageFilter; + using GradientFilterType = + GradientImageFilter; + using GradientGaussianFilterType = + GradientRecursiveGaussianImageFilter; + + typename GaussianFilterType::Pointer gaussianFilter = GaussianFilterType::New(); + typename GradientFilterType::Pointer gradientFilter = GradientFilterType::New(); + typename GradientGaussianFilterType::Pointer gradientGaussianFilter = GradientGaussianFilterType::New(); + + gaussianFilter->SetSigma(this->m_NoiseScale); + gradientGaussianFilter->SetSigma(this->m_NoiseScale); + + using OuterFilterType = + UnaryFunctorImageFilter; + typename OuterFilterType::Pointer outerFilter = OuterFilterType::New(); + + if (this->m_UseGradientRecursiveGaussianImageFilter) + { + gradientGaussianFilter->SetInput(selectionFilter->GetOutput()); + outerFilter->SetInput(gradientGaussianFilter->GetOutput()); + } + else + { + gaussianFilter->SetInput(selectionFilter->GetOutput()); + gradientFilter->SetInput(gaussianFilter->GetOutput()); + outerFilter->SetInput(gradientFilter->GetOutput()); + } + + using AddFilterType = AddImageFilter; + typename AddFilterType::Pointer addFilter = AddFilterType::New(); + addFilter->InPlaceOn(); + addFilter->SetInput1(output); + addFilter->SetInput2(outerFilter->GetOutput()); + addFilter->Update(); + output = addFilter->GetOutput(); + + this->UpdateProgress(index / float(Self::PixelType::Dimension + 1)); + } + this->m_IntermediateResult = output; +} + + +template +void +StructureTensorImageFilter::GenerateData() +{ + this->IntermediateFilter(Dispatch::is_specialized>()); + + using GaussianFilterType = RecursiveGaussianImageFilter; + typename GaussianFilterType::Pointer gaussianFilter = GaussianFilterType::New(); + gaussianFilter->SetInput(m_IntermediateResult); + gaussianFilter->SetSigma(m_FeatureScale); + + if (!m_RescaleForUnitMaximumTrace) + { + m_PostRescaling = 1.; + gaussianFilter->Update(); + this->GraftOutput(gaussianFilter->GetOutput()); + return; + } + + // *** Rescaling for normalization of largest trace *** + + using TraceFilterType = UnaryFunctorImageFilter; + typename TraceFilterType::Pointer traceFilter = TraceFilterType::New(); + traceFilter->SetInput(gaussianFilter->GetOutput()); + + using MaximumCalculatorType = MinimumMaximumImageCalculator; + typename MaximumCalculatorType::Pointer maximumCalculator = MaximumCalculatorType::New(); + maximumCalculator->SetImage(traceFilter->GetOutput()); + + using ScaleFilterType = UnaryFunctorImageFilter; + typename ScaleFilterType::Pointer scaleFilter = ScaleFilterType::New(); + scaleFilter->SetInput(gaussianFilter->GetOutput()); + + traceFilter->Update(); + maximumCalculator->ComputeMaximum(); + m_PostRescaling = 1. / maximumCalculator->GetMaximum(); + scaleFilter->GetFunctor().scaling = m_PostRescaling; + scaleFilter->Update(); + this->GraftOutput(scaleFilter->GetOutput()); +} + +} // end namespace itk + +#endif diff --git a/Modules/Filtering/AnisotropicDiffusionLBR/include/itkUnaryFunctorWithIndexImageFilter.h b/Modules/Filtering/AnisotropicDiffusionLBR/include/itkUnaryFunctorWithIndexImageFilter.h new file mode 100644 index 00000000000..305baf25273 --- /dev/null +++ b/Modules/Filtering/AnisotropicDiffusionLBR/include/itkUnaryFunctorWithIndexImageFilter.h @@ -0,0 +1,98 @@ +/*========================================================================= + * + * 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. + * + *=========================================================================*/ + +// +// Created by Jean-Marie Mirebeau on 07/03/2014. +// +// + +#ifndef itkUnaryFunctorWithIndexImageFilter_h +#define itkUnaryFunctorWithIndexImageFilter_h + +#include "itkImageRegionConstIteratorWithIndex.h" +#include "itkImageScanlineIterator.h" +#include "itkAnisotropicDiffusionLBRMacro.h" +#include "itkImageToImageFilter.h" + +namespace itk +{ + +/** \class UnaryFunctorWithIndexImageFilter + * + * \brief A (simplification of) UnaryFunctorImageFilter, which provides the pixel index to the functor. + * + * \ingroup AnisotropicDiffusionLBR + */ +template +class UnaryFunctorWithIndexImageFilter : public ImageToImageFilter +{ +public: + ITK_DISALLOW_COPY_AND_MOVE(UnaryFunctorWithIndexImageFilter); + + using Self = UnaryFunctorWithIndexImageFilter; + using Superclass = ImageToImageFilter; + using Pointer = SmartPointer; + using ConstPointer = SmartPointer; + + /** Method for creation through the object factory. */ + itkNewMacro(Self); + + /** Run-time type information (and related methods). */ + itkOverrideGetNameOfClassMacro(UnaryFunctorWithIndexImageFilter); + + using InputImageType = TInputImage; + using OutputImageType = TOutputImage; + using FunctorType = TFunctor; + + GetSetFunctorMacro(Functor, FunctorType); + +protected: + UnaryFunctorWithIndexImageFilter() { this->DynamicMultiThreadingOn(); } + + FunctorType m_Functor; + + using InputRegionType = typename InputImageType::RegionType; + using OutputRegionType = typename OutputImageType::RegionType; + + void + DynamicThreadedGenerateData(const OutputRegionType & region) override + { + if (region.GetSize()[0] == 0) + { + return; + } + + ImageRegionConstIteratorWithIndex inputIt(this->GetInput(), region); + ImageScanlineIterator outputIt(this->GetOutput(), region); + + for (inputIt.GoToBegin(), outputIt.GoToBegin(); !outputIt.IsAtEnd();) + { + while (!outputIt.IsAtEndOfLine()) + { + outputIt.Set(m_Functor(inputIt.Value(), inputIt.GetIndex())); + ++inputIt; + ++outputIt; + } + outputIt.NextLine(); + } + } +}; + +} // end namespace itk + +#endif diff --git a/Modules/Filtering/AnisotropicDiffusionLBR/itk-module.cmake b/Modules/Filtering/AnisotropicDiffusionLBR/itk-module.cmake new file mode 100644 index 00000000000..b8f62526703 --- /dev/null +++ b/Modules/Filtering/AnisotropicDiffusionLBR/itk-module.cmake @@ -0,0 +1,23 @@ +set( + DOCUMENTATION + "This module provides coherence-enhancing (CED) and edge-enhancing +(EED) anisotropic diffusion filters built on the +Lattice-Basis-Reduction (LBR) stencil scheme by Jean-Marie Mirebeau. +See the Doxygen on \\\\ref AnisotropicDiffusionLBRImageFilter and +\\\\ref CoherenceEnhancingDiffusionImageFilter for the algorithm and +citations, and the module README for in-tree vs archived-upstream +scope." +) + +itk_module( + AnisotropicDiffusionLBR + DEPENDS + ITKCommon + ITKIOImageBase + ITKImageGradient + TEST_DEPENDS + ITKTestKernel + DESCRIPTION "${DOCUMENTATION}" + EXCLUDE_FROM_DEFAULT + # Header only library, no ENABLE_SHARED +) diff --git a/Modules/Filtering/AnisotropicDiffusionLBR/test/Baseline/Cos3D_cCED.vtk.md5 b/Modules/Filtering/AnisotropicDiffusionLBR/test/Baseline/Cos3D_cCED.vtk.md5 new file mode 100644 index 00000000000..3df854dded2 --- /dev/null +++ b/Modules/Filtering/AnisotropicDiffusionLBR/test/Baseline/Cos3D_cCED.vtk.md5 @@ -0,0 +1 @@ +df95fdb0657f7f8472bdc16c73c5bed0 \ No newline at end of file diff --git a/Modules/Filtering/AnisotropicDiffusionLBR/test/Baseline/FingerPrint_I_20.png.md5 b/Modules/Filtering/AnisotropicDiffusionLBR/test/Baseline/FingerPrint_I_20.png.md5 new file mode 100644 index 00000000000..3b5dbe8354c --- /dev/null +++ b/Modules/Filtering/AnisotropicDiffusionLBR/test/Baseline/FingerPrint_I_20.png.md5 @@ -0,0 +1 @@ +e4e5e233b434ea4c85059d7c62f15554 \ No newline at end of file diff --git a/Modules/Filtering/AnisotropicDiffusionLBR/test/Baseline/FingerPrint_cCED_20.png.md5 b/Modules/Filtering/AnisotropicDiffusionLBR/test/Baseline/FingerPrint_cCED_20.png.md5 new file mode 100644 index 00000000000..750c29db86b --- /dev/null +++ b/Modules/Filtering/AnisotropicDiffusionLBR/test/Baseline/FingerPrint_cCED_20.png.md5 @@ -0,0 +1 @@ +07435f1d44aeb66fd98e642945437662 \ No newline at end of file diff --git a/Modules/Filtering/AnisotropicDiffusionLBR/test/Baseline/FingerPrint_cEED_20.png.md5 b/Modules/Filtering/AnisotropicDiffusionLBR/test/Baseline/FingerPrint_cEED_20.png.md5 new file mode 100644 index 00000000000..76c2e516e80 --- /dev/null +++ b/Modules/Filtering/AnisotropicDiffusionLBR/test/Baseline/FingerPrint_cEED_20.png.md5 @@ -0,0 +1 @@ +5e752e3fa0e46bb530fedc94c7794c73 \ No newline at end of file diff --git a/Modules/Filtering/AnisotropicDiffusionLBR/test/Baseline/Lena_Detail_I_2.png.md5 b/Modules/Filtering/AnisotropicDiffusionLBR/test/Baseline/Lena_Detail_I_2.png.md5 new file mode 100644 index 00000000000..b80ed7eb732 --- /dev/null +++ b/Modules/Filtering/AnisotropicDiffusionLBR/test/Baseline/Lena_Detail_I_2.png.md5 @@ -0,0 +1 @@ +ff36663855e6794712b081689aac70e5 \ No newline at end of file diff --git a/Modules/Filtering/AnisotropicDiffusionLBR/test/Baseline/Lena_Detail_cCED_2.png.md5 b/Modules/Filtering/AnisotropicDiffusionLBR/test/Baseline/Lena_Detail_cCED_2.png.md5 new file mode 100644 index 00000000000..f6a2c0d79ae --- /dev/null +++ b/Modules/Filtering/AnisotropicDiffusionLBR/test/Baseline/Lena_Detail_cCED_2.png.md5 @@ -0,0 +1 @@ +b741b80ce65e20c59f286244f621344a \ No newline at end of file diff --git a/Modules/Filtering/AnisotropicDiffusionLBR/test/Baseline/Lena_Detail_cEED_2.png.md5 b/Modules/Filtering/AnisotropicDiffusionLBR/test/Baseline/Lena_Detail_cEED_2.png.md5 new file mode 100644 index 00000000000..d529c7f1aa4 --- /dev/null +++ b/Modules/Filtering/AnisotropicDiffusionLBR/test/Baseline/Lena_Detail_cEED_2.png.md5 @@ -0,0 +1 @@ +c606fea9c82019d1b4e80d351b803d92 \ No newline at end of file diff --git a/Modules/Filtering/AnisotropicDiffusionLBR/test/Baseline/Oscillations1_CED.png.md5 b/Modules/Filtering/AnisotropicDiffusionLBR/test/Baseline/Oscillations1_CED.png.md5 new file mode 100644 index 00000000000..b5cb7d21602 --- /dev/null +++ b/Modules/Filtering/AnisotropicDiffusionLBR/test/Baseline/Oscillations1_CED.png.md5 @@ -0,0 +1 @@ +50da37ff706c93536c0f33390da4287c \ No newline at end of file diff --git a/Modules/Filtering/AnisotropicDiffusionLBR/test/Baseline/Oscillations1_cCED.png.md5 b/Modules/Filtering/AnisotropicDiffusionLBR/test/Baseline/Oscillations1_cCED.png.md5 new file mode 100644 index 00000000000..35dfbe815ba --- /dev/null +++ b/Modules/Filtering/AnisotropicDiffusionLBR/test/Baseline/Oscillations1_cCED.png.md5 @@ -0,0 +1 @@ +21e83dc09f4c58a44eeb676e49ec3d99 \ No newline at end of file diff --git a/Modules/Filtering/AnisotropicDiffusionLBR/test/Baseline/PacMan_I.png.md5 b/Modules/Filtering/AnisotropicDiffusionLBR/test/Baseline/PacMan_I.png.md5 new file mode 100644 index 00000000000..6a1cf5d0156 --- /dev/null +++ b/Modules/Filtering/AnisotropicDiffusionLBR/test/Baseline/PacMan_I.png.md5 @@ -0,0 +1 @@ +c5b358267defea8babcfebbc66c9fa8b \ No newline at end of file diff --git a/Modules/Filtering/AnisotropicDiffusionLBR/test/Baseline/PacMan_cCED.png.md5 b/Modules/Filtering/AnisotropicDiffusionLBR/test/Baseline/PacMan_cCED.png.md5 new file mode 100644 index 00000000000..223ef256a52 --- /dev/null +++ b/Modules/Filtering/AnisotropicDiffusionLBR/test/Baseline/PacMan_cCED.png.md5 @@ -0,0 +1 @@ +fd8d652016508d93ee861c1db83f3ed4 \ No newline at end of file diff --git a/Modules/Filtering/AnisotropicDiffusionLBR/test/Baseline/PacMan_cEED.png.md5 b/Modules/Filtering/AnisotropicDiffusionLBR/test/Baseline/PacMan_cEED.png.md5 new file mode 100644 index 00000000000..a08c69480fd --- /dev/null +++ b/Modules/Filtering/AnisotropicDiffusionLBR/test/Baseline/PacMan_cEED.png.md5 @@ -0,0 +1 @@ +edf293e2cce2eae1df4f8598e8179641 \ No newline at end of file diff --git a/Modules/Filtering/AnisotropicDiffusionLBR/test/Baseline/Triangle_EED.png.md5 b/Modules/Filtering/AnisotropicDiffusionLBR/test/Baseline/Triangle_EED.png.md5 new file mode 100644 index 00000000000..09f2a9ab526 --- /dev/null +++ b/Modules/Filtering/AnisotropicDiffusionLBR/test/Baseline/Triangle_EED.png.md5 @@ -0,0 +1 @@ +f823b62e9135a37c7438fa07a9e54096 \ No newline at end of file diff --git a/Modules/Filtering/AnisotropicDiffusionLBR/test/Baseline/Triangle_cEED.png.md5 b/Modules/Filtering/AnisotropicDiffusionLBR/test/Baseline/Triangle_cEED.png.md5 new file mode 100644 index 00000000000..ff3bfe4ac73 --- /dev/null +++ b/Modules/Filtering/AnisotropicDiffusionLBR/test/Baseline/Triangle_cEED.png.md5 @@ -0,0 +1 @@ +12473a0cb8d3afa0f8d7eb4f61e6216b \ No newline at end of file diff --git a/Modules/Filtering/AnisotropicDiffusionLBR/test/Baseline/VectorField_Circle_cEED.vtk.md5 b/Modules/Filtering/AnisotropicDiffusionLBR/test/Baseline/VectorField_Circle_cEED.vtk.md5 new file mode 100644 index 00000000000..77458bed992 --- /dev/null +++ b/Modules/Filtering/AnisotropicDiffusionLBR/test/Baseline/VectorField_Circle_cEED.vtk.md5 @@ -0,0 +1 @@ +23c4495de1a746648418abc144972e92 \ No newline at end of file diff --git a/Modules/Filtering/AnisotropicDiffusionLBR/test/Baseline/mrbrain_cEED.vtk.md5 b/Modules/Filtering/AnisotropicDiffusionLBR/test/Baseline/mrbrain_cEED.vtk.md5 new file mode 100644 index 00000000000..9e91116b83f --- /dev/null +++ b/Modules/Filtering/AnisotropicDiffusionLBR/test/Baseline/mrbrain_cEED.vtk.md5 @@ -0,0 +1 @@ +db41c262fba84a75eb399e1e154a5974 \ No newline at end of file diff --git a/Modules/Filtering/AnisotropicDiffusionLBR/test/CMakeLists.txt b/Modules/Filtering/AnisotropicDiffusionLBR/test/CMakeLists.txt new file mode 100644 index 00000000000..12b772071f5 --- /dev/null +++ b/Modules/Filtering/AnisotropicDiffusionLBR/test/CMakeLists.txt @@ -0,0 +1,275 @@ +itk_module_test() + +set(AnisotropicDiffusionLBRTests CoherenceEnhancingDiffusionTest.cxx) + +createtestdriver(AnisotropicDiffusionLBR "${AnisotropicDiffusionLBR-Test_LIBRARIES}" "${AnisotropicDiffusionLBRTests}") + +set(TestOutput ${ITK_TEST_OUTPUT_DIR}) + +################ fig:PacMan ################# + +itk_add_test( + NAME PacMan_cEED + COMMAND + AnisotropicDiffusionLBRTestDriver + --compare + DATA{Baseline/PacMan_cEED.png} + ${TestOutput}/PacMan_cEED.png + CoherenceEnhancingDiffusionTest + DATA{Input/PacMan.png} + ${TestOutput}/PacMan_cEED.png + 20 + 0.05 + cEED + 3 +) + +itk_add_test( + NAME PacMan_cCED + COMMAND + AnisotropicDiffusionLBRTestDriver + --compare + DATA{Baseline/PacMan_cCED.png} + ${TestOutput}/PacMan_cCED.png + CoherenceEnhancingDiffusionTest + DATA{Input/PacMan.png} + ${TestOutput}/PacMan_cCED.png + 20 + 0.05 + cCED + 3 +) + +itk_add_test( + NAME PacMan_I + COMMAND + AnisotropicDiffusionLBRTestDriver + --compare + DATA{Baseline/PacMan_I.png} + ${TestOutput}/PacMan_I.png + CoherenceEnhancingDiffusionTest + DATA{Input/PacMan.png} + ${TestOutput}/PacMan_I.png + 20 + 0.05 + Isotropic + 3 +) + +############# fig : Fingerprint ############# + +itk_add_test( + NAME FingerPrint_cEED + COMMAND + AnisotropicDiffusionLBRTestDriver + --compare + DATA{Baseline/FingerPrint_cEED_20.png} + ${TestOutput}/FingerPrint_cEED_20.png + CoherenceEnhancingDiffusionTest + DATA{Input/FingerPrint.png} + ${TestOutput}/FingerPrint_cEED_20.png + 20 + 0.02 + cEED +) + +itk_add_test( + NAME FingerPrint_cCED + COMMAND + AnisotropicDiffusionLBRTestDriver + --compare + DATA{Baseline/FingerPrint_cCED_20.png} + ${TestOutput}/FingerPrint_cCED_20.png + CoherenceEnhancingDiffusionTest + DATA{Input/FingerPrint.png} + ${TestOutput}/FingerPrint_cCED_20.png + 20 + 0.02 + cCED +) + +itk_add_test( + NAME FingerPrint_I + COMMAND + AnisotropicDiffusionLBRTestDriver + --compare + DATA{Baseline/FingerPrint_I_20.png} + ${TestOutput}/FingerPrint_I_20.png + CoherenceEnhancingDiffusionTest + DATA{Input/FingerPrint.png} + ${TestOutput}/FingerPrint_I_20.png + 20 + 0.02 + Isotropic +) + +################ fig : Cos3D ################ + +itk_add_test( + NAME Cos3D_cCED + COMMAND + AnisotropicDiffusionLBRTestDriver + --compare + DATA{Baseline/Cos3D_cCED.vtk} + ${TestOutput}/Cos3D_cCED.vtk + CoherenceEnhancingDiffusionTest + DATA{Input/Cos3D_Noisy.vtk} + ${TestOutput}/Cos3D_cCED.vtk + 10 + 0.02 + cCED + 4 + 10 +) + +################ fig : Brain ################ + +itk_add_test( + NAME mrbrain_cEED + COMMAND + AnisotropicDiffusionLBRTestDriver + --compare + DATA{Baseline/mrbrain_cEED.vtk} + ${TestOutput}/mrbrain_cEED.vtk + CoherenceEnhancingDiffusionTest + DATA{Input/mrbrain_noisy.vtk} + ${TestOutput}/mrbrain_cEED.vtk + 5 + 0.003 + cEED +) + +################# fig : Lena ################ + +itk_add_test( + NAME Lena_Detail_cCED_2 + COMMAND + AnisotropicDiffusionLBRTestDriver + --compare + DATA{Baseline/Lena_Detail_cCED_2.png} + ${TestOutput}/Lena_Detail_cCED_2.png + CoherenceEnhancingDiffusionTest + DATA{Input/Lena_Detail.png} + ${TestOutput}/Lena_Detail_cCED_2.png + 2 + 0.003 + cCED + 0.5 + 2 + 4 +) + +itk_add_test( + NAME Lena_Detail_cEED_2 + COMMAND + AnisotropicDiffusionLBRTestDriver + --compare + DATA{Baseline/Lena_Detail_cEED_2.png} + ${TestOutput}/Lena_Detail_cEED_2.png + CoherenceEnhancingDiffusionTest + DATA{Input/Lena_Detail.png} + ${TestOutput}/Lena_Detail_cEED_2.png + 2 + 0.003 + cEED + 0.5 + 2 + 4 +) + +itk_add_test( + NAME Lena_Detail_I_2 + COMMAND + AnisotropicDiffusionLBRTestDriver + --compare + DATA{Baseline/Lena_Detail_I_2.png} + ${TestOutput}/Lena_Detail_I_2.png + CoherenceEnhancingDiffusionTest + DATA{Input/Lena_Detail.png} + ${TestOutput}/Lena_Detail_I_2.png + 2 + 0.003 + Isotropic + 0.5 + 2 + 4 +) + +############## fig : Vector field ########### + +itk_add_test( + NAME VectorField_Circle_cEED + COMMAND + AnisotropicDiffusionLBRTestDriver + --compare + DATA{Baseline/VectorField_Circle_cEED.vtk} + ${TestOutput}/VectorField_Circle_cEED.vtk + CoherenceEnhancingDiffusionTest + DATA{Input/VectorField_CircleOpposites.vtk} + ${TestOutput}/VectorField_Circle_cEED.vtk + 10 + 0.05 + cEED +) + +############### fig : Triangle ############## + +itk_add_test( + NAME Triangle_cEED + COMMAND + AnisotropicDiffusionLBRTestDriver + --compare + DATA{Baseline/Triangle_cEED.png} + ${TestOutput}/Triangle_cEED.png + CoherenceEnhancingDiffusionTest + DATA{Input/Triangle.png} + ${TestOutput}/Triangle_cEED.png + 5 + 0.03 + cEED +) + +itk_add_test( + NAME Triangle_EED + COMMAND + AnisotropicDiffusionLBRTestDriver + --compare + DATA{Baseline/Triangle_EED.png} + ${TestOutput}/Triangle_EED.png + CoherenceEnhancingDiffusionTest + DATA{Input/Triangle.png} + ${TestOutput}/Triangle_EED.png + 5 + 0.03 + EED +) + +itk_add_test( + NAME Oscillations1_cCED + COMMAND + AnisotropicDiffusionLBRTestDriver + --compare + DATA{Baseline/Oscillations1_cCED.png} + ${TestOutput}/Oscillations1_cCED.png + CoherenceEnhancingDiffusionTest + DATA{Input/Oscillations_Noisy1.png} + ${TestOutput}/Oscillations1_cCED.png + 20 + 0.05 + cCED +) + +itk_add_test( + NAME Oscillations1_CED + COMMAND + AnisotropicDiffusionLBRTestDriver + --compare + DATA{Baseline/Oscillations1_CED.png} + ${TestOutput}/Oscillations1_CED.png + CoherenceEnhancingDiffusionTest + DATA{Input/Oscillations_Noisy1.png} + ${TestOutput}/Oscillations1_CED.png + 20 + 0.05 + CED +) diff --git a/Modules/Filtering/AnisotropicDiffusionLBR/test/CoherenceEnhancingDiffusionTest.cxx b/Modules/Filtering/AnisotropicDiffusionLBR/test/CoherenceEnhancingDiffusionTest.cxx new file mode 100644 index 00000000000..d2465caacf8 --- /dev/null +++ b/Modules/Filtering/AnisotropicDiffusionLBR/test/CoherenceEnhancingDiffusionTest.cxx @@ -0,0 +1,48 @@ +/*========================================================================= + * + * 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 +#include "CoherenceEnhancingDiffusionCommandLine.h" + +int +CoherenceEnhancingDiffusionTest(int argc, char * argv[]) +{ + try + { + CoherenceEnhancingDiffusionCommandLine::Execute(argc, argv); + } + catch (itk::ExceptionObject & e) + { + std::cerr << "ITK Exception : " << e.GetDescription() << std::endl; + return EXIT_FAILURE; + } + return EXIT_SUCCESS; +} + +/* +int LinearAnisotropicDiffusionTest(int argc, char **argv) +{ + try { + LinearAnisotropicDiffusionCommandLine::Execute(argc, argv); + } catch (itk::ExceptionObject& e) { + std::cerr << "ITK Exception : " << e.GetDescription() << std::endl; + return EXIT_FAILURE; + } + return EXIT_SUCCESS; +}*/ diff --git a/Modules/Filtering/AnisotropicDiffusionLBR/test/Input/Cos3D_Noisy.vtk.md5 b/Modules/Filtering/AnisotropicDiffusionLBR/test/Input/Cos3D_Noisy.vtk.md5 new file mode 100644 index 00000000000..c48093eec6b --- /dev/null +++ b/Modules/Filtering/AnisotropicDiffusionLBR/test/Input/Cos3D_Noisy.vtk.md5 @@ -0,0 +1 @@ +3a7d9131a732794fcb4100909cd3fd1c \ No newline at end of file diff --git a/Modules/Filtering/AnisotropicDiffusionLBR/test/Input/FingerPrint.png.md5 b/Modules/Filtering/AnisotropicDiffusionLBR/test/Input/FingerPrint.png.md5 new file mode 100644 index 00000000000..2992b7a556c --- /dev/null +++ b/Modules/Filtering/AnisotropicDiffusionLBR/test/Input/FingerPrint.png.md5 @@ -0,0 +1 @@ +ed7342b4598d44574b2714834b705cad \ No newline at end of file diff --git a/Modules/Filtering/AnisotropicDiffusionLBR/test/Input/Lena_Detail.png.md5 b/Modules/Filtering/AnisotropicDiffusionLBR/test/Input/Lena_Detail.png.md5 new file mode 100644 index 00000000000..2d8222d7362 --- /dev/null +++ b/Modules/Filtering/AnisotropicDiffusionLBR/test/Input/Lena_Detail.png.md5 @@ -0,0 +1 @@ +45a1845c6fa452c7465bebda5bbe9b0f \ No newline at end of file diff --git a/Modules/Filtering/AnisotropicDiffusionLBR/test/Input/Oscillations_Noisy1.png.md5 b/Modules/Filtering/AnisotropicDiffusionLBR/test/Input/Oscillations_Noisy1.png.md5 new file mode 100644 index 00000000000..4797a352a96 --- /dev/null +++ b/Modules/Filtering/AnisotropicDiffusionLBR/test/Input/Oscillations_Noisy1.png.md5 @@ -0,0 +1 @@ +75818e9d765fb6838a8cf5845ac19b9a \ No newline at end of file diff --git a/Modules/Filtering/AnisotropicDiffusionLBR/test/Input/PacMan.png.md5 b/Modules/Filtering/AnisotropicDiffusionLBR/test/Input/PacMan.png.md5 new file mode 100644 index 00000000000..0abec1e6758 --- /dev/null +++ b/Modules/Filtering/AnisotropicDiffusionLBR/test/Input/PacMan.png.md5 @@ -0,0 +1 @@ +d7955368c6f49cbb451d8901aa40add6 \ No newline at end of file diff --git a/Modules/Filtering/AnisotropicDiffusionLBR/test/Input/Triangle.png.md5 b/Modules/Filtering/AnisotropicDiffusionLBR/test/Input/Triangle.png.md5 new file mode 100644 index 00000000000..c5594b80c04 --- /dev/null +++ b/Modules/Filtering/AnisotropicDiffusionLBR/test/Input/Triangle.png.md5 @@ -0,0 +1 @@ +bce40d3af4f491d728aaba8bb8c9ede9 \ No newline at end of file diff --git a/Modules/Filtering/AnisotropicDiffusionLBR/test/Input/VectorField_CircleOpposites.vtk.md5 b/Modules/Filtering/AnisotropicDiffusionLBR/test/Input/VectorField_CircleOpposites.vtk.md5 new file mode 100644 index 00000000000..e664c2ad407 --- /dev/null +++ b/Modules/Filtering/AnisotropicDiffusionLBR/test/Input/VectorField_CircleOpposites.vtk.md5 @@ -0,0 +1 @@ +0a9e85b2b8dfadb4ab25b828a0f23852 \ No newline at end of file diff --git a/Modules/Filtering/AnisotropicDiffusionLBR/test/Input/mrbrain_noisy.vtk.md5 b/Modules/Filtering/AnisotropicDiffusionLBR/test/Input/mrbrain_noisy.vtk.md5 new file mode 100644 index 00000000000..2e105fba437 --- /dev/null +++ b/Modules/Filtering/AnisotropicDiffusionLBR/test/Input/mrbrain_noisy.vtk.md5 @@ -0,0 +1 @@ +ff88f04e75dc478b283ecdf39d8d7687 \ No newline at end of file diff --git a/Modules/Filtering/AnisotropicDiffusionLBR/wrapping/CMakeLists.txt b/Modules/Filtering/AnisotropicDiffusionLBR/wrapping/CMakeLists.txt new file mode 100644 index 00000000000..cd368a77b2c --- /dev/null +++ b/Modules/Filtering/AnisotropicDiffusionLBR/wrapping/CMakeLists.txt @@ -0,0 +1,3 @@ +itk_wrap_module(AnisotropicDiffusionLBR) +set(WRAPPING_SUBMODULE_ORDER itkAnisotropicDiffusionLBRImageFilter) +itk_auto_load_and_end_wrap_submodules() diff --git a/Modules/Filtering/AnisotropicDiffusionLBR/wrapping/itkAnisotropicDiffusionLBRImageFilter.wrap b/Modules/Filtering/AnisotropicDiffusionLBR/wrapping/itkAnisotropicDiffusionLBRImageFilter.wrap new file mode 100644 index 00000000000..a1c612cfc4f --- /dev/null +++ b/Modules/Filtering/AnisotropicDiffusionLBR/wrapping/itkAnisotropicDiffusionLBRImageFilter.wrap @@ -0,0 +1,3 @@ +itk_wrap_class("itk::AnisotropicDiffusionLBRImageFilter" POINTER) + itk_wrap_image_filter("${WRAP_ITK_SCALAR}" 1 "2;3") +itk_end_wrap_class() diff --git a/Modules/Filtering/AnisotropicDiffusionLBR/wrapping/itkCoherenceEnhancingDiffusionImageFilter.wrap b/Modules/Filtering/AnisotropicDiffusionLBR/wrapping/itkCoherenceEnhancingDiffusionImageFilter.wrap new file mode 100644 index 00000000000..e02e4cfd944 --- /dev/null +++ b/Modules/Filtering/AnisotropicDiffusionLBR/wrapping/itkCoherenceEnhancingDiffusionImageFilter.wrap @@ -0,0 +1,3 @@ +itk_wrap_class("itk::CoherenceEnhancingDiffusionImageFilter" POINTER_WITH_SUPERCLASS) + itk_wrap_image_filter("${WRAP_ITK_SCALAR}" 1 "2;3") +itk_end_wrap_class() diff --git a/Modules/Remote/AnisotropicDiffusionLBR.remote.cmake b/Modules/Remote/AnisotropicDiffusionLBR.remote.cmake deleted file mode 100644 index a019883e2ba..00000000000 --- a/Modules/Remote/AnisotropicDiffusionLBR.remote.cmake +++ /dev/null @@ -1,68 +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( - AnisotropicDiffusionLBR - "Anisotropic Non-Linear Diffusion is a powerful image processing technique, - which allows to simultaneously remove the noise and enhance sharp features - in two or three dimensional images. Anisotropic Diffusion is understood here - in the sense of Weickert, meaning that diffusion tensors are anisotropic and - reflect the local orientation of image features. This is in contrast with - the non-linear diffusion filter of Perona and Malik, which only involves - scalar diffusion coefficients, in other words isotropic diffusion tensors. - - In this module, an anisotropic non-linear diffusion technique based on a recent - adaptive scheme making the diffusion stable and requiring limited numerical resources - is available. - - From the Insight Journal Article: - - \"Anisotropic Diffusion in ITK\" - Mirebeau J., Fehrenbach J., Risser L., Tobji S. - The Insight Journal. 2014 January-December. - ttps://doi.org/10.54294/en3833 - " - MODULE_COMPLIANCE_LEVEL 2 - GIT_REPOSITORY https://github.com/InsightSoftwareConsortium/ITKAnisotropicDiffusionLBR.git - GIT_TAG 753f2498af4dbd4a864fe6378da72e92654b7f80 - ) diff --git a/pyproject.toml b/pyproject.toml index 5333ba058f1..26dc2343519 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -52,6 +52,7 @@ cmd = '''cmake -DITK_USE_CCACHE:BOOL=ON -DCMAKE_C_COMPILER_LAUNCHER:STRING=ccache -DCMAKE_CXX_COMPILER_LAUNCHER:STRING=ccache + -DModule_AnisotropicDiffusionLBR:BOOL=ON -DITK_COMPUTER_MEMORY_SIZE:STRING=11 -DModule_StructuralSimilarity:BOOL=ON''' description = "Configure ITK for CI (with ccache compiler launcher)"