Analysis module #include "diplib/analysis.h"
Assorted analysis tools.
Classes

struct dip::
SubpixelLocationResult  Contains the result of the function
dip::SubpixelLocation
. 
class dip::
Distribution  Holds probability density functions and other types of distribution
Aliases

using dip::
SubpixelLocationArray = std::vector<SubpixelLocationResult>  Contains the result of the functions
dip::SubpixelMaxima
anddip::SubpixelMinima
.
Functions

auto dip::
Find (dip::Image const& in, dip::Image const& mask = {}) > dip::CoordinateArray  Finds the coordinates for all nonzero pixels in the scalar image
in
, optionally constrained to the pixels selected bymask
. 
auto dip::
SubpixelLocation (dip::Image const& in, dip::UnsignedArray const& position, dip::String const& polarity = S::MAXIMUM, dip::String const& method = S::PARABOLIC_SEPARABLE) > dip::SubpixelLocationResult  Gets coordinates of a local extremum with subpixel precision

auto dip::
SubpixelMaxima (dip::Image const& in, dip::Image const& mask = {}, dip::String const& method = S::PARABOLIC_SEPARABLE) > dip::SubpixelLocationArray  Gets coordinates of local maxima with subpixel precision

auto dip::
SubpixelMinima (dip::Image const& in, dip::Image const& mask = {}, dip::String const& method = S::PARABOLIC_SEPARABLE) > dip::SubpixelLocationArray  Gets coordinates of local minima with subpixel precision

auto dip::
MeanShift (dip::Image const& meanShiftVectorResult, dip::FloatArray const& start, dip::dfloat epsilon = 1e3) > dip::FloatArray  Finds the coordinates of a local maximum close to
start

auto dip::
MeanShift (dip::Image const& meanShiftVectorResult, dip::FloatCoordinateArray const& startArray, dip::dfloat epsilon = 1e3) > dip::FloatCoordinateArray  Finds the coordinates of local a maximum close to each point in
startArray
. 
void dip::
GaussianMixtureModel (dip::Image const& in, dip::Image& out, dip::uint dimension = 2, dip::uint numberOfGaussians = 2, dip::uint maxIter = 20, dip::StringSet const& flags = {})  Determines the parameters for a Gaussian Mixture Model for every line in the image.

void dip::
CrossCorrelationFT (dip::Image const& in1, dip::Image const& in2, dip::Image& out, dip::String const& in1Representation = S::SPATIAL, dip::String const& in2Representation = S::SPATIAL, dip::String const& outRepresentation = S::SPATIAL, dip::String const& normalize = S::NORMALIZE)  Calculates the crosscorrelation between two images of equal size.

void dip::
AutoCorrelationFT (dip::Image const& in, dip::Image& out, dip::String const& inRepresentation = S::SPATIAL, dip::String const& outRepresentation = S::SPATIAL)  Computes the autocorrelation function.

auto dip::
FindShift (dip::Image const& in1, dip::Image const& in2, dip::String const& method = S::MTS, dip::dfloat parameter = 0, dip::UnsignedArray maxShift = {}) > dip::FloatArray  Estimates the (subpixel) global shift between
in1
andin2
. 
auto dip::
FourierMellinMatch2D (dip::Image const& in1, dip::Image const& in2, dip::Image& out, dip::String const& interpolationMethod = S::LINEAR, dip::String const& correlationMethod = S::PHASE) > dip::FloatArray  Finds the scaling, translation and rotation between two 2D images using the Fourier Mellin transform

void dip::
StructureTensor (dip::Image const& in, dip::Image const& mask, dip::Image& out, dip::FloatArray const& gradientSigmas = {1.0}, dip::FloatArray const& tensorSigmas = {5.0}, dip::String const& method = S::BEST, dip::StringArray const& boundaryCondition = {}, dip::dfloat truncation = 3)  Computes the structure tensor.

void dip::
StructureTensorAnalysis2D (dip::Image const& in, dip::Image* l1 = nullptr, dip::Image* l2 = nullptr, dip::Image* orientation = nullptr, dip::Image* energy = nullptr, dip::Image* anisotropy1 = nullptr, dip::Image* anisotropy2 = nullptr, dip::Image* curvature = nullptr)  Computes useful image parameters from the 2D structure tensor.

void dip::
StructureTensorAnalysis3D (dip::Image const& in, dip::Image* l1 = nullptr, dip::Image* phi1 = nullptr, dip::Image* theta1 = nullptr, dip::Image* l2 = nullptr, dip::Image* phi2 = nullptr, dip::Image* theta2 = nullptr, dip::Image* l3 = nullptr, dip::Image* phi3 = nullptr, dip::Image* theta3 = nullptr, dip::Image* energy = nullptr, dip::Image* cylindrical = nullptr, dip::Image* planar = nullptr)  Computes useful image parameters from the 3D structure tensor.

void dip::
StructureTensorAnalysis (dip::Image const& in, dip::ImageRefArray& out, dip::StringArray const& outputs)  Interface to
dip::StructureTensorAnalysis2D
anddip::StructureTensorAnalysis3D
. 
auto dip::
StructureAnalysis (dip::Image const& in, dip::Image const& mask, std::vector<dfloat> const& scales = {}, dip::String const& feature = "energy", dip::FloatArray const& gradientSigmas = {1.0}, dip::String const& method = S::BEST, dip::StringArray const& boundaryCondition = {}, dip::dfloat truncation = 3) > dip::Distribution  Analyzes the local structure of the image at multiple scales.

void dip::
MonogenicSignal (dip::Image const& in, dip::Image& out, dip::FloatArray const& wavelengths = {3.0,24.0}, dip::dfloat bandwidth = 0.41, dip::String const& inRepresentation = S::SPATIAL, dip::String const& outRepresentation = S::SPATIAL)  Computes the monogenic signal, a multidimensional generalization of the analytic signal.

void dip::
MonogenicSignalAnalysis (dip::Image const& in, dip::ImageRefArray& out, dip::StringArray const& outputs, dip::dfloat noiseThreshold = 0.2, dip::dfloat frequencySpreadThreshold = 0.5, dip::dfloat sigmoidParameter = 10, dip::dfloat deviationGain = 1.5, dip::String const& polarity = S::BOTH)  Computes useful image parameters from the monogenic signal.

