The library infrastructure » Numeric algorithms and constants module
#include "diplib.h"

Functions and constants to be used in numeric computation, unrelated to images.

Contents

Classes

class dip::ThinPlateSpline
Fits a thin plate spline function to a set of points. Useful for interpolation of scattered points.
struct dip::GaussianParameters
Parameters defining a 1D Gaussian. Returned by dip::GaussianMixtureModel.
class dip::StatisticsAccumulator
StatisticsAccumulator computes population statistics by accumulating the first four central moments.
class dip::VarianceAccumulator
VarianceAccumulator computes mean and standard deviation by accumulating the first two central moments.
class dip::FastVarianceAccumulator
FastVarianceAccumulator computes mean and standard deviation by accumulating the sum of sample values and the sum of the square of sample values.
class dip::CovarianceAccumulator
CovarianceAccumulator computes covariance and correlation of pairs of samples by accumulating the first two central moments and cross-moments.
class dip::DirectionalStatisticsAccumulator
DirectionalStatisticsAccumulator computes directional mean and standard deviation by accumulating a unit vector with the input value as angle.
class dip::MinMaxAccumulator
MinMaxAccumulator computes minimum and maximum values of a sequence of values.
class dip::MomentAccumulator
MomentAccumulator accumulates the zeroth order moment, the first order normalized moments, and the second order normalized central moments, in N dimensions.

Enums

enum class dip::Option::Periodicity: uint8{ NOT_PERIODIC, PERIODIC }
Select if the operation is periodic or not. Used in dip::GaussianMixtureModel.

Functions

