Filtering » Linear filters module

Linear smoothing, sharpening and derivative filters.

Contents

Classes

struct dip::OneDimensionalFilter
Describes a 1D filter

Aliases

using dip::OneDimensionalFilterArray = std::vector<OneDimensionalFilter>
An array of 1D filters

Functions

auto dip::SeparateFilter(dip::Image const& filter) -> dip::OneDimensionalFilterArray
Separates a linear filter (convolution kernel) into a set of 1D filters that can be applied using dip::SeparableConvolution.
void dip::SeparableConvolution(dip::Image const& in, dip::Image& out, dip::OneDimensionalFilterArray const& filterArray, dip::StringArray const& boundaryCondition = {}, dip::BooleanArray process = {})
Applies a convolution with a filter kernel (PSF) that is separable.
void dip::ConvolveFT(dip::Image const& in, dip::Image const& filter, dip::Image& out, dip::String const& inRepresentation = S::SPATIAL, dip::String const& filterRepresentation = S::SPATIAL, dip::String const& outRepresentation = S::SPATIAL, dip::StringArray const& boundaryCondition = {})
Applies a convolution with a filter kernel (PSF) by multiplication in the Fourier domain.
void dip::GeneralConvolution(dip::Image const& in, dip::Image const& filter, dip::Image& out, dip::StringArray const& boundaryCondition = {})
Applies a convolution with a filter kernel (PSF) by direct implementation of the convolution sum.
void dip::Convolution(dip::Image const& in, dip::Image const& filter, dip::Image& out, dip::String const& method = S::BEST, dip::StringArray const& boundaryCondition = {})
Applies a convolution with a filter kernel (PSF).
void dip::Uniform(dip::Image const& in, dip::Image& out, dip::Kernel const& kernel = {}, dip::StringArray const& boundaryCondition = {})
Applies a convolution with a kernel with uniform weights, leading to an average (mean) filter.
void dip::GaussFIR(dip::Image const& in, dip::Image& out, dip::FloatArray sigmas = {1.0}, dip::UnsignedArray derivativeOrder = {0}, dip::StringArray const& boundaryCondition = {}, dip::dfloat truncation = 3)
Finite impulse response implementation of the Gaussian filter and its derivatives
void dip::GaussFT(dip::Image const& in, dip::Image& out, dip::FloatArray sigmas = {1.0}, dip::UnsignedArray derivativeOrder = {0}, dip::dfloat truncation = 3, dip::String const& inRepresentation = S::SPATIAL, dip::String const& outRepresentation = S::SPATIAL, dip::StringArray const& boundaryCondition = {})
Fourier implementation of the Gaussian filter and its derivatives
void dip::GaussIIR(dip::Image const& in, dip::Image& out, dip::FloatArray sigmas = {1.0}, dip::UnsignedArray derivativeOrder = {0}, dip::StringArray const& boundaryCondition = {}, dip::UnsignedArray filterOrder = {}, dip::String const& designMethod = S::DISCRETE_TIME_FIT, dip::dfloat truncation = 3)
Infinite impulse response implementation of the Gaussian filter and its derivatives
void dip::Gauss(dip::Image const& in, dip::Image& out, dip::FloatArray sigmas = {1.0}, dip::UnsignedArray derivativeOrder = {0}, dip::String const& method = S::BEST, dip::StringArray const& boundaryCondition = {}, dip::dfloat truncation = 3)
Convolution with a Gaussian kernel and its derivatives
void dip::FiniteDifference(dip::Image const& in, dip::Image& out, dip::UnsignedArray derivativeOrder = {0}, dip::String const& smoothFlag = S::SMOOTH, dip::StringArray const& boundaryCondition = {}, dip::BooleanArray process = {})
Finite difference derivatives
void dip::SobelGradient(dip::Image const& in, dip::Image& out, dip::uint dimension = 0, dip::StringArray const& boundaryCondition = {})
The Sobel derivative filter
void dip::Derivative(dip::Image const& in, dip::Image& out, dip::UnsignedArray derivativeOrder, dip::FloatArray sigmas = {1.0}, dip::String const& method = S::BEST, dip::StringArray const& boundaryCondition = {}, dip::dfloat truncation = 3)
Computes derivatives
void dip::Dx(dip::Image const& in, dip::Image& out, dip::FloatArray sigmas = {1.0})
Computes the first derivative along x, see dip::Derivative.
void dip::Dy(dip::Image const& in, dip::Image& out, dip::FloatArray sigmas = {1.0})
Computes the first derivative along y, see dip::Derivative.
void dip::Dz(dip::Image const& in, dip::Image& out, dip::FloatArray sigmas = {1.0})
Computes the first derivative along z, see dip::Derivative.
void dip::Dxx(dip::Image const& in, dip::Image& out, dip::FloatArray sigmas = {1.0})
Computes the second derivative along x, see dip::Derivative.
void dip::Dyy(dip::Image const& in, dip::Image& out, dip::FloatArray sigmas = {1.0})
Computes the second derivative along y, see dip::Derivative.
void dip::Dzz(dip::Image const& in, dip::Image& out, dip::FloatArray sigmas = {1.0})
Computes the second derivative along z, see dip::Derivative.
void dip::Dxy(dip::Image const& in, dip::Image& out, dip::FloatArray sigmas = {1.0})
Computes the first derivative along x and y, see dip::Derivative.
void dip::Dxz(dip::Image const& in, dip::Image& out, dip::FloatArray sigmas = {1.0})
Computes the first derivative along x and z, see dip::Derivative.
void dip::Dyz(dip::Image const& in, dip::Image& out, dip::FloatArray sigmas = {1.0})
Computes the first derivative along y and y, see dip::Derivative.
void dip::Gradient(dip::Image const& in, dip::Image& out, dip::FloatArray sigmas = {1.0}, dip::String const& method = S::BEST, dip::StringArray const& boundaryCondition = {}, dip::BooleanArray process = {}, dip::dfloat truncation = 3)
Computes the gradient of the image, resulting in an N-vector image, if the input was N-dimensional.
void dip::GradientMagnitude(dip::Image const& in, dip::Image& out, dip::FloatArray sigmas = {1.0}, dip::String const& method = S::BEST, dip::StringArray const& boundaryCondition = {}, dip::BooleanArray process = {}, dip::dfloat truncation = 3)
Computes the gradient magnitude of the image, equivalent to dip::Norm( dip::Gradient( in )).
void dip::GradientDirection(dip::Image const& in, dip::Image& out, dip::FloatArray sigmas = {1.0}, dip::String const& method = S::BEST, dip::StringArray const& boundaryCondition = {}, dip::BooleanArray process = {}, dip::dfloat truncation = 3)
Computes the direction of the gradient of the image, equivalent to dip::Angle( dip::Gradient( in )).
void dip::Curl(dip::Image const& in, dip::Image& out, dip::FloatArray sigmas = {1.0}, dip::String const& method = S::BEST, dip::StringArray const& boundaryCondition = {}, dip::BooleanArray process = {}, dip::dfloat truncation = 3)
Computes the curl (rotation) of the 2D or 3D vector field in.
void dip::Divergence(dip::Image const& in, dip::Image& out, dip::FloatArray sigmas = {1.0}, dip::String const& method = S::BEST, dip::StringArray const& boundaryCondition = {}, dip::BooleanArray process = {}, dip::dfloat truncation = 3)
Computes the divergence of the vector field in.
void dip::Hessian(dip::Image const& in, dip::Image& out, dip::FloatArray sigmas = {1.0}, dip::String const& method = S::BEST, dip::StringArray const& boundaryCondition = {}, dip::BooleanArray process = {}, dip::dfloat truncation = 3)
Computes the Hessian of the image, resulting in a symmetric NxN tensor image, if the input was N-dimensional.
void dip::Laplace(dip::Image const& in, dip::Image& out, dip::FloatArray sigmas = {1.0}, dip::String const& method = S::BEST, dip::StringArray const& boundaryCondition = {}, dip::BooleanArray process = {}, dip::dfloat truncation = 3)
Computes the Laplacian of the image, equivalent to dip::Trace( dip::Hessian( in )), but more efficient.
void dip::Dgg(dip::Image const& in, dip::Image& out, dip::FloatArray sigmas = {1.0}, dip::String const& method = S::BEST, dip::StringArray const& boundaryCondition = {}, dip::BooleanArray process = {}, dip::dfloat truncation = 3)
Computes the second derivative in the gradient direction.
void dip::LaplacePlusDgg(dip::Image const& in, dip::Image& out, dip::FloatArray sigmas = {1.0}, dip::String const& method = S::BEST, dip::StringArray const& boundaryCondition = {}, dip::BooleanArray process = {}, dip::dfloat truncation = 3)
Adds the second derivative in the gradient direction to the Laplacian.
void dip::LaplaceMinusDgg(dip::Image const& in, dip::Image& out, dip::FloatArray sigmas = {1.0}, dip::String const& method = S::BEST, dip::StringArray const& boundaryCondition = {}, dip::BooleanArray process = {}, dip::dfloat truncation = 3)
Subtracts the second derivative in the gradient direction from the Laplacian.
void dip::Sharpen(dip::Image const& in, dip::Image& out, dip::dfloat weight = 1.0, dip::FloatArray sigmas = {1.0}, dip::String const& method = S::BEST, dip::StringArray const& boundaryCondition = {}, dip::dfloat truncation = 3)
Sharpens in by subtracting the Laplacian of the image.
void dip::UnsharpMask(dip::Image const& in, dip::Image& out, dip::dfloat weight = 1.0, dip::FloatArray sigmas = {1.0}, dip::String const& method = S::BEST, dip::StringArray const& boundaryCondition = {}, dip::dfloat truncation = 3)
Sharpens in by subtracting the smoothed image.
void dip::GaborFIR(dip::Image const& in, dip::Image& out, dip::FloatArray sigmas, dip::FloatArray const& frequencies, dip::StringArray const& boundaryCondition = {}, dip::BooleanArray process = {}, dip::dfloat truncation = 3)
Finite impulse response implementation of the Gabor filter
void dip::GaborIIR(dip::Image const& in, dip::Image& out, dip::FloatArray sigmas, dip::FloatArray const& frequencies, dip::StringArray const& boundaryCondition = {}, dip::BooleanArray process = {}, dip::IntegerArray const& filterOrder = {}, dip::dfloat truncation = 3)
Recursive infinite impulse response implementation of the Gabor filter
void dip::Gabor2D(dip::Image const& in, dip::Image& out, dip::FloatArray sigmas = {5.0,5.0}, dip::dfloat frequency = 0.1, dip::dfloat direction = dip::pi, dip::StringArray const& boundaryCondition = {}, dip::dfloat truncation = 3)
2D Gabor filter with direction parameter
void dip::LogGaborFilterBank(dip::Image const& in, dip::Image& out, dip::FloatArray const& wavelengths = {3.0,6.0,12.0,24.0}, dip::dfloat bandwidth = 0.75, dip::uint nOrientations = 6, dip::String const& inRepresentation = S::SPATIAL, dip::String const& outRepresentation = S::SPATIAL)
Applies a log-Gabor filter bank
void dip::NormalizedConvolution(dip::Image const& in, dip::Image const& mask, dip::Image& out, dip::FloatArray const& sigmas = {1.0}, dip::String const& method = S::BEST, dip::StringArray const& boundaryCondition = {S::ADD_ZEROS}, dip::dfloat truncation = 3)
Computes the normalized convolution with a Gaussian kernel: a Gaussian convolution for missing or uncertain data.
void dip::NormalizedDifferentialConvolution(dip::Image const& in, dip::Image const& mask, dip::Image& out, dip::uint dimension = 0, dip::FloatArray const& sigmas = {1.0}, dip::String const& method = S::BEST, dip::StringArray const& boundaryCondition = {S::ADD_ZEROS}, dip::dfloat truncation = 3)
Computes the normalized differential convolution with a Gaussian kernel: a derivative operator for missing or uncertain data.
void dip::MeanShiftVector(dip::Image const& in, dip::Image& out, dip::FloatArray sigmas = {1.0}, dip::String const& method = S::BEST, dip::StringArray const& boundaryCondition = {}, dip::dfloat truncation = 3)
Computes the mean shift vector for each pixel in the image