void dip::
OrientationSpace (dip::Image const& in, dip::Image& out, dip::uint order = 8, dip::dfloat radCenter = 0.1, dip::dfloat radSigma = 0.8, dip::uint orientations = 0)  Creates an orientation space for a 2D image

auto dip::
PairCorrelation (dip::Image const& object, dip::Image const& mask, dip::Random& random, dip::uint probes = 1000000, dip::uint length = 100, dip::String const& sampling = S::RANDOM, dip::StringSet const& options = {}) > dip::Distribution  Estimates the pair correlation function of the different phases in
object
. 
auto dip::
PairCorrelation (dip::Image const& object, dip::Image const& mask, dip::uint probes = 1000000, dip::uint length = 100, dip::String const& sampling = S::RANDOM, dip::StringSet const& options = {}) > dip::Distribution  like above, using a defaultinitialized
dip::Random
object. 
auto dip::
ProbabilisticPairCorrelation (dip::Image const& phases, dip::Image const& mask, dip::Random& random, dip::uint probes = 1000000, dip::uint length = 100, dip::String const& sampling = S::RANDOM, dip::StringSet const& options = {}) > dip::Distribution  Estimates the probabilistic pair correlation function of the different phases in
phases
. 
auto dip::
ProbabilisticPairCorrelation (dip::Image const& object, dip::Image const& mask, dip::uint probes = 1000000, dip::uint length = 100, dip::String const& sampling = S::RANDOM, dip::StringSet const& options = {}) > dip::Distribution  like above, using a defaultinitialized
dip::Random
object. 
auto dip::
Semivariogram (dip::Image const& in, dip::Image const& mask, dip::Random& random, dip::uint probes = 1000000, dip::uint length = 100, dip::String const& sampling = S::RANDOM) > dip::Distribution  Estimates the expected value of half the square difference between field values at a distance
d
. 
auto dip::
Semivariogram (dip::Image const& object, dip::Image const& mask, dip::uint probes = 1000000, dip::uint length = 100, dip::String const& sampling = S::RANDOM) > dip::Distribution  like above, using a defaultinitialized
dip::Random
object. 
auto dip::
ChordLength (dip::Image const& object, dip::Image const& mask, dip::Random& random, dip::uint probes = 100000, dip::uint length = 100, dip::String const& sampling = S::RANDOM) > dip::Distribution  Estimates the chord length distribution of the different phases in
object
. 
auto dip::
ChordLength (dip::Image const& object, dip::Image const& mask, dip::uint probes = 100000, dip::uint length = 100, dip::String const& sampling = S::RANDOM) > dip::Distribution  like above, using a defaultinitialized
dip::Random
object. 
auto dip::
DistanceDistribution (dip::Image const& object, dip::Image const& region, dip::uint length = 100) > dip::Distribution  Computes the distribution of distances to the background of
region
for the different phases inobject
. 
auto dip::
Granulometry (dip::Image const& in, dip::Image const& mask, std::vector<dfloat> const& scales = {}, dip::String const& type = "isotropic", dip::String const& polarity = S::OPENING, dip::StringSet const& options = {}) > dip::Distribution  Computes the granulometric function for an image

auto dip::
FractalDimension (dip::Image const& in, dip::dfloat eta = 0.5) > dip::dfloat  Estimates the fractal dimension of the binary image
in
the sliding box method. 
void dip::
swap (Sample& a, Sample& b)  Swaps two samples, copying the data from
other
to*this
, and that from*this
toother
. Both must have the same number of values.
Operators

auto dip::
operator<< (std::ostream& os, dip::Distribution const& distribution) > std::ostream&  Writes the distribution to a stream
Class documentation
struct dip::SubpixelLocationResult
Contains the result of the function dip::SubpixelLocation
.
Variables  

