Smooth RGB Image Using Min Max Curvature Flow#

Synopsis#

Smooth an RGB image using min/max curvature flow.

Results#

../../../../_images/bafkreibbmonevzco6rqjcqzk6zgjxblw4yby2uwzhjgkwbydfwxj3gpl4e — Output In VTK Window#

Code#

C++#

#include "itkImageAdaptor.h"
#include "itkImageRegionIterator.h"
#include "itkNthElementImageAdaptor.h"
#include "itkMinMaxCurvatureFlowImageFilter.h"
#include "itkComposeImageFilter.h"
#include "itkImage.h"
#include "itkImageFileReader.h"
#include "itkRGBPixel.h"

#include <sstream>

#ifdef ENABLE_QUICKVIEW
#  include "QuickView.h"
#endif

int
main(int argc, char * argv[])
{
  // Verify command line arguments
  if (argc < 2)
  {
    std::cerr << "Usage: " << std::endl;
    std::cerr << argv[0] << " InputImageFile [iterations]" << std::endl;
    return EXIT_FAILURE;
  }

  int iterations = 5;
  if (argc > 2)
  {
    std::stringstream ss(argv[2]);
    ss >> iterations;
  }
  std::string inputFileName = argv[1];

  // Setup types
  using ComponentType = float;
  using PixelType = itk::RGBPixel<ComponentType>;
  using RGBImageType = itk::Image<PixelType, 2>;
  using ImageType = itk::Image<ComponentType, 2>;
  using ImageAdaptorType = itk::NthElementImageAdaptor<RGBImageType, unsigned char>;
  using ComposeType = itk::ComposeImageFilter<ImageType, RGBImageType>;
  using CurvatureFlowType = itk::MinMaxCurvatureFlowImageFilter<ImageAdaptorType, ImageType>;

  const auto input = itk::ReadImage<RGBImageType>(inputFileName);

  // Run the filter for each component
  auto rAdaptor = ImageAdaptorType::New();
  rAdaptor->SelectNthElement(0);
  rAdaptor->SetImage(input);

  auto rCurvatureFilter = CurvatureFlowType::New();
  rCurvatureFilter->SetInput(rAdaptor);
  rCurvatureFilter->SetNumberOfIterations(iterations);
  rCurvatureFilter->Update();

  auto gAdaptor = ImageAdaptorType::New();
  gAdaptor->SelectNthElement(1);
  gAdaptor->SetImage(input);

  auto gCurvatureFilter = CurvatureFlowType::New();
  gCurvatureFilter->SetInput(gAdaptor);
  gCurvatureFilter->SetNumberOfIterations(iterations);
  gCurvatureFilter->Update();

  auto bAdaptor = ImageAdaptorType::New();
  bAdaptor->SelectNthElement(2);
  bAdaptor->SetImage(input);

  auto bCurvatureFilter = CurvatureFlowType::New();
  bCurvatureFilter->SetInput(bAdaptor);
  bCurvatureFilter->SetNumberOfIterations(iterations);
  bCurvatureFilter->Update();

  // compose an RGB image from the three filtered images

  auto compose = ComposeType::New();
  compose->SetInput1(rCurvatureFilter->GetOutput());
  compose->SetInput2(gCurvatureFilter->GetOutput());
  compose->SetInput3(bCurvatureFilter->GetOutput());

#ifdef ENABLE_QUICKVIEW
  QuickView viewer;
  viewer.AddRGBImage(input.GetPointer(), true, itksys::SystemTools::GetFilenameName(inputFileName));

  std::stringstream desc;
  desc << "MinMaxCurvatureFlow iterations = " << iterations;
  viewer.AddRGBImage(compose->GetOutput(), true, desc.str());

  viewer.Visualize();
#endif

  return EXIT_SUCCESS;
}

Classes demonstrated#

template<typename TInputImage, typename TOutputImage> class MinMaxCurvatureFlowImageFilter : public itk::CurvatureFlowImageFilter<TInputImage, TOutputImage>

Denoise an image using min/max curvature flow.

MinMaxCurvatureFlowImageFilter implements a curvature driven image denoising algorithm [108]. Iso-brightness contours in the grayscale input image are viewed as a level set. The level set is then evolved using a curvature-based speed function:

$I_t = F_{\mbox{minmax}} |\nabla I|$

where $F_{\mbox{minmax}} = \max(\kappa,0)$ if $\mbox{Avg}_{\mbox{stencil}}(x)$ is less than or equal to $T_{threshold}$ and $\min(\kappa,0)$ , otherwise. $\kappa$ is the mean curvature of the iso-brightness contour at point $x$ .

In min/max curvature flow, movement is turned on or off depending on the scale of the noise one wants to remove. Switching depends on the average image value of a region of radius $R$ around each point. The choice of $R$ , the stencil radius, governs the scale of the noise to be removed.

The threshold value $T_{threshold}$ is the average intensity obtained in the direction perpendicular to the gradient at point $x$ at the extrema of the local neighborhood.

This filter make use of the multi-threaded finite difference solver hierarchy. Updates are computed using a MinMaxCurvatureFlowFunction object. A zero flux Neumann boundary condition is used when computing derivatives near the data boundary.