ITK  4.8.0 Insight Segmentation and Registration Toolkit
Examples/Filtering/BilateralImageFilter.cxx
/*=========================================================================
*
*
* you may not use this file except in compliance with the License.
* You may obtain a copy of the License at
*
*
* Unless required by applicable law or agreed to in writing, software
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
*
*=========================================================================*/
// 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
//
//
// 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}
//
//
// 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
//
//
// k(\mathbf{x}) = \int_\omega c(\mathbf{x},\mathbf{w})
// s( f(\mathbf{x}),f(\mathbf{w})) d \mathbf{w}
//
//
// The default implementation of this filter uses Gaussian kernels for both
// $c()$ and $s()$. The $c$ kernel can be described as
//
//
// c(\mathbf{x},\mathbf{w}) = e^{(\frac{ {\left|| \mathbf{x} - \mathbf{w} \right||}^2 }{\sigma^2_c} )}
//
//
// 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
//
//
// s(f(\mathbf{x}),f(\mathbf{w})) = e^{(\frac{ {( f(\mathbf{x}) - f(\mathbf{w})}^2 }{\sigma^2_s} )}
//
//
// 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"
// Software Guide : BeginLatex
//
// The first step required to use this filter is to include its header file.
//
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
// 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
typedef unsigned char InputPixelType;
typedef unsigned char OutputPixelType;
typedef itk::Image< InputPixelType, 2 > InputImageType;
typedef itk::Image< OutputPixelType, 2 > OutputImageType;
// Software Guide : EndCodeSnippet
// 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
InputImageType, OutputImageType > FilterType;
FilterType::Pointer filter = FilterType::New();
// Software Guide : EndCodeSnippet
// 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
// 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
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
const unsigned int Dimension = InputImageType::ImageDimension;
double domainSigmas[ Dimension ];
for(unsigned int i=0; i<Dimension; i++)
{
domainSigmas[i] = atof( argv[3] );
}
const double rangeSigma = atof( 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 a intensity rescaler
// filter and then to a writer. Invoking \code{Update()} on the writer
// triggers the execution of both filters.
//
// Software Guide : EndLatex
typedef unsigned char WritePixelType;
typedef itk::Image< WritePixelType, 2 > WriteImageType;
OutputImageType, WriteImageType > RescaleFilterType;
RescaleFilterType::Pointer rescaler = RescaleFilterType::New();
rescaler->SetOutputMinimum( 0 );
rescaler->SetOutputMaximum( 255 );
WriterType::Pointer 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}