dip::FloatArray coordinates  Coordinates of local extremum 
dip::dfloat value  Interpolated value at local extremum 
Function documentation
dip::CoordinateArray
dip::Find (dip::Image const& in,
dip::Image const& mask = {})
Finds the coordinates for all nonzero pixels in the scalar image in
, optionally constrained to the pixels selected by mask
.
in
must be scalar, but can have any data type. mask
, if forged, must be of the same sizes as in
, or be
singleton expandable to that size, and must be binary.
The output array contains the coordinates for all nonzero pixels, in linear index order (that is, the image is traversed rowwise to find the pixels). Note that this function is intended to be used with images that have relatively few nonzero pixels, and there usually is a better alternative than to list the coordinates of all nonzero pixels.
dip::SubpixelLocationResult
dip::SubpixelLocation (dip::Image const& in,
dip::UnsignedArray const& position,
dip::String const& polarity = S::MAXIMUM,
dip::String const& method = S::PARABOLIC_SEPARABLE)
Gets coordinates of a local extremum with subpixel precision
Determines the subpixel location of a local maximum or minimum close to position
. position
should point
to a pixel whose value is larger than its direct neighbors’ (if polarity
is "maximum"
) or smaller than
its direct neighbors’ (polarity
is "minimum"
).
The coordinates
member of the output struct will contain the the subpixel location of this local
extremum. The value
member will contain the interpolated grey value at the location of the extremum.
method
determines which method is used. These are the allowed values:
"linear"
: Computes the center of gravity of 3 pixels around the extremum, in each dimension independently. The value at the extremum is that of the pixel atposition
."parabolic"
: Fits a parabola to a 3pixelwide block around the extremum (for up to 3D only)."parabolic separable"
: Fits a parabola in each dimension independently. The value at the extremum is the maximum/minimum of the values obtained in each of the 1D processes, and thus not equivalent to the grey value obtained by true interpolation."gaussian"
: Fits a Gaussian to a 3pixelwide block around the extremum (for up to 3D only)."gaussian separable"
: Fits a Gaussian in each dimension independently. The value at the extremum is the maximum/minimum of the values obtained in each of the 1D processes, and thus not equivalent to the grey value obtained by true interpolation."integer"
: Doesn’t look for subpixel locations, returning the integer coordinates of the extremum.
The image in
must be scalar and realvalued.
dip::SubpixelLocationArray
dip::SubpixelMaxima (dip::Image const& in,
dip::Image const& mask = {},
dip::String const& method = S::PARABOLIC_SEPARABLE)
Gets coordinates of local maxima with subpixel precision
Detects local maxima in the image, and returns their coordinates, with subpixel precision.
Only pixels where mask
is on will be examined. Local maxima are detected using dip::Maxima
,
then their position is determined more precisely using dip::SubpixelLocation
.
A local maximum can not touch the edge of the image. That is, its integer location must be one pixel away from the edge.
See dip::SubpixelLocation
for the definition of the method
parameter.
dip::SubpixelLocationArray
dip::SubpixelMinima (dip::Image const& in,
dip::Image const& mask = {},
dip::String const& method = S::PARABOLIC_SEPARABLE)
Gets coordinates of local minima with subpixel precision
Detects local minima in the image, and returns their coordinates, with subpixel precision.
Only pixels where mask
is on will be examined. Local minima are detected using dip::Minima
,
then their position is determined more precisely using dip::SubpixelLocation
.
A local minimum can not touch the edge of the image. That is, its integer location must be one pixel away from the edge.
See dip::SubpixelLocation
for the definition of the method
parameter.
dip::FloatArray
dip::MeanShift (dip::Image const& meanShiftVectorResult,
dip::FloatArray const& start,
dip::dfloat epsilon = 1e3)
Finds the coordinates of a local maximum close to start
The mean shift method iteratively finds the center of gravity of a Gaussianshaped neighborhood around start
,
moving start
to this location, until convergence. The location that the process converges to is a local
maximum. Convergence is assumed when the magnitude of the shift is smaller than epsilon
.
To speed up the computation within this function, the output of dip::MeanShiftVector
is used. This filter
precomputes the center of gravity of all neighborhoods in the image. Specify the neighborhood sizes in that
function.
dip::FloatCoordinateArray
dip::MeanShift (dip::Image const& meanShiftVectorResult,
dip::FloatCoordinateArray const& startArray,
dip::dfloat epsilon = 1e3)
Finds the coordinates of local a maximum close to each point in startArray
.
Repeatedly calls the dip::MeanShift
function described above, for each point in startArray
. Duplicate maxima
are not removed, such that out[ii]
is the local maximum arrived to from startArray[ii]
.
void
dip::GaussianMixtureModel (dip::Image const& in,
dip::Image& out,
dip::uint dimension = 2,
dip::uint numberOfGaussians = 2,
dip::uint maxIter = 20,
dip::StringSet const& flags = {})
Determines the parameters for a Gaussian Mixture Model for every line in the image.
in
must be realvalued and scalar. For each image line along dimension
, a Gaussian Mixture Model (GMM)
is fitted to find the numberOfGaussians
most important Gaussian components. In the output image, dimension
indexes the Gaussian component. The components are stored in order of importance, the component with the largest
amplitude first. The parameters to each Gaussian are stored as tensor components, in the order
position, amplitude and sigma.
out
is a doubleprecision floatingpoint image with three tensor components. It has a size of numberOfGaussians
along dimension
, and the same size as in
along all other dimensions.
maxIter
sets how many iterations are run. There is currently no other stopping criterion.
flags
can contain the following values:
"periodic"
specifies that the image is periodic alongdimension
."pixel size"
scales the parameters written to the output image with the pixel size alongdimension
. By default these parameters are in pixels.
void
dip::CrossCorrelationFT (dip::Image const& in1,
dip::Image const& in2,
dip::Image& out,
dip::String const& in1Representation = S::SPATIAL,
dip::String const& in2Representation = S::SPATIAL,
dip::String const& outRepresentation = S::SPATIAL,
dip::String const& normalize = S::NORMALIZE)
Calculates the crosscorrelation between two images of equal size.
The normalize
flag determines whether to compute the regular crosscorrelation ("don't normalize"
)
or to normalize the frequencydomain inputs so that only the phase information is of importance
("normalize"
or "phase"
). The normalization results as a very sharp peak in the spatial domain, which
increases the precision with which a shift between the two images can be determined.
Note that this normalization is not related to what is commonly referred to as “normalized crosscorrelation”,
where the input images are whitened (in the spatial domain) before the crosscorrelations is computed. The
method is instead related to the “phase correlation” as proposed by Kuglin and Hines (1975).
When normalize
is "normalize"
, the computation performed in the Fourier Domain is (as per Luengo & van Vliet, 2000)
out = ( in1 * dip::Conjugate(in2) ) / dip::SquareModulus(in1);
When normalize
is set to "phase"
, the computation is (as per Kuglin & Hines, 1975)
out = ( in1 * dip::Conjugate(in2) ) / ( dip::Modulus(in1) * dip::Modulus(in2) );
(See dip::Conjugate
, dip::SquareModulus
and dip::Modulus
.)
These two approaches are identical if the only difference between in1
and in2
is a shift, except that
"normalize"
is computationally cheaper than "phase"
. If the images are otherwise not identical, then
the "phase"
method might be the better choice. If there exist strong differences, the nonnormalized
method can be the best option.
As elsewhere, the origin is in the middle of the image, on the pixel to the right of
the center in case of an evensized image. Thus, for in1==in2
, only this pixel will be set.
See dip::FindShift
with the "CC"
, "NCC"
or "PC"
methods for localizing this peak.
If in1
or in2
is already Fourier transformed, set in1Representation
or in2Representation
to "frequency"
. Similarly, if outRepresentation
is "frequency"
, the output will not be
inversetransformed, so will be in the frequency domain.
in1
and in2
must be scalar images with the same dimensionality and sizes. Spatialdomain
images must be realvalued, and frequencydomain images are assumed (but not tested) to be
conjugatesymmetric (i.e. be Fourier transforms of realvalued images).
out
will be realvalued if outRepresentation
is "spatial"
.
void
dip::AutoCorrelationFT (dip::Image const& in,
dip::Image& out,
dip::String const& inRepresentation = S::SPATIAL,
dip::String const& outRepresentation = S::SPATIAL)
Computes the autocorrelation function.
Equivalent to calling dip::CrossCorrelationFT
with the same image as the two inputs, and the
normalize
flag set to "don't normalize"
(normalization would yield a very boring output indeed).
As elsewhere, the origin is in the middle of the image, on the pixel to the right of the center in case of an evensized image.
If in
is already Fourier transformed, set inRepresentation
to "frequency"
.
Similarly, if outRepresentation
is "frequency"
, the output will not be
inversetransformed, so will be in the frequency domain.
in
must be a scalar image. Spatialdomain images must be realvalued, and a frequencydomain image
is assumed (but not tested) to be conjugatesymmetric (i.e. the Fourier transform of realvalued image).
out
will be realvalued if outRepresentation
is "spatial"
.
dip::FloatArray
dip::FindShift (dip::Image const& in1,
dip::Image const& in2,
dip::String const& method = S::MTS,
dip::dfloat parameter = 0,
dip::UnsignedArray maxShift = {})
Estimates the (subpixel) global shift between in1
and in2
.
The numbers found represent the shift of in2
with respect to in1
, or equivalently, the position of the
topleft pixel of in2
in the coordinate system of in1
. Shifting in1
by the returned shift aligns the
two images. See dip::Shift
.
There are various methods that can be used to determine the shift, see below.
For the methods that require that the shift be small, first the integer pixel is calculated using cross
correlation, and both images are cropped to the common part. maxShift
can be used to restrict the shift
found. This is useful, for example, when the images contain a repeating pattern. Noise in the images can
cause a shift of several pattern periods to be “optimal” (e.g. the crosscorrelation yields a larger value),
even if a much smaller shift would also explain the differences. In this case, set maxShift
to be slightly
smaller than the pattern period along each dimension.
Valid strings for method
are:

"integer only"
: The crosscorrelation method simply computes the crosscorrelation, usingdip::CrossCorrelationFT
, and then finds the position of the pixel with the largest value. This method works for any number of dimensions.parameter
is ignored. 
"CC"
: The crosscorrelation method computes the crosscorrelation, usingdip::CrossCorrelationFT
, and then usesdip::SubpixelLocation
to find the location of the largest peak with subpixel precision. This method works for any number of dimensions.parameter
is ignored. 
"NCC"
: As"CC"
, but using the normalized crosscorrelation, which makes the peak much sharper (Luengo Hendriks, 1998). This method works for any number of dimensions.parameter
is ignored. See the notes indip::CrossCorrelationFT
regarding the normalization of the crosscorrelation using the"normalize"
flag, which is not as what is commonly referred to as “normalized crosscorrelation”. 
"PC"
: As"CC"
, but using phase correlation. This method works for any number of dimensions.parameter
is ignored. See the notes indip::CrossCorrelationFT
regarding the normalization of the crosscorrelation using the"phase"
flag. 
"CPF"
: The CPF method (see Luengo Hendriks (1998), where it is called FFTS) uses the phase of the crosscorrelation (as calculated bydip::CrossCorrelationFT
) to estimate the shift.parameter
sets the largest frequency used in this estimation. The maximum value that makes sense issqrt(0.5)
. Any larger value will give the same result. Choose smaller values to ignore the higher frequencies, which have a smaller SNR and are more affected by aliasing. Ifparameter
is <= 0, the optimal found for images subsampled by a factor four will be used (parameter = 0.2). This method only supports 2D images. 
"MTS"
: The MTS method (see Luengo Hendriks (1998), where it is called GRS) uses a first order Taylor approximation of the equationin1(t) = in2(ts)
at scaleparameter
. Ifparameter
is <= 0, a scale of 1 will be used. This means that the images will be smoothed with a Gaussian kernel of 1. This method is more accurate than CPF. This method supports images with a dimensionality between 1 and 3. 
"ITER"
: The ITER method is an iterative version of the MTS method. It is known that a single estimation with MTS has a bias due to truncation of the Taylor expansion series (Pham et al., 2005). The bias can be expressed as a polynomial of the subpixel displacements. As a result, if the MTS method is applied iteratively, and the shift is refined after each iteration, the bias eventually becomes negligible. By using just 3 iterations, and noticing thatlog(bias_increment)
is a linear sequence, it is possible to correct for the bias up to O(10^{6}).
Set parameter
to 0 for normal behavior. parameter
in the range (0,0.1] specifies the desired accuracy.
A negative parameter
causes round(parameter)
iterations to be run.
This method supports images with a dimensionality between 1 and 3.
"PROJ"
: The PROJ method computes the shift in each dimension separately, applying the ITER method on the projection of the image onto each axis. It is fast and fairly accurate for high SNR. Should not be used for low SNR.parameter
is passed unchanged to the ITER method. This method supports images with any number of dimensions.
in1
and in2
must be scalar images with the same dimensionality and sizes.
dip::FloatArray
dip::FourierMellinMatch2D (dip::Image const& in1,
dip::Image const& in2,
dip::Image& out,
dip::String const& interpolationMethod = S::LINEAR,
dip::String const& correlationMethod = S::PHASE)
Finds the scaling, translation and rotation between two 2D images using the Fourier Mellin transform
The output array represents a 2x3 affine transform matrix (in columnmajor order) that can be used as input
to dip::AffineTransform
to transform in2
to match in1
. out
is the in2
image transformed in this way.
The two images must be realvalued, scalar and twodimensional, and have the same sizes.
interpolationMethod
has a restricted set of options: "linear"
, "3cubic"
, or "nearest"
.
See Interpolation methods for their definition. If in
is binary, interpolationMethod
will be
ignored, nearest neighbor interpolation will be used.
correlationMethod
determines the normalization applied in the cross correlation. It can be
"don't normalize"
, "normalize"
or "phase"
. See dip::CrossCorrelationFT
for an explanation of
these flags.
Example:
dip::Image in1, in2, out; // assign to in1 and in2 two similar images auto T = FourierMellinMatch2D( in1, in2, out ); // out is a transformed version of in2, as per dip::AffineTransform( in2, out, T );
void
dip::StructureTensor (dip::Image const& in,
dip::Image const& mask,
dip::Image& out,
dip::FloatArray const& gradientSigmas = {1.0},
dip::FloatArray const& tensorSigmas = {5.0},
dip::String const& method = S::BEST,
dip::StringArray const& boundaryCondition = {},
dip::dfloat truncation = 3)
Computes the structure tensor.
The structure tensor is a tensor image that contains, at each pixel, information about the local image structure. The eigenvalues of the structure tensor are larger when there are stronger gradients locally. That is, if all eigenvalues are small, the image is locally uniform. If one eigenvalue is large, then there is a unique line or edge orientation (in 2D), or a planelike edge or structure (in 3D). In 3D, if two eigenvalues are large then there is a linelike structure. The associated eigenvalues indicate the orientation of this structure. See the literature references below for more information.
in
must be a scalar, realvalued image. out
will be a symmetric NxN tensor image, where N is the number
of dimensions in in
. Out is computed by:
dip::Image g = dip::Gradient( in, gradientSigmas ); dip::Image out = dip::Gauss( g * dip::Transpose( g ), tensorSigmas );
If mask
is given (not a raw image), then it is interpreted as confidence weights for the input pixels. It
should have values between 0 and 1 (or be binary). Normalized differential convolution is used to compute the
gradients. This method applies Gaussian gradients taking missing and uncertain values into account. See
dip::NormalizedDifferentialConvolution
.
See dip::Gauss
for the meaning of the parameters method
, boundaryCondition
and truncation
.
The functions dip::StructureTensorAnalysis2D
and dip::StructureTensorAnalysis3D
can be used with the output
of this function to obtain useful image parameters.
void
dip::StructureTensorAnalysis2D (dip::Image const& in,
dip::Image* l1 = nullptr,
dip::Image* l2 = nullptr,
dip::Image* orientation = nullptr,
dip::Image* energy = nullptr,
dip::Image* anisotropy1 = nullptr,
dip::Image* anisotropy2 = nullptr,
dip::Image* curvature = nullptr)
Computes useful image parameters from the 2D structure tensor.
in
must be a 2D, symmetric 2x2 tensor image obtained from dip::StructureTensor
. This function takes a pointer
to output images, instead of taking them by reference. Set pointers to nullptr
if you do not want the given
output computed. Use this function as follows:
dip::Image st = dip::StructureTensor( img ); dip::Image energy, orientation; dip::StructureTensorAnalysis2D( st, nullptr, nullptr, &orientation, &energy );
(note how the last two parameters were not given, they default to nullptr
). The code above computes both the
orientation and energy values of the structure tensor.
The output images will be reallocated to be the same size as the input image. They will be scalar and of a floatingpoint type.
The output images are defined as follows:
Image  Description 