template<typename T>
auto dip::maximum_gauss_truncation() -> dip::dfloat constexpr
Maximum meaningful truncation value for a Gaussian. Larger truncation values will lead to differences of more than one machine epsilon between the middle and the ends of the Gaussian. T must be a floating-point type.
auto dip::gcd(dip::uint a, dip::uint b) -> dip::uint constexpr
Compute the greatest common denominator of two positive integers.
template<typename T, <SFINAE> = 0>
auto dip::div_ceil(T lhs, T rhs) -> T constexpr
Integer division, unsigned, return ceil.
template<typename T, <SFINAE> = 0>
auto dip::div_floor(T lhs, T rhs) -> T constexpr
Integer division, unsigned, return floor.
template<typename T, typename <SFINAE> = T>
auto dip::div_round(T lhs, T rhs) -> T constexpr
Integer division, return rounded.
auto dip::modulo(dip::uint value, dip::uint period) -> dip::uint constexpr
Integer modulo, result is always positive, as opposed to % operator.
auto dip::modulo(dip::sint value, dip::sint period) -> dip::sint constexpr
Integer modulo, result is always positive, as opposed to % operator.
template<typename T, typename <SFINAE>>
auto dip::floor_cast(T v) -> dip::sint constexpr
Fast floor operation, without checks, returning a dip::sint.
template<typename T, typename <SFINAE>>
auto dip::ceil_cast(T v) -> dip::sint constexpr
Fast ceil operation, without checks, returning a dip::sint.
template<typename T, typename <SFINAE>>
auto dip::round_cast(T v) -> dip::sint
Fast round operation, without checks, returning a dip::sint.
template<typename T, bool inverse = false, typename <SFINAE>>
auto dip::consistent_round(T v) -> dip::sint constexpr
Consistent rounding, without checks, returning a dip::sint.
template<typename T>
auto dip::abs(T value) -> T constexpr
constexpr version of std::abs. Prefer std::abs outside of constexpr functions.
template<typename T>
auto dip::clamp(T const& v, T const& lo, T const& hi) -> T const& constexpr
Clamps a value between a min and max value (a.k.a. clip, saturate, etc.).
auto dip::pow10(dip::sint power) -> dip::dfloat constexpr
Computes integer powers of 10, assuming the power is relatively small.
auto dip::ApproximatelyEquals(dip::dfloat lhs, dip::dfloat rhs, dip::dfloat tolerance = 1e-6) -> bool constexpr
Approximate floating-point equality: abs(lhs-rhs)/lhs <= tolerance.
auto dip::LengthUnicode(dip::String const& string) -> dip::uint
Counts the length of a (UTF-8 encoded) Unicode string.
auto dip::BesselJ0(dip::dfloat x) -> dip::dfloat
Computes the Bessel function J of the order 0 (with around 7 digits of precision).
auto dip::BesselJ1(dip::dfloat x) -> dip::dfloat
Computes the Bessel function J of the order 1 (with around 7 digits of precision).
auto dip::BesselJN(dip::dfloat x, dip::uint n) -> dip::dfloat
Computes the Bessel function J of the order n (with around 7 digits of precision).
auto dip::BesselY0(dip::dfloat x) -> dip::dfloat
Computes the Bessel function Y of the order 0 (with around 7 digits of precision).
auto dip::BesselY1(dip::dfloat x) -> dip::dfloat
Computes the Bessel function Y of the order 1 (with around 7 digits of precision).
auto dip::BesselYN(dip::dfloat x, dip::uint n) -> dip::dfloat
Computes the Bessel function Y of the order n (with around 7 digits of precision).
auto dip::Sinc(dip::dfloat x) -> dip::dfloat
Computes the sinc function.
auto dip::Phi(dip::dfloat x) -> dip::dfloat
Computes phi, the integral of the PDF of a Normal distribution with unit variance and zero mean from minus infinity to x.
auto dip::Phi(dip::dfloat x, dip::dfloat m, dip::dfloat s) -> dip::dfloat
Computes phi, the integral of the PDF of a Normal distribution with standard deviation s and mean m from minus infinity to x.
auto dip::HypersphereSurface(dip::uint n, dip::dfloat r) -> dip::dfloat constexpr
Computes the surface area of an n-dimensional hypersphere with radius r.
auto dip::HypersphereVolume(dip::uint n, dip::dfloat r) -> dip::dfloat constexpr
Computes the volume of an n-dimensional hypersphere with radius r.
void dip::SymmetricEigenDecomposition(dip::uint n, dip::ConstSampleIterator input, dip::SampleIterator lambdas, dip::SampleIterator vectors = nullptr)
Finds the eigenvalues and eigenvectors of a symmetric, real-valued matrix.
void dip::SymmetricEigenDecomposition2(dip::ConstSampleIterator input, dip::SampleIterator lambdas, dip::SampleIterator vectors = nullptr)
Finds the eigenvalues and eigenvectors of a 2x2 symmetric, real-valued matrix.
void dip::SymmetricEigenDecomposition3(dip::ConstSampleIterator input, dip::SampleIterator lambdas, dip::SampleIterator vectors = nullptr)
Finds the eigenvalues and eigenvectors of a 3x3 symmetric, real-valued matrix.
void dip::LargestEigenvector(dip::uint n, dip::ConstSampleIterator input, dip::SampleIterator vector)
Finds the largest eigenvector of a symmetric, real-valued matrix.
void dip::SmallestEigenvector(dip::uint n, dip::ConstSampleIterator input, dip::SampleIterator vector)
Finds the smallest eigenvector of a symmetric, real-valued matrix.
void dip::SymmetricEigenDecompositionPacked(dip::uint n, dip::ConstSampleIterator input, dip::SampleIterator lambdas, dip::SampleIterator vectors = nullptr)
Finds the eigenvalues and eigenvectors of a symmetric, real-valued matrix, where only the unique values are given.
void dip::EigenDecomposition(dip::uint n, dip::ConstSampleIterator input, dip::SampleIterator lambdas, dip::SampleIterator vectors = nullptr)
Finds the eigenvalues and eigenvectors of a square, real-valued matrix.
void dip::EigenDecomposition(dip::uint n, dip::ConstSampleIterator input, dip::SampleIterator lambdas, dip::SampleIterator vectors = nullptr)
Finds the eigenvalues and eigenvectors of a square, complex-valued matrix.
template<typename T>
auto dip::Sum(dip::uint n, dip::ConstSampleIterator input) -> T
Computes the sum of the values of a vector.
template<typename T>
auto dip::SumAbsSquare(dip::uint n, dip::ConstSampleIterator input) -> dip::FloatType
Computes the sum of the square of the values of a vector.
template<typename T>
auto dip::Product(dip::uint n, dip::ConstSampleIterator input) -> T
Computes the product of the values of a vector.
template<typename T>
auto dip::Norm(dip::uint n, dip::ConstSampleIterator input) -> dip::FloatType
Computes the norm of a vector.
template<typename T>
auto dip::SquareNorm(dip::uint n, dip::ConstSampleIterator input) -> dip::FloatType
Computes the square norm of a vector.
auto dip::Determinant(dip::uint n, dip::ConstSampleIterator input) -> dip::dfloat
Computes the determinant of a square, real-valued matrix.
auto dip::Determinant(dip::uint n, dip::ConstSampleIterator input) -> dip::dcomplex
Computes the determinant of a square, complex-valued matrix.
template<typename T>
auto dip::DeterminantDiagonal(dip::uint n, dip::ConstSampleIterator input) -> T
Computes the determinant of a diagonal matrix.
template<typename T>
auto dip::Trace(dip::uint n, dip::ConstSampleIterator input) -> T
Computes the trace of a square matrix.
template<typename T>
auto dip::TraceDiagonal(dip::uint n, dip::ConstSampleIterator input) -> T
Computes the trace of a diagonal matrix.
void dip::SingularValueDecomposition(dip::uint m, dip::uint n, dip::ConstSampleIterator input, dip::SampleIterator output, dip::SampleIterator U = nullptr, dip::SampleIterator V = nullptr)
Computes the “thin” singular value decomposition of a real-valued matrix
void dip::SingularValueDecomposition(dip::uint m, dip::uint n, dip::ConstSampleIterator input, dip::SampleIterator output, dip::SampleIterator U = nullptr, dip::SampleIterator V = nullptr)
Computes the “thin” singular value decomposition of a complex-valued matrix
void dip::Inverse(dip::uint n, dip::ConstSampleIterator input, dip::SampleIterator output)
Computes the inverse of a square, real-valued matrix.
void dip::Inverse(dip::uint n, dip::ConstSampleIterator input, dip::SampleIterator output)
Computes the inverse of a square, complex-valued matrix.
void dip::PseudoInverse(dip::uint m, dip::uint n, dip::ConstSampleIterator input, dip::SampleIterator output, dip::dfloat tolerance = 1e-7)
Computes the Moore-Penrose pseudo-inverse of a real-valued matrix, using the Jacobi SVD decomposition.
void dip::PseudoInverse(dip::uint m, dip::uint n, dip::ConstSampleIterator input, dip::SampleIterator output, dip::dfloat tolerance = 1e-7)
Computes the Moore-Penrose pseudo-inverse of a complex-valued matrix, using the Jacobi SVD decomposition.
auto dip::Rank(dip::uint m, dip::uint n, dip::ConstSampleIterator input) -> dip::uint
Computes the rank of a real-valued matrix.
auto dip::Rank(dip::uint m, dip::uint n, dip::ConstSampleIterator input) -> dip::uint
Computes the rank of a complex-valued matrix.
void dip::Solve(dip::uint m, dip::uint n, dip::ConstSampleIterator A, dip::ConstSampleIterator b, dip::SampleIterator output)
Solves a system of real-valued equations, using the Jacobi SVD decomposition.
auto dip::GaussianMixtureModel(dip::ConstSampleIterator data, dip::SampleIterator responsibilities, dip::uint size, dip::uint numberOfGaussians, dip::uint maxIter = 20, dip::Option::Periodicity periodicity = Option::Periodicity::NOT_PERIODIC) -> std::vector<GaussianParameters>
Determines the parameters for a Gaussian Mixture Model.

