Skip to content

Commit

Permalink
ENH: refactor LaplacianSharpeningImageFilter (#4052)
Browse files Browse the repository at this point in the history
  • Loading branch information
@ntustison
ntustison committed Jun 14, 2023
1 parent d3b5e8b commit c93f068
Show file tree
Hide file tree
Showing 5 changed files with 114 additions and 188 deletions.
8 changes: 5 additions & 3 deletions Examples/Filtering/LaplacianSharpeningImageFilter.cxx
View file
Expand Up @@ -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<PixelType, Dimension>;
using CharImageType = itk::Image<CharPixelType, Dimension>;

using ReaderType = itk::ImageFileReader<CharImageType>;
using ReaderType = itk::ImageFileReader<ImageType>;
using WriterType = itk::ImageFileWriter<CharImageType>;

using RescaleFilter =
itk::RescaleIntensityImageFilter<CharImageType, CharImageType>;
itk::RescaleIntensityImageFilter<ImageType, CharImageType>;

using LaplacianSharpeningFilter =
itk::LaplacianSharpeningImageFilter<CharImageType, CharImageType>;
itk::LaplacianSharpeningImageFilter<ImageType, ImageType>;


// Setting the IO
Expand Down
24 changes: 13 additions & 11 deletions Modules/Filtering/ImageFeature/include/itkLaplacianSharpeningImageFilter.h
View file
Expand Up @@ -66,6 +66,7 @@ class ITK_TEMPLATE_EXPORT LaplacianSharpeningImageFilter : public ImageToImageFi
using RealType = typename NumericTraits<OutputPixelType>::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 */
Expand All @@ -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. */
Expand All @@ -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<InputImageDimension, ImageDimension>));
itkConceptMacro(InputPixelTypeIsFloatingPointCheck, (Concept::IsFloatingPoint<InputPixelType>));
itkConceptMacro(OutputPixelTypeIsFloatingPointCheck, (Concept::IsFloatingPoint<OutputPixelType>));
// 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
Expand Down
255 changes: 84 additions & 171 deletions Modules/Filtering/ImageFeature/include/itkLaplacianSharpeningImageFilter.hxx
View file
Expand Up @@ -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 <typename TInputImage, typename TOutputImage>
void
LaplacianSharpeningImageFilter<TInputImage, TOutputImage>::PrintSelf(std::ostream & os, Indent indent) const
LaplacianSharpeningImageFilter<TInputImage, TOutputImage>::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 <typename TInputImage, typename TOutputImage>
void
LaplacianSharpeningImageFilter<TInputImage, TOutputImage>::GenerateInputRequestedRegion()
LaplacianSharpeningImageFilter<TInputImage, TOutputImage>::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<TInputImage *>(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<RealType, ImageDimension> oper;
oper.CreateOperator();
typename StatisticsImageFilter<InputImageType>::Pointer inputCalculator =
StatisticsImageFilter<InputImageType>::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<RealType>(inputCalculator->GetMinimum());
auto inputMaximum = static_cast<RealType>(inputCalculator->GetMaximum());
auto inputMean = static_cast<RealType>(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 <typename TInputImage, typename TOutputImage>
void
LaplacianSharpeningImageFilter<TInputImage, TOutputImage>::GenerateData()
{
// Create the Laplacian operator
LaplacianOperator<RealType, ImageDimension> 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<RealType, ImageDimension>;
using NOIF = NeighborhoodOperatorImageFilter<InputImageType, RealImageType>;
ZeroFluxNeumannBoundaryCondition<InputImageType> nbc;

auto filter = NOIF::New();
filter->OverrideBoundaryCondition(static_cast<typename NOIF::ImageBoundaryConditionPointerType>(&nbc));
using LaplacianImageFilter = LaplacianImageFilter<InputImageType, RealImageType>;
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<RealImageType>::Pointer filteredCalculator =
MinimumMaximumImageCalculator<RealImageType>::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<RealType>(filteredCalculator->GetMinimum());
RealType filteredScale = static_cast<RealType>(filteredCalculator->GetMaximum()) - filteredShift;

// determine how the data will need to scaled to be properly combined
typename MinimumMaximumImageCalculator<InputImageType>::Pointer inputCalculator =
MinimumMaximumImageCalculator<InputImageType>::New();
typename MinimumMaximumImageCalculator<RealImageType>::Pointer filteredCalculator =
MinimumMaximumImageCalculator<RealImageType>::New();
// Combine the input image and the laplacian image

inputCalculator->SetImage(this->GetInput());
inputCalculator->SetRegion(this->GetOutput()->GetRequestedRegion());
inputCalculator->Compute();
using AddImageFilterType = AddImageFilter<RealImageType>;
using MultiplyImageFilterType = MultiplyImageFilter<RealImageType>;
using SubtractImageFilterType = SubtractImageFilter<InputImageType, RealImageType, RealImageType>;

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<RealType>(inputCalculator->GetMinimum());
inputScale =
static_cast<RealType>(inputCalculator->GetMaximum()) - static_cast<RealType>(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<RealImageType> it(filter->GetOutput(), filter->GetOutput()->GetRequestedRegion());
ImageRegionConstIterator<InputImageType> 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<RealImageType>::Pointer enhancedCalculator =
StatisticsImageFilter<RealImageType>::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<RealType>(inIt.Get());
value = invalue - value;
it.Set(value);
auto enhancedMean = static_cast<RealType>(enhancedCalculator->GetMean());

inputSum += invalue;
enhancedSum += value;
++it;
++inIt;
}
RealType inputMean = inputSum / static_cast<RealType>(this->GetOutput()->GetRequestedRegion().GetNumberOfPixels());
RealType enhancedMean =
enhancedSum / static_cast<RealType>(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<OutputPixelType>(inputMinimum);
auto castInputMaximum = static_cast<OutputPixelType>(inputMaximum);

ImageRegionIterator<OutputImageType> 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<OutputPixelType>(value));
}

++outIt;
++it;
}
// Shift and window the output

using IntensityWindowingFilterType = IntensityWindowingImageFilter<RealImageType, OutputImageType>;
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 <typename TInputImage, typename TOutputImage>
void
LaplacianSharpeningImageFilter<TInputImage, TOutputImage>::PrintSelf(std::ostream & os, Indent indent) const
{
Superclass::PrintSelf(os, indent);

os << indent << "UseImageSpacing: " << (m_UseImageSpacing ? "On" : "Off") << std::endl;
}

} // end namespace itk

#endif

0 comments on commit c93f068

Please sign in to comment.