l1 
The largest eigenvalue. 
l2 
The smallest eigenvalue. 
orientation 
Orientation. Lies in the interval (π/2, π/2). 
energy 
Sum of the two eigenvalues l1 and l2 . 
anisotropy1 
Measure for local anisotropy: ( l1  l2 ) / ( l1 + l2 ) . 
anisotropy2 
Measure for local anisotropy: 1  l2 / l1 , where l1 > 0 . 
curvature 
Magnitude of the curvature (1/bending radius). 
Curvature sign cannot be computed because the structure tensor does not distinguish the direction of the gradient.
Note that l1
and l2
will both reference data within the same data segment, and therefore will likely not
have normal strides.
For a 3D structure tensor analysis, see the function dip::StructureTensorAnalysis3D
.
A different interface to this function is available in dip::StructureTensorAnalysis
.
Note that eigenvalues and eigenvectors can also be computed using dip::Eigenvalues
and dip::EigenDecomposition
.
void
dip::StructureTensorAnalysis3D (dip::Image const& in,
dip::Image* l1 = nullptr,
dip::Image* phi1 = nullptr,
dip::Image* theta1 = nullptr,
dip::Image* l2 = nullptr,
dip::Image* phi2 = nullptr,
dip::Image* theta2 = nullptr,
dip::Image* l3 = nullptr,
dip::Image* phi3 = nullptr,
dip::Image* theta3 = nullptr,
dip::Image* energy = nullptr,
dip::Image* cylindrical = nullptr,
dip::Image* planar = nullptr)
Computes useful image parameters from the 3D structure tensor.
in
must be a 3D, symmetric 3x3 tensor image obtained from dip::StructureTensor
. This function takes a pointer
to output images, instead of taking them by reference. Set pointers to nullptr
if you do not want the given
output computed. Use this function as follows:
dip::Image st = dip::StructureTensor( img ); dip::Image energy; dip::StructureTensorAnalysis3D( st, nullptr, nullptr, nullptr, nullptr, nullptr, nullptr, nullptr, nullptr, nullptr, &energy );
(note how the last two parameters were not given, they default to nullptr
). The code above computes the
energy value of the structure tensor.
The output images will be reallocated to be the same size as the input image. They will be scalar and of a floatingpoint type.
The output images are defined as follows:
Image  Description 

l1 
The largest eigenvalue. 
phi1 
First component of the orientation of the eigenvector for l1 . 
theta1 
Second component of the orientation of the eigenvector for l1 . 
l2 
The middle eigenvalue. 
phi2 
First component of the orientation of the eigenvector for l2 . 
theta2 
Second component of the orientation of the eigenvector for l2 . 
l3 
The smallest eigenvalue. 
phi3 
First component of the orientation of the eigenvector for l3 . 
theta3 
Second component of the orientation of the eigenvector for l3 . 
energy 
Sum of the three eigenvalues l1 , l2 and l3 . 
cylindrical 
Measure for local anisotropy: ( l2  l3 ) / ( l2 + l3 ) . 
planar 
Measure for local anisotropy: ( l1  l2 ) / ( l1 + l2 ) . 
For a 2D structure tensor analysis, see the function dip::StructureTensorAnalysis2D
.
A different interface to this function is available in dip::StructureTensorAnalysis
.
Note that eigenvalues and eigenvectors can also be computed using dip::Eigenvalues
and dip::EigenDecomposition
.
void
dip::StructureTensorAnalysis (dip::Image const& in,
dip::ImageRefArray& out,
dip::StringArray const& outputs)
Interface to dip::StructureTensorAnalysis2D
and dip::StructureTensorAnalysis3D
.
in
is as in dip::StructureTensorAnalysis2D
or dip::StructureTensorAnalysis3D
. That is, either a 2D
symmetric 2x2 tensor image or a 3D symmetric 3x3 tensor image, realvalued, as obtained from dip::StructureTensor
.
out
is an array with references to output images. It should have exactly as many elements as outputs
.
outputs
is an array with one or more of the following strings, indicating which outputs are needed:
 For 2D inputs:
"l1"
,"l2"
,"orientation"
,"energy"
,"anisotropy1"
,"anisotropy2"
,"curvature"
.  For 3D inputs:
"l1"
,"phi1"
,"theta1"
,"l2"
,"phi2"
,"theta2"
,"l3"
,"phi3"
,"theta3"
,"energy"
,"cylindrical"
,"planar"
.
The order of the strings in outputs
indicates the order they will be written to the out
array.
See the functions dip::StructureTensorAnalysis2D
and dip::StructureTensorAnalysis3D
for more information
on these outputs.
dip::Distribution
dip::StructureAnalysis (dip::Image const& in,
dip::Image const& mask,
std::vector<dfloat> const& scales = {},
dip::String const& feature = "energy",
dip::FloatArray const& gradientSigmas = {1.0},
dip::String const& method = S::BEST,
dip::StringArray const& boundaryCondition = {},
dip::dfloat truncation = 3)
Analyzes the local structure of the image at multiple scales.
Computes the structure tensor, smoothed with each of the values in scales
multiplied by gradientSigmas
,
determines the feature feature
from it, and averages across the image for each scale. This leads to a series
of data points showing how the selected feature changes across scales. scales
and gradientSigmas
are
given in pixels, the image’s pixel size is not taken into account.
The reason that each element of scales
is multiplied by gradientSigmas
is to allow analysis of
nonisotropic images. The gradientSigmas
corrects for the nonisotropy, allowing scales
to be a scalar
value for each scale.
scales
defaults to a series of 10 values geometrically spaced by sqrt(2)
, and starting at 1.0. Units are
pixels.
feature
can be any of the strings allowed by dip::StructureTensorAnalysis
, but "energy"
, "anisotropy1"
,
and "anisotropy2"
(in 2D), and "energy"
, "cylindrical"
, and "planar"
(in 2D) make the most sense.
See dip::Gauss
for the meaning of the parameters method
, boundaryCondition
and truncation
.
in
must be scalar and realvalued, and either 2D or 3D. mask
must be of the same size, and limits the
region over which the structure tensor feature is averaged.
See dip::StructureTensor
for more information about the structure tensor.
void
dip::MonogenicSignal (dip::Image const& in,
dip::Image& out,
dip::FloatArray const& wavelengths = {3.0,24.0},
dip::dfloat bandwidth = 0.41,
dip::String const& inRepresentation = S::SPATIAL,
dip::String const& outRepresentation = S::SPATIAL)
Computes the monogenic signal, a multidimensional generalization of the analytic signal.
The monogenic signal of an ndimensional image has n+1 components. The first component has been filtered with the even (symmetric) filter component, and the other n components have been filtered with the odd (antisymmetric) filter components, where each of those filter components is antisymmetric along a different dimension.
This function splits the frequency spectrum of the image into wavelengths.size()
components. The radial
component of a logGabor filter bank will be used. These radial frequency filters are defined by
wavelengths
(in pixels) and bandwidth
. See dip::LogGaborFilterBank
for details.
The filters are always applied 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.
in
must be scalar and realvalued if given in the spatial domain. If in
is in the frequency domain, it is
expected to be conjugate symmetric, and thus have a realvalued inverse transform (the imaginary part of the
inverse transform will be discarded).
Out will be a tensor image with wavelengths.size()
tensor columns and n+1 tensor rows. The data type will be
singleprecision float for spatialdomain output, or singleprecision complex for frequencydomain output.
void
dip::MonogenicSignalAnalysis (dip::Image const& in,
dip::ImageRefArray& out,
dip::StringArray const& outputs,
dip::dfloat noiseThreshold = 0.2,
dip::dfloat frequencySpreadThreshold = 0.5,
dip::dfloat sigmoidParameter = 10,
dip::dfloat deviationGain = 1.5,
dip::String const& polarity = S::BOTH)
Computes useful image parameters from the monogenic signal.
in
is a tensor image produced by dip::MonogenicSignal
, for which in.TensorRows() == in.Dimensionality() + 1
,
and in.TensorColumns() > 1
. In the case of a single column (a single wavelength was given), the congruency output
is not possible. in
is realvalued (singleprecision float) and in the spatial domain.
out
is an array with references to output images. It should have exactly as many elements as outputs
.
outputs
is an array with one or more of the following strings, indicating which outputs are needed:
Image  Description 

"congruency" 
Phase congruency, a measure of line/edge significance 
"orientation" 
Line/edge orientation (in the range [π/2, π/2]), 2D only 
"phase" 
Phase angle (π/2 is a white line, 0 is an edge, π/2 is a black line) 
"energy" 
Raw phase congruency energy 
"symmetry" 
Phase symmetry, a contrast invariant measure of symmetry 
"symenergy" 
Raw phase symmetry energy 
The order of the strings in outputs
indicates the order they will be written to the out
array.
The output images will be reallocated to be the same size as the input image. They will be scalar and of a floatingpoint type.
noiseThreshold
indicates the noise threshold to use on the energy images when computing phase congruency
or symmetry.
Two different algorithms for phase congruency are implemented:

Kovesi’s method works for images in 2D only, and requires several scales (>2) to be present in the monogenic signal. This method will use the following input arguments:
frequencySpreadThreshold
is a threshold that avoids high phase congruency values if the frequency spread is not large enough.sigmoidParameter
is the parameter to the sigmoid function that weighs congruency with the frequency spread.deviationGain
determines how much the calculated phase deviation should be magnified. Larger values create a sharper response to edges, but also reduces this response’s magnitude. Sensible values are in the range [1,2].

