diff --git a/Examples/Filtering/LaplacianSharpeningImageFilter.cxx b/Examples/Filtering/LaplacianSharpeningImageFilter.cxx index 874d7747479..333cfaef067 100644 --- a/Examples/Filtering/LaplacianSharpeningImageFilter.cxx +++ b/Examples/Filtering/LaplacianSharpeningImageFilter.cxx @@ -34,19 +34,21 @@ main(int argc, char * argv[]) const char * inputFilename = argv[1]; const char * outputFilename = argv[2]; + using PixelType = float; using CharPixelType = unsigned char; constexpr unsigned int Dimension = 2; + using ImageType = itk::Image; using CharImageType = itk::Image; - using ReaderType = itk::ImageFileReader; + using ReaderType = itk::ImageFileReader; using WriterType = itk::ImageFileWriter; using RescaleFilter = - itk::RescaleIntensityImageFilter; + itk::RescaleIntensityImageFilter; using LaplacianSharpeningFilter = - itk::LaplacianSharpeningImageFilter; + itk::LaplacianSharpeningImageFilter; // Setting the IO diff --git a/Modules/Filtering/ImageFeature/include/itkLaplacianSharpeningImageFilter.h b/Modules/Filtering/ImageFeature/include/itkLaplacianSharpeningImageFilter.h index 39bc0d7815f..cc9b7a4f552 100644 --- a/Modules/Filtering/ImageFeature/include/itkLaplacianSharpeningImageFilter.h +++ b/Modules/Filtering/ImageFeature/include/itkLaplacianSharpeningImageFilter.h @@ -66,6 +66,7 @@ class ITK_TEMPLATE_EXPORT LaplacianSharpeningImageFilter : public ImageToImageFi using RealType = typename NumericTraits::RealType; using InputPixelType = typename TInputImage::PixelType; using InputInternalPixelType = typename TInputImage::InternalPixelType; + static constexpr unsigned int InputImageDimension = TInputImage::ImageDimension; static constexpr unsigned int ImageDimension = TOutputImage::ImageDimension; /** Image type alias support */ @@ -83,17 +84,6 @@ class ITK_TEMPLATE_EXPORT LaplacianSharpeningImageFilter : public ImageToImageFi /** Method for creation through the object factory. */ itkNewMacro(Self); - /** LaplacianSharpeningImageFilter needs a larger input requested - * region than the output requested region (larger in the direction - * of the derivative). As such, LaplacianSharpeningImageFilter - * needs to provide an implementation for - * GenerateInputRequestedRegion() in order to inform the pipeline - * execution model. - * - * \sa ImageToImageFilter::GenerateInputRequestedRegion() */ - void - GenerateInputRequestedRegion() override; - /** Enable/Disable using the image spacing information in * calculations. Use this option if you want derivatives in * physical space. Default is UseImageSpacingOn. */ @@ -104,11 +94,23 @@ class ITK_TEMPLATE_EXPORT LaplacianSharpeningImageFilter : public ImageToImageFi itkSetMacro(UseImageSpacing, bool); itkGetConstMacro(UseImageSpacing, bool); +#ifdef ITK_USE_CONCEPT_CHECKING + // Begin concept checking + itkConceptMacro(SameDimensionCheck, (Concept::SameDimension)); + itkConceptMacro(InputPixelTypeIsFloatingPointCheck, (Concept::IsFloatingPoint)); + itkConceptMacro(OutputPixelTypeIsFloatingPointCheck, (Concept::IsFloatingPoint)); + // End concept checking +#endif + protected: LaplacianSharpeningImageFilter() { m_UseImageSpacing = true; } ~LaplacianSharpeningImageFilter() override = default; + /** This filter must produce all of its output at once. */ + void + EnlargeOutputRequestedRegion(DataObject * output) override; + /** Standard pipeline method. While this class does not implement a * ThreadedGenerateData(), its GenerateData() delegates all * calculations to an NeighborhoodOperatorImageFilter. Since the diff --git a/Modules/Filtering/ImageFeature/include/itkLaplacianSharpeningImageFilter.hxx b/Modules/Filtering/ImageFeature/include/itkLaplacianSharpeningImageFilter.hxx index 69121d08ad0..f69d564c760 100644 --- a/Modules/Filtering/ImageFeature/include/itkLaplacianSharpeningImageFilter.hxx +++ b/Modules/Filtering/ImageFeature/include/itkLaplacianSharpeningImageFilter.hxx @@ -18,217 +18,130 @@ #ifndef itkLaplacianSharpeningImageFilter_hxx #define itkLaplacianSharpeningImageFilter_hxx -#include "itkNeighborhoodOperatorImageFilter.h" -#include "itkLaplacianOperator.h" -#include "itkProgressAccumulator.h" +#include "itkLaplacianImageFilter.h" #include "itkMinimumMaximumImageCalculator.h" -#include "itkImageRegionIterator.h" +#include "itkStatisticsImageFilter.h" +#include "itkAddImageFilter.h" +#include "itkMultiplyImageFilter.h" +#include "itkSubtractImageFilter.h" +#include "itkIntensityWindowingImageFilter.h" namespace itk { + template void -LaplacianSharpeningImageFilter::PrintSelf(std::ostream & os, Indent indent) const +LaplacianSharpeningImageFilter::EnlargeOutputRequestedRegion(DataObject * output) { - Superclass::PrintSelf(os, indent); + Superclass::EnlargeOutputRequestedRegion(output); - os << indent << "UseImageSpacing: " << (m_UseImageSpacing ? "On" : "Off") << std::endl; + if (output != nullptr) + { + output->SetRequestedRegionToLargestPossibleRegion(); + } } template void -LaplacianSharpeningImageFilter::GenerateInputRequestedRegion() +LaplacianSharpeningImageFilter::GenerateData() { - // call the superclass' implementation of this method. This should - // copy the output requested region to the input requested region - Superclass::GenerateInputRequestedRegion(); - // get pointers to the input and output - InputImagePointer inputPtr = const_cast(this->GetInput()); + auto localInput = TInputImage::New(); + localInput->Graft(this->GetInput()); - if (!inputPtr) - { - return; - } + // Calculate needed input image statistics - // Build an operator so that we can determine the kernel size - LaplacianOperator oper; - oper.CreateOperator(); + typename StatisticsImageFilter::Pointer inputCalculator = + StatisticsImageFilter::New(); - // get a copy of the input requested region (should equal the output - // requested region) - typename TInputImage::RegionType inputRequestedRegion; - inputRequestedRegion = inputPtr->GetRequestedRegion(); + inputCalculator->SetInput(localInput); + inputCalculator->Update(); - // pad the input requested region by the operator radius - inputRequestedRegion.PadByRadius(oper.GetRadius()); + auto inputMinimum = static_cast(inputCalculator->GetMinimum()); + auto inputMaximum = static_cast(inputCalculator->GetMaximum()); + auto inputMean = static_cast(inputCalculator->GetMean()); - // crop the input requested region at the input's largest possible region - if (inputRequestedRegion.Crop(inputPtr->GetLargestPossibleRegion())) - { - inputPtr->SetRequestedRegion(inputRequestedRegion); - return; - } - else - { - // Couldn't crop the region (requested region is outside the largest - // possible region). Throw an exception. - - // store what we tried to request (prior to trying to crop) - inputPtr->SetRequestedRegion(inputRequestedRegion); - - // build an exception - InvalidRequestedRegionError e(__FILE__, __LINE__); - e.SetLocation(ITK_LOCATION); - e.SetDescription("Requested region is (at least partially) outside the largest possible region."); - e.SetDataObject(inputPtr); - throw e; - } -} + auto inputShift = inputMinimum; + auto inputScale = inputMaximum - inputMinimum; -template -void -LaplacianSharpeningImageFilter::GenerateData() -{ - // Create the Laplacian operator - LaplacianOperator oper; - double s[ImageDimension]; - for (unsigned int i = 0; i < ImageDimension; ++i) - { - if (this->GetInput()->GetSpacing()[i] == 0.0) - { - itkExceptionMacro("Image spacing cannot be zero"); - } - else - { - s[i] = 1.0 / this->GetInput()->GetSpacing()[i]; - } - } - oper.SetDerivativeScalings(s); - oper.CreateOperator(); + // Calculate the Laplacian filtered image - // do calculations in floating point using RealImageType = Image; - using NOIF = NeighborhoodOperatorImageFilter; - ZeroFluxNeumannBoundaryCondition nbc; - auto filter = NOIF::New(); - filter->OverrideBoundaryCondition(static_cast(&nbc)); + using LaplacianImageFilter = LaplacianImageFilter; + typename LaplacianImageFilter::Pointer laplacianFilter = LaplacianImageFilter::New(); + laplacianFilter->SetInput(localInput); + laplacianFilter->SetUseImageSpacing(m_UseImageSpacing); + laplacianFilter->Update(); - // Create a process accumulator for tracking the progress of this minipipeline - auto progress = ProgressAccumulator::New(); - progress->SetMiniPipelineFilter(this); + // Calculate needed laplacian filtered image statistics - // Register the filter with the with progress accumulator using - // equal weight proportion - progress->RegisterInternalFilter(filter, 0.8f); + typename MinimumMaximumImageCalculator::Pointer filteredCalculator = + MinimumMaximumImageCalculator::New(); - // - // set up the mini-pipeline - // - filter->SetOperator(oper); - filter->SetInput(this->GetInput()); - filter->GetOutput()->SetRequestedRegion(this->GetOutput()->GetRequestedRegion()); + filteredCalculator->SetImage(laplacianFilter->GetOutput()); + filteredCalculator->SetRegion(laplacianFilter->GetOutput()->GetRequestedRegion()); + filteredCalculator->Compute(); - // execute the mini-pipeline - filter->Update(); + RealType filteredShift = static_cast(filteredCalculator->GetMinimum()); + RealType filteredScale = static_cast(filteredCalculator->GetMaximum()) - filteredShift; - // determine how the data will need to scaled to be properly combined - typename MinimumMaximumImageCalculator::Pointer inputCalculator = - MinimumMaximumImageCalculator::New(); - typename MinimumMaximumImageCalculator::Pointer filteredCalculator = - MinimumMaximumImageCalculator::New(); + // Combine the input image and the laplacian image - inputCalculator->SetImage(this->GetInput()); - inputCalculator->SetRegion(this->GetOutput()->GetRequestedRegion()); - inputCalculator->Compute(); + using AddImageFilterType = AddImageFilter; + using MultiplyImageFilterType = MultiplyImageFilter; + using SubtractImageFilterType = SubtractImageFilter; - filteredCalculator->SetImage(filter->GetOutput()); - filteredCalculator->SetRegion(this->GetOutput()->GetRequestedRegion()); - filteredCalculator->Compute(); + auto addImageFilter1 = AddImageFilterType::New(); + addImageFilter1->SetInput(laplacianFilter->GetOutput()); + addImageFilter1->SetConstant2(-filteredShift); - RealType inputShift, inputScale, filteredShift, filteredScale; - inputShift = static_cast(inputCalculator->GetMinimum()); - inputScale = - static_cast(inputCalculator->GetMaximum()) - static_cast(inputCalculator->GetMinimum()); + auto multiplyImageFilter = MultiplyImageFilterType::New(); + multiplyImageFilter->SetInput(addImageFilter1->GetOutput()); + multiplyImageFilter->SetConstant(inputScale / filteredScale); - filteredShift = filteredCalculator->GetMinimum(); // no need to cast - filteredScale = filteredCalculator->GetMaximum() - filteredCalculator->GetMinimum(); + auto addImageFilter2 = AddImageFilterType::New(); + addImageFilter2->SetInput(multiplyImageFilter->GetOutput()); + addImageFilter2->SetConstant2(inputShift); - ImageRegionIterator it(filter->GetOutput(), filter->GetOutput()->GetRequestedRegion()); - ImageRegionConstIterator inIt(this->GetInput(), this->GetOutput()->GetRequestedRegion()); + auto subtractImageFilter = SubtractImageFilterType::New(); + subtractImageFilter->SetInput1(localInput); + subtractImageFilter->SetInput2(addImageFilter2->GetOutput()); - // combine the input and laplacian images - RealType value, invalue; - RealType inputSum = 0.0; - RealType enhancedSum = 0.0; - while (!it.IsAtEnd()) - { - value = it.Get(); // laplacian value + // Calculate needed combined image statistics - // rescale to [0,1] - value = (value - filteredShift) / filteredScale; + typename StatisticsImageFilter::Pointer enhancedCalculator = + StatisticsImageFilter::New(); - // rescale to the input dynamic range - value = value * inputScale + inputShift; + enhancedCalculator->SetInput(subtractImageFilter->GetOutput()); + enhancedCalculator->Update(); - // combine the input and laplacian image (note that we subtract - // the laplacian due to the signs in our laplacian kernel). - invalue = static_cast(inIt.Get()); - value = invalue - value; - it.Set(value); + auto enhancedMean = static_cast(enhancedCalculator->GetMean()); - inputSum += invalue; - enhancedSum += value; - ++it; - ++inIt; - } - RealType inputMean = inputSum / static_cast(this->GetOutput()->GetRequestedRegion().GetNumberOfPixels()); - RealType enhancedMean = - enhancedSum / static_cast(this->GetOutput()->GetRequestedRegion().GetNumberOfPixels()); - - // update progress - this->UpdateProgress(0.9); - - // copy and cast the output - typename TOutputImage::Pointer output = this->GetOutput(); - output->SetBufferedRegion(output->GetRequestedRegion()); - output->Allocate(); - - RealType inputMinimum = inputCalculator->GetMinimum(); - RealType inputMaximum = inputCalculator->GetMaximum(); - auto castInputMinimum = static_cast(inputMinimum); - auto castInputMaximum = static_cast(inputMaximum); - - ImageRegionIterator outIt(output, output->GetRequestedRegion()); - it.GoToBegin(); - while (!outIt.IsAtEnd()) - { - value = it.Get(); - - // adjust value to make the mean intensities before and after match - value = value - enhancedMean + inputMean; - - if (value < inputMinimum) - { - outIt.Set(castInputMinimum); - } - else if (value > inputMaximum) - { - outIt.Set(castInputMaximum); - } - else - { - outIt.Set(static_cast(value)); - } - - ++outIt; - ++it; - } + // Shift and window the output + + using IntensityWindowingFilterType = IntensityWindowingImageFilter; + auto intensityWindowingFilter = IntensityWindowingFilterType::New(); + intensityWindowingFilter->SetInput(subtractImageFilter->GetOutput()); + intensityWindowingFilter->SetOutputMinimum(inputMinimum); + intensityWindowingFilter->SetOutputMaximum(inputMaximum); + intensityWindowingFilter->SetWindowMinimum(inputMinimum - (inputMean - enhancedMean)); + intensityWindowingFilter->SetWindowMaximum(inputMaximum - (inputMean - enhancedMean)); + intensityWindowingFilter->GraftOutput(this->GetOutput()); + intensityWindowingFilter->Update(); - // update progress - this->UpdateProgress(1.0); + this->GraftOutput(intensityWindowingFilter->GetOutput()); } + +template +void +LaplacianSharpeningImageFilter::PrintSelf(std::ostream & os, Indent indent) const +{ + Superclass::PrintSelf(os, indent); + + os << indent << "UseImageSpacing: " << (m_UseImageSpacing ? "On" : "Off") << std::endl; +} + } // end namespace itk #endif diff --git a/Modules/Filtering/ImageFeature/test/itkLaplacianSharpeningImageFilterTest.cxx b/Modules/Filtering/ImageFeature/test/itkLaplacianSharpeningImageFilterTest.cxx index 404a629547f..8bb6d8dd5fa 100644 --- a/Modules/Filtering/ImageFeature/test/itkLaplacianSharpeningImageFilterTest.cxx +++ b/Modules/Filtering/ImageFeature/test/itkLaplacianSharpeningImageFilterTest.cxx @@ -17,6 +17,7 @@ *=========================================================================*/ #include "itkLaplacianSharpeningImageFilter.h" +#include "itkCastImageFilter.h" #include "itkImageFileReader.h" #include "itkImageFileWriter.h" #include "itkTestingMacros.h" @@ -34,7 +35,7 @@ itkLaplacianSharpeningImageFilterTest(int argc, char * argv[]) return EXIT_FAILURE; } - using PixelType = unsigned char; + using PixelType = float; constexpr unsigned int Dimension = 2; using ImageType = itk::Image; @@ -62,8 +63,16 @@ itkLaplacianSharpeningImageFilterTest(int argc, char * argv[]) ITK_TRY_EXPECT_NO_EXCEPTION(laplacianSharpeningFilter->Update()); - itk::WriteImage(laplacianSharpeningFilter->GetOutput(), argv[2]); + using OutputPixelType = unsigned char; + using OutputImageType = itk::Image; + using CasterType = itk::CastImageFilter; + + auto caster = CasterType::New(); + caster->SetInput(laplacianSharpeningFilter->GetOutput()); + caster->Update(); + + itk::WriteImage(caster->GetOutput(), argv[2]); std::cout << "Test finished." << std::endl; return EXIT_SUCCESS; diff --git a/Modules/Filtering/ImageFeature/wrapping/itkLaplacianSharpeningImageFilter.wrap b/Modules/Filtering/ImageFeature/wrapping/itkLaplacianSharpeningImageFilter.wrap index 375775f503c..68f19ffff3f 100644 --- a/Modules/Filtering/ImageFeature/wrapping/itkLaplacianSharpeningImageFilter.wrap +++ b/Modules/Filtering/ImageFeature/wrapping/itkLaplacianSharpeningImageFilter.wrap @@ -1,3 +1,3 @@ itk_wrap_class("itk::LaplacianSharpeningImageFilter" POINTER) - itk_wrap_image_filter("${WRAP_ITK_SCALAR}" 2) + itk_wrap_image_filter("${WRAP_ITK_REAL}" 2) itk_end_wrap_class()