Class documentation

struct dip::OneDimensionalFilter

Describes a 1D filter

The weights are in filter. If isComplex, then the values in filter are interpreted as real/imaginary pairs. In this case, filter must have an even length, with each two consecutive elements representing a single filter weight. The filter.data() pointer can thus be cast to dip::dcomplex.

The origin is placed either at the index given by origin if it’s non-negative, or at index filter.size() / 2 if origin is negative. Note that filter.size() / 2 is either the middle pixel if the filter is odd in length, or the pixel to the right of the center if it is even in length:

size of filter origin origin location
any 1 x 0 x x x x
any 5 x x x x x 0
any odd value -1 x x 0 x x
any even value -1 x x x 0 x x

Note that, if positive, origin must be an index to one of the samples in the filter array: origin < filter.size().

symmetry indicates the filter shape: "general" (or an empty string) indicates no symmetry. "even" indicates even symmetry, "odd" indicates odd symmetry, and "conj" indicates complex conjugate symmetry. In these three cases, the filter represents the left half of the full filter, with the rightmost element at the origin (and not repeated). The full filter is thus always odd in size. "d-even", "d-odd" and "d-conj" are similar, but duplicate the rightmost element, yielding an even-sized filter. The origin for the symmetric filters is handled identically to the general filter case.