Felsberg’s method works in any number of dimensions, and requires exactly two scales to be present in the monogenic signal. This method ignores the three input arguments specified above.
Phase symmetry can be computed with any number of dimensions and any number of scales. The polarity
input
argument will be used. It is the string "white"
, "black"
or "both"
, indicating which symmetry features
to find.
The other outputs are computed as part of the computation of phase congruency or phase symmetry. These outputs can be requested without requesting the two main parameters, but might require some of the input arguments discussed above to be set.
void
dip::OrientationSpace (dip::Image const& in,
dip::Image& out,
dip::uint order = 8,
dip::dfloat radCenter = 0.1,
dip::dfloat radSigma = 0.8,
dip::uint orientations = 0)
Creates an orientation space for a 2D image
out
is a 3D image, where the 3^{rd} dimension represents orientation. Structures (lines and edges) in
the image are separated into different orientation planes. This transformation separates, for example, intersecting
rings (such as the logo of the Olympic Games) into nontouching objects. This is accomplished using an orientationselective
quadrature filter. Because it’s a quadrature filter, the output is complexvalued, with the real component responding
to lines and the imaginary component responding to edges.
order
determines the angular selectivity. Increasing this number makes the filter longer, and therefore it responds
to a narrower range of orientations. radCenter
determines the overall scaling of the filter. Increasing this
value makes the filter smaller. The radSigma
parameter controls the bandwidth of the filter. Reducing this value
makes the filter more specific to a frequency range (the scale selected by radCenter
).
By default, the orientation space is sampled at order * 2 + 1
different orientations (this is the number
of filters applied, and the number of samples along the 3^{rd} dimension). orientation
, when positive,
determines the number of orientation used. This parameter can be used to oversample the orientation space.
dip::Distribution
dip::PairCorrelation (dip::Image const& object,
dip::Image const& mask,
dip::Random& random,
dip::uint probes = 1000000,
dip::uint length = 100,
dip::String const& sampling = S::RANDOM,
dip::StringSet const& options = {})
Estimates the pair correlation function of the different phases in object
.
If object
is a binary image, the image is a regarded as a twophase image.
In case object
is of an unsigned integer type, the image is regarded as a labeled image,
with each integer value encoding a phase (note that 0 is included here).
Optionally a mask
image can be provided to select which pixels in object
should be used to estimate the pair correlation.
probes
specifies how many random point pairs should be drawn to estimate the function.
length
specifies the maximum correlation length in pixels.
The correlation function can be estimated using a random sampling method (sampling
is
"random"
), or a grid sampling method ("grid"
). For grid sampling, probes
can be 0,
in which case all possible pairs along all image axes are used. Otherwise, the grid covers
a subset of points, with all possible pairs along all image axes; the actual number of pairs
is approximate. In general, the grid method needs a lot more probes to be precise, but
it is faster for a given number of probes because it uses sequential memory access.
options
can contain one or more of the following strings:
"covariance"
: Compute covariance instead of correlation."normalize volume"
: Normalizes the distribution using the volume fraction"normalize volume^2"
: Normalizes the distribution using the square of the volume fraction
Without any of the option strings given, the output has N
values, with N
the number of phases
in the image. Element i
gives the probability of two points at a given distance both hitting
phase number i
. When normalized using the "normalize volume"
option, the probability given is
that of the second point hitting phase i
, given that the first point hits phase i
. This normalization
is computed by dividing all values (distances) for element i
by the volume fraction of phase i
(which is given by the estimated probability at distance equal to 0). When normalized using the
"normalize volume^2"
option, the values given are expected to tends towards 1, as the square of the
volume fraction of phase i
is the probability of two random points hitting phase i
.
These values are read from the output distribution
by distribution[distance].Y(i)
, and the
corresponding distance in physical units is given by distribution[distance].X()
. If object
does not
have a pixel size, its pixel sizes are not isotropic, or have no physical units, then the distances given
are in pixels. distance
runs from 0 to length
(inclusive).
With the "covariance"
option, the output distribution has N
×N
values. The element at
(i,j)
gives the probability for two points at a given distance to land on phases i
and j
respectively. The value for (i,j)
will be identical to that for (j,i)
. The values at the diagonal
will correspond to the pair correlation as described above. The rest of the elements are for the
cases where the two points fall in different phases. Normalization breaks the symmetry, since each
column i
is normalized by the volume fraction for phase i
. That is, (i,j)
will give the probability
that the second point, at a given distance, will hit phase j
given that the first point hits phase i
.
dip::Distribution
dip::ProbabilisticPairCorrelation (dip::Image const& phases,
dip::Image const& mask,
dip::Random& random,
dip::uint probes = 1000000,
dip::uint length = 100,
dip::String const& sampling = S::RANDOM,
dip::StringSet const& options = {})
Estimates the probabilistic pair correlation function of the different phases in phases
.
phases
is a realvalued image of a floatingpoint type, with one or more channels (tensor elements).
Each channel represents the probability per pixel of one of the phases.
The function assumes, but does not check, that these values are with the [0,1] range, and add up to
1 or less.
All other parameters are as in dip::PairCorrelation
.
dip::Distribution
dip::Semivariogram (dip::Image const& in,
dip::Image const& mask,
dip::Random& random,
dip::uint probes = 1000000,
dip::uint length = 100,
dip::String const& sampling = S::RANDOM)
Estimates the expected value of half the square difference between field values at a distance d
.
The input image in
must be scalar and realvalued. An isotropic semivariogram is computed as a function
of the distance, which is varied from 0 to length
(in pixels). Therefore, the field in in
is assumed isotropic
and stationary.
By definition, the semivariogram is zero at the origin (at zero distance). If there is no spatial correlation, the semivariogram is constant for nonzero distance. The magnitude of the jump at the origin is called nugget. The value of the semivariogram in the limit of the distance to infinity is the variance of the field; this value is called the sill. The distance at which the sill is reached (or, more robustly, the distance at which 0.95 times the sill is reached) is called range, and provides a meaningful description of the field.
Optionally a mask
image can be provided to select which pixels in in
should be used to estimate the
semivariogram.
probes
specifies how many random point pairs should be drawn to estimate the semivariogram.
length
specifies the maximum pair distance in pixels.
The semivariogram can be estimated using a random sampling method (sampling
is "random"
), or a grid sampling
method ("grid"
). For grid sampling, probes
can be 0, in which case all possible pairs along all image
axes are used. Otherwise, the grid covers a subset of points, with all possible pairs along all image axes;
the actual number of pairs is approximate. In general, the grid method needs a lot more probes to be precise,
but it is faster for a given number of probes because it uses sequential memory access.
The output dip::Distribution
has x samples given in physical units. However, if in
does not
have a pixel size, its pixel sizes are not isotropic, or have no physical units, then the distances given
are in pixels.
dip::Distribution
dip::ChordLength (dip::Image const& object,
dip::Image const& mask,
dip::Random& random,
dip::uint probes = 100000,
dip::uint length = 100,
dip::String const& sampling = S::RANDOM)
Estimates the chord length distribution of the different phases in object
.
If object
is a binary image, the image is a regarded as a twophase image.
In case object
is of an unsigned integer type, the image is regarded as a labeled image,
with each integer value encoding a phase (note that 0 is included here).
Optionally a mask
image can be provided to select which pixels in object
should be used for the estimation.
probes
specifies how many random lines should be drawn to estimate the distribution.
length
specifies the maximum chord length in pixels.
The chord length distribution can be estimated using a random sampling method (sampling
is
"random"
), or a grid sampling method ("grid"
). For grid sampling, probes
can be 0,
in which case all lines along all image axes are used. Otherwise, the grid covers
a subset of lines; the actual number of lines is approximate. In general, the grid method
needs a lot more probes to be precise, but it is faster for a given number of probes because
it uses sequential memory access.
The output has N
values, with N
the number of phases in the image. Element i
gives the
probability of a chord of length len
pixels within phase number i
. The chord length is the
length of an intersection of a line with the phase.
These values are read from the output distribution
by distribution[len].Y(i)
, and the
corresponding length in physical units is given by distribution[len].X()
. If object
does not
have a pixel size, its pixel sizes are not isotropic, or have no physical units, then the distances given
are in pixels. len
runs from 0 to length
(inclusive).
dip::Distribution
dip::DistanceDistribution (dip::Image const& object,
dip::Image const& region,
dip::uint length = 100)
Computes the distribution of distances to the background of region
for the different phases in object
.
If object
is a binary image, the image is a regarded as a twophase image.
In case object
is of an unsigned integer type, the image is regarded as a labeled image,
with each integer value encoding a phase (note that 0 is included here).
The region
image is a binary image, the distances to its background will be computed.
If region
is a labeled image, the distances to the zerovalued pixels will be computed.
These distances will take the pixels sizes of region
into account. If region
does not
have pixels sizes, those of object
will be used instead.
The output dip::Distribution
has x value increments given by the pixel size. That is,
the distribution has one sample per pixel. length
is the length of the distribution, and
thus limits the distances taken into account.
dip::Distribution
dip::Granulometry (dip::Image const& in,
dip::Image const& mask,
std::vector<dfloat> const& scales = {},
dip::String const& type = "isotropic",
dip::String const& polarity = S::OPENING,
dip::StringSet const& options = {})
Computes the granulometric function for an image
The granulometry yields a volumeweighted, greyvalue–weighted, cumulative distribution of object sizes. It can be used to obtain a size distribution of the image without attempting to segment and separate individual objects. It is computed by a series of openings or closings at different scales. The result at each scale is integrated (summed). The obtained series of values is scaled such that the value for scale 0 is 0, and for scale infinity is 1.
The derivative of the cumulative distribution is a volumeweighted and greyvalue–weighted distribution of
object sizes. See dip::Distribution::Differentiate
.
Greyvalue–weighted means that objects with a larger greyvalue contrast will be weighted more heavily. Ensuring a uniform greyvalue contrast to prevent this characteristic from affecting the estimated size distribution.
Volumeweighted means that objects are weighted by their volume (area in 2D). By dividing the distribution (note: not the cumulative distribution) by the volume corresponding to each scale, it is possible to approximate a countbased distribution.
This function implements various granulometries, to be specified through the parameters type
, polarity
and options
. The following type
values specify the shapes of the structuring element (SE), which determines
the measurement type:
"isotropic"
: An isotropic SE leads to a size distribution dictated by the width of objects."length"
: A line SE leads to a size distribution dictated by the length of objects. We use (constrained) path openings or closings (see Luengo, 2010).
The polarity
flag determines whether it is white objects on a black background ("opening"
) or black objects
on a white background ("closing"
) that are being analyzed.
The options
parameter can contain a set of flags that modify how the operations are applied. The allowed flags
differ depending on the type
flag.

