Examples/Filtering/BilateralImageFilter.cxx

/*=========================================================================
 *
 *  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.
 *
 *=========================================================================*/
 
//  Software Guide : BeginCommandLineArgs
//    INPUTS:  {BrainProtonDensitySlice.png}
//    OUTPUTS: {BilateralImageFilterOutput.png}
//    ARGUMENTS:    6 5
//  Software Guide : EndCommandLineArgs
 
//  Software Guide : BeginLatex
//
//  The \doxygen{BilateralImageFilter} performs smoothing by using both
//  domain and range neighborhoods.  Pixels that are close to a pixel in the
//  image domain and similar to a pixel in the image range are used to
//  calculate the filtered value. Two Gaussian kernels (one in the image
//  domain and one in the image range) are used to smooth the image. The
//  result is an image that is smoothed in homogeneous regions yet has edges
//  preserved. The result is similar to anisotropic diffusion but the
//  implementation is non-iterative.  Another benefit to bilateral filtering
//  is that any distance metric can be used for kernel smoothing the image
//  range. Bilateral filtering is capable of reducing the noise in an image
//  by an order of magnitude while maintaining edges.  The bilateral operator
//  used here was described by Tomasi and Manduchi (\emph{Bilateral Filtering
//  for Gray and Color Images}. IEEE ICCV. 1998.)
//
//  The filtering operation can be described by the following equation
//
//  \begin{equation}
//  h(\mathbf{x}) = k(\mathbf{x})^{-1} \int_\omega f(\mathbf{w})
//  c(\mathbf{x},\mathbf{w}) s( f(\mathbf{x}),f(\mathbf{w})) d \mathbf{w}
//  \end{equation}
//
//  where $\mathbf{x}$ holds the coordinates of a $ND$ point, $f(\mathbf{x})$
//  is the input image and $h(\mathbf{x})$ is the output image. The
//  convolution kernels $c()$ and $s()$ are associated with the spatial and
//  intensity domain respectively. The $ND$ integral is computed over
//  $\omega$ which is a neighborhood of the pixel located at
//  $\mathbf{x}$. The normalization factor $k(\mathbf{x})$ is computed as
//
//  \begin{equation}
//  k(\mathbf{x}) = \int_\omega c(\mathbf{x},\mathbf{w})
//  s( f(\mathbf{x}),f(\mathbf{w})) d \mathbf{w}
//  \end{equation}
//
//  The default implementation of this filter uses Gaussian kernels for both
//  $c()$ and $s()$. The $c$ kernel can be described as
//
//  \begin{equation}
//  c(\mathbf{x},\mathbf{w}) = e^{(- \frac{ {\left|| \mathbf{x} - \mathbf{w}
//  \right||}^2
//  }{\sigma^2_c} )} \end{equation}
//
//  where $\sigma_c$ is provided by the user and defines how close pixel
//  neighbors should be in order to be considered for the computation of the
//  output value.  The $s$ kernel is given by
//
//  \begin{equation}
//  s(f(\mathbf{x}),f(\mathbf{w})) = e^{(- \frac{ {( f(\mathbf{x}) -
//  f(\mathbf{w}) )}^2
//  }{\sigma^2_s} )} \end{equation}
//
//  where $\sigma_s$ is provided by the user and defines how close the
//  neighbor's intensity be in order to be considered for the computation of
//  the output value.
//
//  \index{itk::BilateralImageFilter}
//
//  Software Guide : EndLatex
 
 
#include "itkImage.h"
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
#include "itkRescaleIntensityImageFilter.h"
 
//  Software Guide : BeginLatex
//
//  The first step required to use this filter is to include its header file.
//
//  \index{itk::BilateralImageFilter!header}
//
//  Software Guide : EndLatex
 