Operators

auto dip::operator+(dip::StatisticsAccumulator lhs, dip::StatisticsAccumulator const& rhs) -> dip::StatisticsAccumulator
Combine two accumulators
auto dip::operator+(dip::VarianceAccumulator lhs, dip::VarianceAccumulator const& rhs) -> dip::VarianceAccumulator
Combine two accumulators
auto dip::operator+(dip::FastVarianceAccumulator lhs, dip::FastVarianceAccumulator const& rhs) -> dip::FastVarianceAccumulator
Combine two accumulators
auto dip::operator+(dip::DirectionalStatisticsAccumulator lhs, dip::DirectionalStatisticsAccumulator const& rhs) -> dip::DirectionalStatisticsAccumulator
Combine two accumulators
auto dip::operator+(dip::MinMaxAccumulator lhs, dip::MinMaxAccumulator const& rhs) -> dip::MinMaxAccumulator
Combine two accumulators
auto dip::operator+(dip::MomentAccumulator lhs, dip::MomentAccumulator const& rhs) -> dip::MomentAccumulator
Combine two accumulators

Variables

dip::dfloat const dip::pi constexpr
The constant π.
dip::dfloat const dip::nan constexpr
A NaN value.
dip::dfloat const dip::infinity constexpr
Infinity.

Class documentation