The following table summarizes the result of using various symmetry values. The filter array in all cases has n elements represented in this example as [a,b,c].

symmetry resulting array resulting array length
"general" [a,b,c] n
"even" [a,b,c,b,a] 2n - 1
"odd" [a,b,c,-b,-a] 2n - 1
"conj" [a,b,c,b*,a*] 2n - 1
"d-even" [a,b,c,c,b,a] 2n
"d-odd" [a,b,c,-c,-b,-a] 2n
"d-conj" [a,b,c,c*,b*,a*] 2n

The convolution is applied to each tensor component separately, which is always the correct behavior for linear filters.

Variables
std::vector<dfloat> filter Filter weights.
dip::sint origin Origin of the filter if non-negative.
dip::String symmetry Filter shape: "" == "general", "even", "odd", "conj", "d-even", "d-odd" or "d-conj".
bool isComplex If true, filter contains complex data.

Function documentation

dip::OneDimensionalFilterArray dip::SeparateFilter(dip::Image const& filter)

Separates a linear filter (convolution kernel) into a set of 1D filters that can be applied using dip::SeparableConvolution.

If filter does not represent a separable kernel, the output dip::OneDimensionalFilterArray object is empty (it’s empty method returns true, and it’s size method return 0).

void dip::SeparableConvolution(dip::Image const& in, dip::Image& out, dip::OneDimensionalFilterArray const& filterArray, dip::StringArray const& boundaryCondition = {}, dip::BooleanArray process = {})

Applies a convolution with a filter kernel (PSF) that is separable.

filter is an array with exactly one dip::OneDimensionalFilter element for each dimension of in. Alternatively, it can have a single element, which will be used unchanged for each dimension. For the dimensions that are not processed (process is false for those dimensions), the filter array can have nonsensical data or a zero-length filter weights array. Any filter array that is zero size or the equivalent of {1} will not be applied either.

boundaryCondition indicates how the boundary should be expanded in each dimension. See dip::BoundaryCondition.

process indices which dimensions to process, and can be {} to indicate all dimensions are to be processed.

void dip::ConvolveFT(dip::Image const& in, dip::Image const& filter, dip::Image& out, dip::String const& inRepresentation = S::SPATIAL, dip::String const& filterRepresentation = S::SPATIAL, dip::String const& outRepresentation = S::SPATIAL, dip::StringArray const& boundaryCondition = {})

Applies a convolution with a filter kernel (PSF) by multiplication in the Fourier domain.

filter is an image, and must be equal in size or smaller than in. If both in and filter are real, out will be real too, otherwise it will have a complex type.

As elsewhere, the origin of filter is in the middle of the image, on the pixel to the right of the center in case of an even-sized image.

If in or filter is already Fourier transformed, set inRepresentation or filterRepresentation to "frequency". Similarly, if outRepresentation is "frequency", the output will not be inverse-transformed, so will be in the frequency domain. These three values are "spatial" by default. If any of these three values is "frequency", then out will be complex, no checks are made to see if the inputs in frequency domain have the complex conjugate symmetry required for the result to be real-valued. Use dip::Image::Real if you expect the output to be real-valued in this case.

boundaryCondition indicates how the boundary should be expanded in each dimension. See dip::BoundaryCondition for a description of each option. It is ignored unless inRepresentation, filterRepresentation and outRepresentation are all "spatial". If the array is empty (the default), then a periodic boundary condition is imposed. This is the natural boundary condition for the method, the image will be Fourier transformed as-is. For other boundary conditions, the image will be padded before the transform is applied. The padding will extend the image by at least half the size of filter in all dimensions, and the padding will make the image size a multiple of small integers so that it is cheaper to compute the Fourier transform (see dip::OptimalFourierTransformSize). The output image will be cropped to the size of the input.

void dip::GeneralConvolution(dip::Image const& in, dip::Image const& filter, dip::Image& out, dip::StringArray const& boundaryCondition = {})

Applies a convolution with a filter kernel (PSF) by direct implementation of the convolution sum.

filter is an image, and must be equal in size or smaller than in.

As elsewhere, the origin of filter is in the middle of the image, on the pixel to the right of the center in case of an even-sized image.

Note that this is a really expensive way to compute the convolution for any filter that has more than a small amount of non-zero values. It is always advantageous to try to separate your filter into a set of 1D filters (see dip::SeparateFilter and dip::SeparableConvolution). If this is not possible, use dip::ConvolveFT with larger filters to compute the convolution in the Fourier domain.

Also, if all non-zero filter weights have the same value, dip::Uniform implements a more efficient algorithm. If filter is a binary image, dip::Uniform is called.

boundaryCondition indicates how the boundary should be expanded in each dimension. See dip::BoundaryCondition.