For
"isotropic"
granulometries:"reconstruction"
: uses openings or closings by reconstruction instead of structural openings or closings. This leads to objects not being broken up in the same way. Objects need to be clearly separated spatially for this to work."shifted"
: uses subpixel shifted isotropic structuring elements. This allows a finer sampling of the scale axis (see Luengo et al., 2007). Ignored for images with more than 3 dimensions."interpolate"
: interpolates by a factor up to 8x for smaller scales, attempting to avoid SE diameters smaller than 8. This improves precision of the result for small scales (see Luengo et al., 2007)."subsample"
: subsamples for larger scales, such that the largest SE diameter is 64. This speeds up computation, at the expense of precision.

For
"length"
granulometries:"unconstrained"
: by default, we use constrained path openings or closings, which improves the precision of the measurement, but is a little bit more expensive (see Luengo, 2010). This option causes the use of normal path openings or closings."robust"
: applies path openings or closings in such a way that they are less sensitive to noise.
scales
defaults to a series of 12 values geometrically spaced by sqrt(2)
, and starting at sqrt(2)
. scales
are in pixels, the image’s pixel size is not taken into account.
in
must be scalar and realvalued. mask
must have the same sizes, and limits the region
in which objects are measured.
dip::dfloat
dip::FractalDimension (dip::Image const& in,
dip::dfloat eta = 0.5)
Estimates the fractal dimension of the binary image in
the sliding box method.
The sliding box method is an enhancement of the classical box counting method that counts many more boxes at each scale, and therefore is more precise. By sliding the box one pixel at a time, it is also not affected by partial boxes (i.e. the boxes at the right and bottom edge of the image that do not fit within the image domain). The counts are computed in an efficient manner, which makes it similar in complexity to counting only the set of tessellating boxes.
The smallest scale used is a box size of 1, and the largest scale is at most half the smallest image size (i.e.
min(width,height)/2
. In between, scales grow exponentially with a factor 1+eta
. Thus, if eta
is 1,
then each scale uses boxes double the size of the previous scale, and if eta
is smaller then the steps are
smaller and more scales are generated.
The image in
must be scalar and binary, and typically is applied to an edge map of objects.
void
dip::swap (Sample& a,
Sample& b)
#include "diplib/distribution.h"
Swaps two samples, copying the data from other
to *this
, and that from *this
to other
.
Both must have the same number of values.
std::ostream&
dip::operator<< (std::ostream& os,
dip::Distribution const& distribution)
#include "diplib/distribution.h"
Writes the distribution to a stream