struct dip::GaussianParameters

Parameters defining a 1D Gaussian. Returned by dip::GaussianMixtureModel.

Variables
dip::dfloat position The location of the origin, in pixels
dip::dfloat amplitude The amplitude (value at the origin)
dip::dfloat sigma The sigma (width)

Enum documentation

enum class dip::Option::Periodicity: uint8

Select if the operation is periodic or not. Used in dip::GaussianMixtureModel.

Enumerators
NOT_PERIODIC = 0 The operation is not periodic
PERIODIC = 1 The operation is periodic, left and right ends of the data are contiguous

Function documentation

template<typename T, bool inverse = false, typename <SFINAE>>
dip::sint dip::consistent_round(T v) constexpr

Consistent rounding, without checks, returning a dip::sint.

This rounding is consistent in that half-way cases are rounded in the same direction for positive and negative values. The inverse template parameter indicates the direction for these cases. By default, it matches std::round for positive values.

void dip::SymmetricEigenDecomposition(dip::uint n, dip::ConstSampleIterator input, dip::SampleIterator lambdas, dip::SampleIterator vectors = nullptr)

Finds the eigenvalues and eigenvectors of a symmetric, real-valued matrix.

input is a pointer to n*n values, in column-major order; only the lower triangle will be used.

lambdas is a pointer to space for n values, which will be written sorted by magnitude, largest to smallest.

vectors is a pointer to space for n*n values and will receive the n eigenvectors. The eigenvectors can be accessed at &vectors[ 0 ], &vectors[ n ], &vectors[ 2*n ], etc. If vectors is nullptr, no eigenvectors are computed.

void dip::LargestEigenvector(dip::uint n, dip::ConstSampleIterator input, dip::SampleIterator vector)

Finds the largest eigenvector of a symmetric, real-valued matrix.

input is a pointer to n*n values, in column-major order; only the lower triangle will be used.

vector is a pointer to space for n values, and will receive the eigenvector corresponding to the largest eigenvalue by magnitude. The full decomposition as in dip::SymmetricEigenDecomposition is computed, but only one eigenvector is written to the output.

void dip::SmallestEigenvector(dip::uint n, dip::ConstSampleIterator input, dip::SampleIterator vector)

Finds the smallest eigenvector of a symmetric, real-valued matrix.

input is a pointer to n*n values, in column-major order; only the lower triangle will be used.

vector is a pointer to space for n values, and will receive the eigenvector corresponding to the smallest eigenvalue by magnitude. The full decomposition as in dip::SymmetricEigenDecomposition is computed, but only one eigenvector is written to the output.

void dip::SymmetricEigenDecompositionPacked(dip::uint n, dip::ConstSampleIterator input, dip::SampleIterator lambdas, dip::SampleIterator vectors = nullptr)

Finds the eigenvalues and eigenvectors of a symmetric, real-valued matrix, where only the unique values are given.

Calls dip::SymmetricEigenDecomposition after copying over the input values to a temporary buffer.

input is a pointer to n*(n+1)/2 values, stored in the same order as symmetric tensors are stored in an image (see dip::Tensor::Shape). That is, fist are the main diagonal elements, then the elements above the diagonal, column-wise. This translates to:

  • 2D: xx, yy, xy
  • 3D: xx, yy, zz, xy, xz, yz
  • 4D: xx, yy, zz, tt, xy, xz, yz, xt, yt, zt
  • etc.