// Software Guide : BeginCodeSnippet
#include "itkBilateralImageFilter.h"
// Software Guide : EndCodeSnippet
 
 
int
main(int argc, char * argv[])
{
  if (argc < 5)
  {
    std::cerr << "Usage: " << std::endl;
    std::cerr << argv[0]
              << "  inputImageFile  outputImageFile  domainSigma  rangeSigma"
              << std::endl;
    return EXIT_FAILURE;
  }
 
 
  //  Software Guide : BeginLatex
  //
  //  The image types are instantiated using pixel type and dimension.
  //
  //  Software Guide : EndLatex
 
  // Software Guide : BeginCodeSnippet
  using InputPixelType = unsigned char;
  using OutputPixelType = unsigned char;
 
  using InputImageType = itk::Image<InputPixelType, 2>;
  using OutputImageType = itk::Image<OutputPixelType, 2>;
  // Software Guide : EndCodeSnippet
 
 
  using ReaderType = itk::ImageFileReader<InputImageType>;
 
 
  //  Software Guide : BeginLatex
  //
  //  The bilateral filter type is now instantiated using both the input
  //  image and the output image types and the filter object is created.
  //
  //  \index{itk::BilateralImageFilter!instantiation}
  //  \index{itk::BilateralImageFilter!New()}
  //  \index{itk::BilateralImageFilter!Pointer}
  //
  //  Software Guide : EndLatex
 
  // Software Guide : BeginCodeSnippet
  using FilterType =
    itk::BilateralImageFilter<InputImageType, OutputImageType>;
  auto filter = FilterType::New();
  // Software Guide : EndCodeSnippet
 
 
  auto reader = ReaderType::New();
  reader->SetFileName(argv[1]);
 
 
  //  Software Guide : BeginLatex
  //
  //  The input image can be obtained from the output of another
  //  filter. Here, an image reader is used as a source.
  //
  //  Software Guide : EndLatex
 
  // Software Guide : BeginCodeSnippet
  filter->SetInput(reader->GetOutput());
  // Software Guide : EndCodeSnippet
 
 
  //  Software Guide : BeginLatex
  //
  //  The Bilateral filter requires two parameters. First, we must specify the
  //  standard deviation $\sigma$ to be used for the Gaussian kernel on image
  //  intensities. Second, the set of $\sigma$s to be used along each
  //  dimension in the space domain. This second parameter is supplied as an
  //  array of \code{float} or \code{double} values. The array dimension
  //  matches the image dimension. This mechanism makes it possible to enforce
  //  more coherence along some directions. For example, more smoothing can be
  //  done along the $X$ direction than along the $Y$ direction.
  //
  //  In the following code example, the $\sigma$ values are taken from the
  //  command line.  Note the use of \code{ImageType::ImageDimension} to get
  //  access to the image dimension at compile time.
  //
  //  Software Guide : EndLatex
 
  // Software Guide : BeginCodeSnippet
  constexpr unsigned int Dimension = InputImageType::ImageDimension;
  double                 domainSigmas[Dimension];
  for (double & domainSigma : domainSigmas)
  {
    domainSigma = std::stod(argv[3]);
  }
  const double rangeSigma = std::stod(argv[4]);
  // Software Guide : EndCodeSnippet
 
 
  //  Software Guide : BeginLatex
  //
  //  The filter parameters are set with the methods \code{SetRangeSigma()}
  //  and \code{SetDomainSigma()}.
  //
  //  \index{itk::BilateralImageFilter!SetRangeSigma()}
  //  \index{itk::BilateralImageFilter!SetDomainSigma()}
  //  \index{SetDomainSigma()!itk::BilateralImageFilter}
  //  \index{SetRangeSigma()!itk::BilateralImageFilter}
  //
  //  Software Guide : EndLatex
 
 
  // Software Guide : BeginCodeSnippet
  filter->SetDomainSigma(domainSigmas);
  filter->SetRangeSigma(rangeSigma);
  // Software Guide : EndCodeSnippet
 
 
  //  Software Guide : BeginLatex
  //
  //  The output of the filter is connected here to an intensity rescaler
  //  filter and then to a writer. Invoking \code{Update()} on the writer
  //  triggers the execution of both filters.
  //
  //  Software Guide : EndLatex
 
  using WritePixelType = unsigned char;
  using WriteImageType = itk::Image<WritePixelType, 2>;
  using RescaleFilterType =
    itk::RescaleIntensityImageFilter<OutputImageType, WriteImageType>;
  auto rescaler = RescaleFilterType::New();
  rescaler->SetOutputMinimum(0);
  rescaler->SetOutputMaximum(255);
 
  using WriterType = itk::ImageFileWriter<WriteImageType>;
  auto writer = WriterType::New();
  writer->SetFileName(argv[2]);
 
  // Software Guide : BeginCodeSnippet
  rescaler->SetInput(filter->GetOutput());
  writer->SetInput(rescaler->GetOutput());
  writer->Update();
  // Software Guide : EndCodeSnippet
 
 
  //  Software Guide : BeginLatex
  //
  // \begin{figure}
  // \center
  // \includegraphics[width=0.44\textwidth]{BrainProtonDensitySlice}
  // \includegraphics[width=0.44\textwidth]{BilateralImageFilterOutput}
  // \itkcaption[BilateralImageFilter output]{Effect of the
  // BilateralImageFilter on a slice from a MRI proton density image  of the
  // brain.} \label{fig:BilateralImageFilterInputOutput} \end{figure}
  //
  //  Figure \ref{fig:BilateralImageFilterInputOutput} illustrates the effect
  //  of this filter on a MRI proton density image of the brain. In this
  //  example the filter was run with a range $\sigma$ of $5.0$ and a domain
  //  $\sigma$ of $6.0$.  The figure shows how homogeneous regions are
  //  smoothed and edges are preserved.
  //
  //  \relatedClasses
  //  \begin{itemize}
  //  \item \doxygen{GradientAnisotropicDiffusionImageFilter}
  //  \item \doxygen{CurvatureAnisotropicDiffusionImageFilter}
  //  \item \doxygen{CurvatureFlowImageFilter}
  //  \end{itemize}
  //
  //  Software Guide : EndLatex
 
  return EXIT_SUCCESS;
}