void dip::Convolution(dip::Image const& in, dip::Image const& filter, dip::Image& out, dip::String const& method = S::BEST, dip::StringArray const& boundaryCondition = {})

Applies a convolution with a filter kernel (PSF).

Calls either dip::SeparableConvolution, dip::ConvolveFT or dip::GeneralConvolution depending on method and the properties of filter. method can be one of:

  • "separable": Attempts to separate filter into 1D kernels using dip::SeparateFilter, and applies dip::SeparableConvolution if successful. Throws an exception if the filter is not separable.
  • "fourier": Calls dip::ConvolveFT.
  • "direct": Calls dip::GeneralConvolution.
  • "best": Uses the method that is most efficient given the sizes of in and filter, and whether filter is separable or not. It estimates the cost of each of the methods using simple models that have been fitted to timing data generated on one specific computer. These costs might not match actual costs on other machines, but form a suitable default. For applications where performance is critical, it is recommended to time the operations on the target machine, and explicitly select the best algorithm.

boundaryCondition indicates how the boundary should be expanded in each dimension. See dip::BoundaryCondition.

When calling dip::ConvolveFT, it never leaves boundaryCondition empty, to force the function to pad the image and use the same boundary condition that other methods would use. This ensures that the function doesn’t produce different results for a different choice of method. To prevent padding, call dip::ConvolveFT directly.

void dip::Uniform(dip::Image const& in, dip::Image& out, dip::Kernel const& kernel = {}, dip::StringArray const& boundaryCondition = {})

Applies a convolution with a kernel with uniform weights, leading to an average (mean) filter.

The size and shape of the kernel is given by kernel, which you can define through a default shape with corresponding sizes, or through a binary image. See dip::Kernel.

boundaryCondition indicates how the boundary should be expanded in each dimension. See dip::BoundaryCondition.

void dip::GaussFIR(dip::Image const& in, dip::Image& out, dip::FloatArray sigmas = {1.0}, dip::UnsignedArray derivativeOrder = {0}, dip::StringArray const& boundaryCondition = {}, dip::dfloat truncation = 3)

Finite impulse response implementation of the Gaussian filter and its derivatives

Convolves the image with a 1D Gaussian kernel along each dimension. For each dimension, provide a value in sigmas and derivativeOrder. The zeroth-order derivative is a plain smoothing, no derivative is computed. Derivatives with order up to 3 can be computed with this function. For higher-order derivatives, use dip::GaussFT.

The value of sigma determines the smoothing effect. For values smaller than about 0.8, the result is an increasingly poor approximation to the Gaussian filter. Use dip::GaussFT for very small sigmas. Conversely, for very large sigmas it is more efficient to use dip::GaussIIR, which runs in a constant time with respect to the sigma. Dimensions where sigma is 0 or negative are not processed, even if the derivative order is non-zero.

For the smoothing filter (derivativeOrder is 0), the size of the kernel is given by 2 * std::ceil( truncation * sigma ) + 1. The default value for truncation is 3, which assures a good approximation of the Gaussian kernel without unnecessary expense. It is possible to reduce computation slightly by decreasing this parameter, but it is not recommended. For derivatives, the value of truncation is increased by 0.5 * derivativeOrder.

boundaryCondition indicates how the boundary should be expanded in each dimension. See dip::BoundaryCondition.

void dip::GaussFT(dip::Image const& in, dip::Image& out, dip::FloatArray sigmas = {1.0}, dip::UnsignedArray derivativeOrder = {0}, dip::dfloat truncation = 3, dip::String const& inRepresentation = S::SPATIAL, dip::String const& outRepresentation = S::SPATIAL, dip::StringArray const& boundaryCondition = {})

Fourier implementation of the Gaussian filter and its derivatives

Convolves the image with a Gaussian kernel by multiplication in the Fourier domain. For each dimension, provide a value in sigmas and derivativeOrder. The value of sigma determines the smoothing effect. The zeroth-order derivative is a plain smoothing, no derivative is computed.

The values of sigmas are translated to the Fourier domain, and a Fourier-domain Gaussian is computed. Frequencies above std::ceil(( truncation + 0.5 * derivativeOrder ) * FDsigma ) are set to 0. It is a relatively minute computational difference if truncation were to be infinity, so it is not worth while to try to speed up the operation by decreasing truncation.

Dimensions where sigma is 0 or negative are not smoothed. Note that it is possible to compute a derivative without smoothing in the Fourier domain.

If in is already Fourier transformed, set inRepresentation to "frequency". Similarly, if outRepresentation is "frequency", the output will not be inverse-transformed, and so will be in the frequency domain. These two values are "spatial" by default. If any of these values is "frequency", then out will be complex, no checks are made to see if the inputs in frequency domain have the complex conjugate symmetry required for the result to be real-valued. Use dip::Image::Real if you expect the output to be real-valued in this case.

boundaryCondition indicates how the boundary should be expanded in each dimension. See dip::BoundaryCondition. The default empty boundary condition indicates no boundary extension is to be applied, the convolution will be circular (periodic boundary condition) as is natural with the DFT convolution. Specifying a boundary condition will cause the input image to be padded to a good DFT size (a product of small integers, see dip::OptimalFourierTransformSize) that is large enough to prevent visible effects of the circular convolution. Thus, specifying "periodic" as the boundary condition could, depending on the sizes of the image, speed up the operation compared to leaving the boundary condition empty.

If inRepresentation is "frequency", then boundaryCondition is ignored.

If outRepresentation is "frequency", then the padding caused by boundaryCondition will affect the output size. Leave boundaryCondition empty in this case if a predictable output size is needed.

void dip::GaussIIR(dip::Image const& in, dip::Image& out, dip::FloatArray sigmas = {1.0}, dip::UnsignedArray derivativeOrder = {0}, dip::StringArray const& boundaryCondition = {}, dip::UnsignedArray filterOrder = {}, dip::String const& designMethod = S::DISCRETE_TIME_FIT, dip::dfloat truncation = 3)

Infinite impulse response implementation of the Gaussian filter and its derivatives