See dip::SymmetricEigenDecomposition for information on lambdas and vectors.

void dip::EigenDecomposition(dip::uint n, dip::ConstSampleIterator input, dip::SampleIterator lambdas, dip::SampleIterator vectors = nullptr)

Finds the eigenvalues and eigenvectors of a square, real-valued matrix.

input is a pointer to n*n values, in column-major order.

lambdas is a pointer to space for n values, sorted by magnitude, largest to smallest

vectors is a pointer to space for n*n values and will receive the n eigenvectors. The eigenvectors can be accessed at &vectors[ 0 ], &vectors[ n ], &vectors[ 2*n ], etc. If vectors is nullptr, no eigenvectors are computed.

void dip::EigenDecomposition(dip::uint n, dip::ConstSampleIterator input, dip::SampleIterator lambdas, dip::SampleIterator vectors = nullptr)

Finds the eigenvalues and eigenvectors of a square, complex-valued matrix.

input is a pointer to n*n values, in column-major order.

lambdas is a pointer to space for n values, sorted by magnitude, largest to smallest

vectors is a pointer to space for n*n values and will receive the n eigenvectors. The eigenvectors can be accessed at &vectors[ 0 ], &vectors[ n ], &vectors[ 2*n ], etc. If vectors is nullptr, no eigenvectors are computed.

template<typename T>
T dip::Sum(dip::uint n, dip::ConstSampleIterator input)

Computes the sum of the values of a vector.

input is a pointer to n values.

template<typename T>
dip::FloatType dip::SumAbsSquare(dip::uint n, dip::ConstSampleIterator input)

Computes the sum of the square of the values of a vector.

input is a pointer to n values.

template<typename T>
T dip::Product(dip::uint n, dip::ConstSampleIterator input)

Computes the product of the values of a vector.

input is a pointer to n values.

template<typename T>
dip::FloatType dip::Norm(dip::uint n, dip::ConstSampleIterator input)

Computes the norm of a vector.

input is a pointer to n values.

template<typename T>
dip::FloatType dip::SquareNorm(dip::uint n, dip::ConstSampleIterator input)

Computes the square norm of a vector.

input is a pointer to n values.

dip::dfloat dip::Determinant(dip::uint n, dip::ConstSampleIterator input)

Computes the determinant of a square, real-valued matrix.

input is a pointer to n*n values, in column-major order.

dip::dcomplex dip::Determinant(dip::uint n, dip::ConstSampleIterator input)

Computes the determinant of a square, complex-valued matrix.

input is a pointer to n*n values, in column-major order.

template<typename T>
T dip::DeterminantDiagonal(dip::uint n, dip::ConstSampleIterator input)

Computes the determinant of a diagonal matrix.

input is a pointer to n values, representing the matrix’s main diagonal.

template<typename T>
T dip::Trace(dip::uint n, dip::ConstSampleIterator input)

Computes the trace of a square matrix.

input is a pointer to n*n values, in column-major order.

template<typename T>
T dip::TraceDiagonal(dip::uint n, dip::ConstSampleIterator input)

Computes the trace of a diagonal matrix.

input is a pointer to n values, representing the matrix’s main diagonal.

void dip::SingularValueDecomposition(dip::uint m, dip::uint n, dip::ConstSampleIterator input, dip::SampleIterator output, dip::SampleIterator U = nullptr, dip::SampleIterator V = nullptr)

Computes the “thin” singular value decomposition of a real-valued matrix

input is a pointer to m*n values, in column-major order.

output is a pointer to p values, where p = std::min( m, n ). It contains the singular values of input, sorted in decreasing order.

U and V are pointers to m*p and n*p values, respectively. The left and right singular vectors will be written to then. If either of them is nullptr, neither is computed, and only output is filled.

SingularValueDecomposition uses the two-sided Jacobi SVD decomposition algorithm. This is efficient for small matrices only.

void dip::SingularValueDecomposition(dip::uint m, dip::uint n, dip::ConstSampleIterator input, dip::SampleIterator output, dip::SampleIterator U = nullptr, dip::SampleIterator V = nullptr)

Computes the “thin” singular value decomposition of a complex-valued matrix

input is a pointer to m*n values, in column-major order.

output is a pointer to p values, where p = std::min( m, n ). It contains the singular values of input, sorted in decreasing order.

U and V are pointers to m*p and n*p values, respectively. The left and right singular vectors will be written to then. If either of them is nullptr, neither is computed, and only output is filled.

SingularValueDecomposition uses the two-sided Jacobi SVD decomposition algorithm. This is efficient for small matrices only.

void dip::Inverse(dip::uint n, dip::ConstSampleIterator input, dip::SampleIterator output)

Computes the inverse of a square, real-valued matrix.

input and output are pointers to n*n values, in column-major order.

The result is undetermined if the matrix is not invertible.

void dip::Inverse(dip::uint n, dip::ConstSampleIterator input, dip::SampleIterator output)

Computes the inverse of a square, complex-valued matrix.

input and output are pointers to n*n values, in column-major order.

The result is undetermined if the matrix is not invertible.

void dip::PseudoInverse(dip::uint m, dip::uint n, dip::ConstSampleIterator input, dip::SampleIterator output, dip::dfloat tolerance = 1e-7)

Computes the Moore-Penrose pseudo-inverse of a real-valued matrix, using the Jacobi SVD decomposition.

input is a pointer to m*n values, in column-major order.

output is a pointer to n*m values, in column-major order.

tolerance is an appropriate tolerance. Singular values smaller than tolerance * max(n,m) * p, with p the largest singular value, will be set to zero in the inverse.

void dip::PseudoInverse(dip::uint m, dip::uint n, dip::ConstSampleIterator input, dip::SampleIterator output, dip::dfloat tolerance = 1e-7)

Computes the Moore-Penrose pseudo-inverse of a complex-valued matrix, using the Jacobi SVD decomposition.

input and output are pointers to m*n values, in column-major order.

output is a pointer to n*m values, in column-major order.

tolerance is an appropriate tolerance. Singular values smaller than tolerance * max(n,m) * p, with p the largest singular value, will be set to zero in the inverse.

dip::uint dip::Rank(dip::uint m, dip::uint n, dip::ConstSampleIterator input)

Computes the rank of a real-valued matrix.

input is a pointer to m*n values, in column-major order.

dip::uint dip::Rank(dip::uint m, dip::uint n, dip::ConstSampleIterator input)

Computes the rank of a complex-valued matrix.

input is a pointer to m*n values, in column-major order.

void dip::Solve(dip::uint m, dip::uint n, dip::ConstSampleIterator A, dip::ConstSampleIterator b, dip::SampleIterator output)

Solves a system of real-valued equations, using the Jacobi SVD decomposition.

Solves \(A x = b\) , where A is an mxn matrix (stored in column-major order), and b is a vector with m values. The unknown x will have n values, and will be written to output.

std::vector<GaussianParameters> dip::GaussianMixtureModel(dip::ConstSampleIterator data, dip::SampleIterator responsibilities, dip::uint size, dip::uint numberOfGaussians, dip::uint maxIter = 20, dip::Option::Periodicity periodicity = Option::Periodicity::NOT_PERIODIC)

Determines the parameters for a Gaussian Mixture Model.

data is an iterator (or pointer) to the first of size samples of a GMM (not random samples drawn from such a distribution, but rather samples of a function representing the distribution). numberOfGaussians Gaussians will be fitted to it using the Expectation Maximization (EM) procedure.

The parameters are initialized deterministically, the means are distributed equally over the domain, the sigma are all set to the distance between means, and the amplitude are set to 1.

responsibilities optionally points to a buffer of size size * numberOfGaussians that will be used internally. If set to nullptr or a default-initialized iterator, a buffer will be allocated internally. Use this parameter when repeatedly calling this function to avoid memory allocations.

maxIter sets how many iterations are run. There is currently no other stopping criterion. periodicity determines if the data is considered periodic or not.

The output is sorted by amplitude, most important component first.