Convolves the image with an IIR 1D Gaussian kernel along each dimension. For each dimension, provide a value in sigmas and derivativeOrder. The zeroth-order derivative is a plain smoothing, no derivative is computed. Derivatives with order up to 4 can be computed with this function. For higher-order derivatives, use dip::GaussFT.

The value of sigma determines the smoothing effect. For smaller values, the result is an increasingly poor approximation to the Gaussian filter. This function is efficient only for very large sigmas. Dimensions where sigma is 0 or negative are not processed, even if the derivative order is non-zero.

boundaryCondition indicates how the boundary should be expanded in each dimension. See dip::BoundaryCondition.

The filterOrder and designMethod determine how the filter is implemented. By default, designMethod is “discrete time fit”. This is the method described in van Vliet et al. (1998). filterOrder can be between 1 and 5, with 3 producing good results, and increasing order producing better results. When computing derivatives, a higher filterOrder is necessary. By default, filterOrder is 3 + derivativeOrder, capped at 5. The alternative designMethod is “forward backward”. This is the method described in Young and van Vliet (1995). Here filterOrder can be between 3 and 5.

void dip::Gauss(dip::Image const& in, dip::Image& out, dip::FloatArray sigmas = {1.0}, dip::UnsignedArray derivativeOrder = {0}, dip::String const& method = S::BEST, dip::StringArray const& boundaryCondition = {}, dip::dfloat truncation = 3)

Convolution with a Gaussian kernel and its derivatives

Convolves the image with a Gaussian kernel. For each dimension, provide a value in sigmas and derivativeOrder. The value of sigma determines the smoothing effect. The zeroth-order derivative is a plain smoothing, no derivative is computed. Dimensions where sigma is 0 or negative are not smoothed. Only the “FT” method can compute the derivative along a dimension where sigma is zero or negative.

How the convolution is computed depends on the value of method:

  • "FIR": Finite impulse response implementation, see dip::GaussFIR.
  • "IIR": Infinite impulse response implementation, see dip::GaussIIR.
  • "FT": Fourier domain implementation, see dip::GaussFT.
  • "best": Picks the best method, according to the values of sigmas and derivativeOrder:
    • if any derivativeOrder is larger than 3, use the FT method,
    • else if any sigmas is smaller than 0.8, use the FT method,
    • else if any sigmas is larger than 10, use the IIR method,
    • else use the FIR method.

boundaryCondition indicates how the boundary should be expanded in each dimension. See dip::BoundaryCondition.

void dip::FiniteDifference(dip::Image const& in, dip::Image& out, dip::UnsignedArray derivativeOrder = {0}, dip::String const& smoothFlag = S::SMOOTH, dip::StringArray const& boundaryCondition = {}, dip::BooleanArray process = {})

Finite difference derivatives

Computes derivatives using the finite difference method. Set a derivativeOrder for each dimension. Derivatives of order up to 2 can be computed with this function. The zeroth-order derivative implies either a smoothing is applied (smoothFlag == "smooth") or the dimension is not processed at all.

The smoothing filter is [1,2,1]/4 (as in the Sobel filter), the first order derivative is [1,0,-1]/2 (central difference), and the second order derivative is [1,-2,1] (which is the composition of twice the non-central difference [1,-1]). Thus, computing the first derivative twice does not yield the same result as computing the second derivative directly.

boundaryCondition indicates how the boundary should be expanded in each dimension. See dip::BoundaryCondition.

Set process to false for those dimensions that should not be filtered.

void dip::SobelGradient(dip::Image const& in, dip::Image& out, dip::uint dimension = 0, dip::StringArray const& boundaryCondition = {})

The Sobel derivative filter

This function applies the generalization of the Sobel derivative filter to arbitrary dimensions. Along the dimension dimension, the central difference is computed, and along all other dimensions, the triangular smoothing filter [1,2,1]/4 is applied.

boundaryCondition indicates how the boundary should be expanded in each dimension. See dip::BoundaryCondition.

This function calls dip::FiniteDifference.

void dip::Derivative(dip::Image const& in, dip::Image& out, dip::UnsignedArray derivativeOrder, dip::FloatArray sigmas = {1.0}, dip::String const& method = S::BEST, dip::StringArray const& boundaryCondition = {}, dip::dfloat truncation = 3)

Computes derivatives

This function provides an interface to the various derivative filters in DIPlib.

For each dimension, provide a value in sigmas and derivativeOrder. The value of sigma determines the smoothing effect. The zeroth-order derivative is a plain smoothing, no derivative is computed. If method is "best", "gaussfir" or "gaussiir", dimensions where sigma is 0 or negative are not processed, even if the derivative order is non-zero. That is, sigma must be positive for the dimension(s) where the derivative is to be computed.

method indicates which derivative filter is used:

  • "best": A Gaussian derivative, see dip::Gauss.
  • "gaussfir": The FIR implementation of the Gaussian derivative, see dip::GaussFIR.
  • "gaussiir": The IIR implementation of the Gaussian derivative, see dip::GaussIIR.
  • "gaussft": The FT implementation of the Gaussian derivative, see dip::GaussFT.
  • "finitediff": A finite difference derivative, see dip::FiniteDifference.

A finite difference derivative is an approximation to the derivative operator on the discrete grid. In contrast, convolving an image with the derivative of a Gaussian provides the exact derivative of the image convolved with a Gaussian:

\[ \frac{\partial G}{\partial x} \ast f = \frac{\partial}{\partial x}(G \ast f) \]

Thus (considering the regularization provided by the Gaussian smoothing is beneficial) it is always better to use Gaussian derivatives than finite difference derivatives.

boundaryCondition indicates how the boundary should be expanded in each dimension. See dip::BoundaryCondition.

void dip::Gradient(dip::Image const& in, dip::Image& out, dip::FloatArray sigmas = {1.0}, dip::String const& method = S::BEST, dip::StringArray const& boundaryCondition = {}, dip::BooleanArray process = {}, dip::dfloat truncation = 3)

Computes the gradient of the image, resulting in an N-vector image, if the input was N-dimensional.

Each tensor component corresponds to the first derivative along the given dimension: out[ 0 ] is the derivative along x (dimension with index 0), out[ 1 ] is the derivative along y (dimension with index 1), etc.

The input image must be scalar.

Set process to false for those dimensions along which no derivative should be taken. For example, if in is a 3D image, and process is {true,false,false}, then out will be a scalar image, containing only the derivative along the x axis.

By default uses Gaussian derivatives in the computation. Set method = "finitediff" for finite difference approximations to the gradient. See dip::Derivative for more information on the other parameters.

void dip::GradientMagnitude(dip::Image const& in, dip::Image& out, dip::FloatArray sigmas = {1.0}, dip::String const& method = S::BEST, dip::StringArray const& boundaryCondition = {}, dip::BooleanArray process = {}, dip::dfloat truncation = 3)

Computes the gradient magnitude of the image, equivalent to dip::Norm( dip::Gradient( in )).

For non-scalar images, applies the operation to each image channel. See dip::Gradient for information on the parameters.

void dip::GradientDirection(dip::Image const& in, dip::Image& out, dip::FloatArray sigmas = {1.0}, dip::String const& method = S::BEST, dip::StringArray const& boundaryCondition = {}, dip::BooleanArray process = {}, dip::dfloat truncation = 3)

Computes the direction of the gradient of the image, equivalent to dip::Angle( dip::Gradient( in )).

The input image must be scalar. For a 2D gradient, the output is scalar also, containing the angle of the gradient to the x-axis. For a 3D gradient, the output has two tensor components, containing the azimuth and inclination. See dip::Angle for an explanation.

See dip::Gradient for information on the parameters.

void dip::Curl(dip::Image const& in, dip::Image& out, dip::FloatArray sigmas = {1.0}, dip::String const& method = S::BEST, dip::StringArray const& boundaryCondition = {}, dip::BooleanArray process = {}, dip::dfloat truncation = 3)

Computes the curl (rotation) of the 2D or 3D vector field in.

Curl is defined as by \(\mathrm{curl}\,\mathbf{f} = \nabla \times \mathbf{f}\) , for a 3-vector \(\mathbf{f}=(f_x, f_y, f_z)\) (the vector image in), resulting in a 3-vector with components:

\begin{eqnarray*} (\mathrm{curl}\,\mathbf{f})_x &=& \frac{\partial}{\partial y} f_z - \frac{\partial}{\partial z} f_y \, , \\ (\mathrm{curl}\,\mathbf{f})_y &=& \frac{\partial}{\partial z} f_x - \frac{\partial}{\partial x} f_z \, , \\ (\mathrm{curl}\,\mathbf{f})_z &=& \frac{\partial}{\partial x} f_y - \frac{\partial}{\partial y} f_x \, . \end{eqnarray*}

For the 2D case, \(f_z\) is assumed to be zero, and only the z-component of the curl is computed, yielding a scalar output.

in is expected to be a 2D or 3D image with a 2-vector or a 3-vector tensor representation, respectively. However, the image can have more dimensions if they are excluded from processing through process. See dip::Gradient for information on the parameters.

void dip::Divergence(dip::Image const& in, dip::Image& out, dip::FloatArray sigmas = {1.0}, dip::String const& method = S::BEST, dip::StringArray const& boundaryCondition = {}, dip::BooleanArray process = {}, dip::dfloat truncation = 3)

Computes the divergence of the vector field in.

Divergence is defined as

\[ \mathrm{div}\,\mathbf{f} = \nabla \cdot \mathbf{f} = \frac{\partial}{\partial x}f_x + \frac{\partial}{\partial y}f_y + \frac{\partial}{\partial z}f_z \; , \]

with \(\mathbf{f}=(f_x, f_y, f_z)\) the vector image in. This concept naturally extends to any number of dimensions.

in is expected to have as many dimensions as tensor components. However, the image can have more dimensions if they are excluded from processing through process. See dip::Gradient for information on the parameters.

void dip::Hessian(dip::Image const& in, dip::Image& out, dip::FloatArray sigmas = {1.0}, dip::String const& method = S::BEST, dip::StringArray const& boundaryCondition = {}, dip::BooleanArray process = {}, dip::dfloat truncation = 3)

Computes the Hessian of the image, resulting in a symmetric NxN tensor image, if the input was N-dimensional.

The Hessian of input image \(f\) is given by \(\mathbf{H} = \nabla \nabla^T f\) , with tensor components

\[ \mathbf{H}_{i,j} = \frac{\partial^2}{\partial u_i \partial u_j} f \; . \]

Each tensor component corresponds to one of the second-order derivatives. Note that \(H\) is a symmetric matrix (order of differentiation does not matter). Duplicate entries are not stored in the symmetric tensor image.

Image dimensions for which process is false do not participate in the set of dimensions that form the Hessian matrix. Thus, a 5D image with only two dimensions selected by the process array will yield a 2-by-2 Hessian matrix.

By default this function uses Gaussian derivatives in the computation. Set method = "finitediff" for finite difference approximations to the gradient. See dip::Derivative for more information on the other parameters.

The input image must be scalar.

void dip::Laplace(dip::Image const& in, dip::Image& out, dip::FloatArray sigmas = {1.0}, dip::String const& method = S::BEST, dip::StringArray const& boundaryCondition = {}, dip::BooleanArray process = {}, dip::dfloat truncation = 3)

Computes the Laplacian of the image, equivalent to dip::Trace( dip::Hessian( in )), but more efficient.

The Laplacian of input image \(f\) is written as \(\nabla\cdot\nabla f = \nabla^2 f = \Delta f\) , and given by

\[ \Delta f = \sum_i \frac{\partial^2}{\partial u_i^2} f \; . \]

See dip::Gradient for information on the parameters.

If method is “finitediff”, it does not add second order derivatives, but instead computes a convolution with a 3x3(x3x…) kernel where all elements are -1 and the middle element is \(3^d - 1\) (with \(d\) the number of image dimensions). That is, the kernel sums to 0. For a 2D image, this translates to the well-known kernel:

\[ \begin{bmatrix} -1 & -1 & -1 \\ -1 & 8 & -1 \\ -1 & -1 & -1 \end{bmatrix} \]

void dip::Dgg(dip::Image const& in, dip::Image& out, dip::FloatArray sigmas = {1.0}, dip::String const& method = S::BEST, dip::StringArray const& boundaryCondition = {}, dip::BooleanArray process = {}, dip::dfloat truncation = 3)

Computes the second derivative in the gradient direction.

The second derivative in the gradient direction is computed by Raleigh quotient of the Hessian matrix and the gradient vector:

\[ f_{gg} = \frac{ \nabla^T \! f \; \nabla \nabla^T \! f \; \nabla f } { \nabla^T \! f \; \nabla f } \; . \]

This function is equivalent to:

Image g = dip::Gradient( in, ... );
Image H = dip::Hessian( in, ... );
Image Dgg = dip::Transpose( g ) * H * g;
Dgg /= dip::Transpose( g ) * g;

See dip::Derivative for how derivatives are computed, and the meaning of the parameters. See dip::Gradient or dip::Hessian for the meaning of the process parameter

void dip::LaplacePlusDgg(dip::Image const& in, dip::Image& out, dip::FloatArray sigmas = {1.0}, dip::String const& method = S::BEST, dip::StringArray const& boundaryCondition = {}, dip::BooleanArray process = {}, dip::dfloat truncation = 3)

Adds the second derivative in the gradient direction to the Laplacian.

This function computes dip::Laplace( in ) + dip::Dgg( in ), but avoiding computing the second derivatives twice.

The zero-crossings of the result correspond to the edges in the image, just as they do for the individual Laplace and Dgg operators. However, the localization is improved by an order of magnitude with respect to the individual operators.

See dip::Laplace and dip::Dgg for more information.

void dip::LaplaceMinusDgg(dip::Image const& in, dip::Image& out, dip::FloatArray sigmas = {1.0}, dip::String const& method = S::BEST, dip::StringArray const& boundaryCondition = {}, dip::BooleanArray process = {}, dip::dfloat truncation = 3)

Subtracts the second derivative in the gradient direction from the Laplacian.

This function computes dip::Laplace( in ) - dip::Dgg( in ), but avoiding computing the second derivatives twice.

For two-dimensional images, this is equivalent to the second order derivative in the direction perpendicular to the gradient direction.

See dip::Laplace and dip::Dgg for more information.

void dip::Sharpen(dip::Image const& in, dip::Image& out, dip::dfloat weight = 1.0, dip::FloatArray sigmas = {1.0}, dip::String const& method = S::BEST, dip::StringArray const& boundaryCondition = {}, dip::dfloat truncation = 3)

Sharpens in by subtracting the Laplacian of the image.

The actual operation applied is:

out = in - dip::Laplace( in ) * weight;

See dip::Laplace and dip::Gradient for information on the parameters.

void dip::UnsharpMask(dip::Image const& in, dip::Image& out, dip::dfloat weight = 1.0, dip::FloatArray sigmas = {1.0}, dip::String const& method = S::BEST, dip::StringArray const& boundaryCondition = {}, dip::dfloat truncation = 3)

Sharpens in by subtracting the smoothed image.

The actual operation applied is:

out = in * ( 1+weight ) - dip::Gauss( in ) * weight;

See dip::Gauss and dip::Gradient for information on the parameters.

void dip::GaborFIR(dip::Image const& in, dip::Image& out, dip::FloatArray sigmas, dip::FloatArray const& frequencies, dip::StringArray const& boundaryCondition = {}, dip::BooleanArray process = {}, dip::dfloat truncation = 3)

Finite impulse response implementation of the Gabor filter

Convolves the image with an FIR 1D Gabor kernel along each dimension. For each dimension, provide a value in sigmas and frequencies. The value of sigma determines the amount of local averaging. For smaller values, the result is more precise spatially, but less selective of frequencies. Dimensions where sigma is 0 or negative are not processed.

Frequencies are in the range [0, 0.5), with 0.5 being the frequency corresponding to a period of 2 pixels.

The output is complex-valued. Typically, the magnitude is the interesting part of the result.

boundaryCondition indicates how the boundary should be expanded in each dimension. See dip::BoundaryCondition.

Set process to false for those dimensions that should not be filtered. This is equivalent to setting sigmas to 0 for those dimensions.

This function is relatively slow compared to dip::GaborIIR, even for small sigmas. Prefer to use the IIR implementation.

void dip::GaborIIR(dip::Image const& in, dip::Image& out, dip::FloatArray sigmas, dip::FloatArray const& frequencies, dip::StringArray const& boundaryCondition = {}, dip::BooleanArray process = {}, dip::IntegerArray const& filterOrder = {}, dip::dfloat truncation = 3)

Recursive infinite impulse response implementation of the Gabor filter

Convolves the image with an IIR 1D Gabor kernel along each dimension. For each dimension, provide a value in sigmas and frequencies. The value of sigma determines the amount of local averaging. For smaller values, the result is more precise spatially, but less selective of frequencies. Dimensions where sigma is 0 or negative are not processed.

Frequencies are in the range [0, 0.5), with 0.5 being the frequency corresponding to a period of 2 pixels.

The output is complex-valued. Typically, the magnitude is the interesting part of the result.

boundaryCondition indicates how the boundary should be expanded in each dimension. See dip::BoundaryCondition.

Set process to false for those dimensions that should not be filtered. This is equivalent to setting sigmas to 0 for those dimensions.

void dip::Gabor2D(dip::Image const& in, dip::Image& out, dip::FloatArray sigmas = {5.0,5.0}, dip::dfloat frequency = 0.1, dip::dfloat direction = dip::pi, dip::StringArray const& boundaryCondition = {}, dip::dfloat truncation = 3)

2D Gabor filter with direction parameter

Convolves the 2D image with a Gabor kernel. This is a convenience wrapper around dip::GaborIIR. The value of sigma determines the amount of local averaging, and can be different for each dimension. For smaller values, the result is more precise spatially, but less selective of frequencies.

frequency is in the range [0, 0.5), with 0.5 being the frequency corresponding to a period of 2 pixels. direction is the filter direction, in the range [0, 2π].

The output is complex-valued. Typically, the magnitude is the interesting part of the result.

boundaryCondition indicates how the boundary should be expanded in each dimension. See dip::BoundaryCondition.

To use cartesian frequency coordinates, see dip::GaborIIR.

void dip::LogGaborFilterBank(dip::Image const& in, dip::Image& out, dip::FloatArray const& wavelengths = {3.0,6.0,12.0,24.0}, dip::dfloat bandwidth = 0.75, dip::uint nOrientations = 6, dip::String const& inRepresentation = S::SPATIAL, dip::String const& outRepresentation = S::SPATIAL)

Applies a log-Gabor filter bank

A log-Gabor filter is a Gabor filter computed on the logarithm of the frequency, leading to a shorter tail of the Gaussian, in the frequency domain, towards the lower frequencies. The origin (DC component) is thus never included in the filter. This gives it better scale selection properties than the traditional Gabor filter.

This function generates a filter bank with wavelengths.size() times nOrientations filters. The width of the filters in the angular axis is determined by the number of orientations used, and their locations are always equally distributed over π radian, starting at 0. The radial location (scales) of the filters is determined by wavelengths (in pixels), which determines the center for each scale filter. The widths of the filters in this direction are determined by the bandwidth parameter; the default value of 0.75 corresponds approximately to one octave, 0.55 to two octaves, and 0.41 to three octaves.

wavelengths.size() and nOrientations must be at least 1. If nOrientations is 1, no orientation filtering is applied, the filters become purely real. These filters can be defined for images of any dimensionality. For more than one orientation, the filters are complex-valued in the spatial domain, and can only be created for 2D images. See dip::MonogenicSignal for a generalization to arbitrary dimensionality.

If in is not forged, its sizes will be used to generate the filters, which will be returned. Thus, this is identical to (but slightly cheaper than) using a delta pulse image as input.

The filters are always generated in the frequency domain. If outRepresentation is "spatial", the inverse Fourier transform will be applied to bring the result back to the spatial domain. Otherwise, it should be "frequency", and no inverse transform will be applied. Likewise, inRepresentation specifies whether in has already been converted to the frequency domain or not.

Out will be a tensor image with wavelengths.size() tensor rows and nFrequencyScales tensor columns. The data type will be either single-precision float or single-precision complex, depending on the selected parameters.

void dip::NormalizedConvolution(dip::Image const& in, dip::Image const& mask, dip::Image& out, dip::FloatArray const& sigmas = {1.0}, dip::String const& method = S::BEST, dip::StringArray const& boundaryCondition = {S::ADD_ZEROS}, dip::dfloat truncation = 3)

Computes the normalized convolution with a Gaussian kernel: a Gaussian convolution for missing or uncertain data.

The normalized convolution is a convolution that handles missing or uncertain data. mask is an image, expected to be in the range [0,1], that indicates the confidence in each of the values of in. Missing values are indicated by setting the corresponding value in mask to 0.

The normalized convolution is then Convolution( in * mask ) / Convolution( mask ).

This function applies convolutions with a Gaussian kernel, using dip::Gauss. See that function for the meaning of the parameters. boundaryCondition defaults to "add zeros", the normalized convolution then takes pixels outside of the image domain as missing values.

void dip::NormalizedDifferentialConvolution(dip::Image const& in, dip::Image const& mask, dip::Image& out, dip::uint dimension = 0, dip::FloatArray const& sigmas = {1.0}, dip::String const& method = S::BEST, dip::StringArray const& boundaryCondition = {S::ADD_ZEROS}, dip::dfloat truncation = 3)

Computes the normalized differential convolution with a Gaussian kernel: a derivative operator for missing or uncertain data.

The normalized convolution is a convolution that handles missing or uncertain data. mask is an image, expected to be in the range [0,1], that indicates the confidence in each of the values of in. Missing values are indicated by setting the corresponding value in mask to 0.

The normalized differential convolution is defined here as the derivative of the normalized convolution with a Gaussian kernel:

\[ \frac{\partial}{\partial x} \frac{(f \, m) \ast g}{m \ast g} = \frac{(f \, m) \ast \frac{\partial}{\partial x} g}{m \ast g} - \frac{(f \, m) \ast g}{m \ast g} \frac{m \ast \frac{\partial}{\partial x} g}{m \ast g} \; . \]

\(\ast\) is the convolution operator, \(f\) is in, \(m\) is mask, and \(g\) is the Gaussian kernel

The derivative is computed along dimension.

This function uses dip::Gauss. See that function for the meaning of the parameters. boundaryCondition defaults to "add zeros", the normalized convolution then takes pixels outside of the image domain as missing values.

void dip::MeanShiftVector(dip::Image const& in, dip::Image& out, dip::FloatArray sigmas = {1.0}, dip::String const& method = S::BEST, dip::StringArray const& boundaryCondition = {}, dip::dfloat truncation = 3)

Computes the mean shift vector for each pixel in the image

The output image is a vector image, indicating the step to take to move the window center to its center of mass. Repeatedly following the vector will lead to a local maximum of the image in. in must be scalar and real-valued.

The mean shift at a given location \(x\) is then given by

\[ s = \frac{ \sum_i{(x-x_i) w(x-x_i) f(x_i)} }{ \sum_i{w(x-x_i) f(x_i)} } = \frac{ \left ( -x w \right) \ast f }{ w \ast f } \; , \]

where \(f\) is the image, \(w\) is a windowing function, and \(\ast\) indicates convolution.

We use a Gaussian window with sizes given by sigmas. A Gaussian window causes slower convergence than a uniform window, but yields a smooth trajectory and more precise results (according to Comaniciu and Meer, 2002). It also allows us to rewrite the above (with \(g_\sigma\) the Gaussian window with parameter \(\sigma\) ) as

\[ s = \frac{ \left ( -x g_\sigma \right) \ast f }{ g_\sigma \ast f } = \frac{ \left ( \sigma^2 \nabla g_\sigma \right) \ast f }{ g_\sigma \ast f } \; . \]

Thus, we can write this filter as dip::Gradient(in, sigmas) / dip::Gauss(in, sigmas) * sigmas * sigmas. See dip::Derivative for more information on the parameters. Do not use method = "finitediff", as it will lead to nonsensical results.