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

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

• Reference

## 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 central normalized moments, in `N` dimensions.

## Enums

enum class dip::Option::Periodicity: int{
Select if the operation is periodic or not. Used in `dip::GaussianMixtureModel`.

## Functions

template<typename T>
) -> 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>>
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.
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(
Computes the Bessel function J of the order 0.
auto dip::BesselJ1(
Computes the Bessel function J of the order 1.
auto dip::BesselJN(
Computes the Bessel function J of the order `n`.
auto dip::BesselY0(
Computes the Bessel function Y of the order 0.
auto dip::BesselY1(
Computes the Bessel function Y of the order 1.
auto dip::BesselYN(
Computes the Bessel function Y of the order `n`.
auto dip::Sinc(
Computes the sinc function.
auto dip::Phi(
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(
Computes phi, the integral of the PDF of a Normal distribution with standard deviation `s` and mean `m` from minus infinity to `x`.
dip::uint n, dip::dfloat r) -> dip::dfloat constexpr
Computes the surface area of an `n`-dimensional hypersphere with radius `r`.
dip::uint n, dip::dfloat r) -> dip::dfloat constexpr
Computes the volume of an `n`-dimensional hypersphere with radius `r`.
dip::uint n, ConstSampleIterator<dip::dfloat> input, SampleIterator<dip::dfloat> lambdas, SampleIterator<dip::dfloat> vectors = nullptr)
Finds the eigenvalues and eigenvectors of a symmetric, real-valued matrix.
ConstSampleIterator<dip::dfloat> input, SampleIterator<dip::dfloat> lambdas, SampleIterator<dip::dfloat> vectors = nullptr)
Finds the eigenvalues and eigenvectors of a 2x2 symmetric, real-valued matrix.
ConstSampleIterator<dip::dfloat> input, SampleIterator<dip::dfloat> lambdas, SampleIterator<dip::dfloat> vectors = nullptr)
Finds the eigenvalues and eigenvectors of a 3x3 symmetric, real-valued matrix.
Finds the largest eigenvector of a symmetric, real-valued matrix.
Finds the smallest eigenvector of a symmetric, real-valued matrix.
dip::uint n, ConstSampleIterator<dip::dfloat> input, SampleIterator<dip::dfloat> lambdas, SampleIterator<dip::dfloat> vectors = nullptr)
Finds the eigenvalues and eigenvectors of a symmetric, real-valued matrix, where only the unique values are given.
dip::uint n, ConstSampleIterator<dip::dfloat> input, SampleIterator<dip::dcomplex> lambdas, SampleIterator<dip::dcomplex> vectors = nullptr)
Finds the eigenvalues and eigenvectors of a square, real-valued matrix.
dip::uint n, ConstSampleIterator<dip::dcomplex> input, SampleIterator<dip::dcomplex> lambdas, SampleIterator<dip::dcomplex> vectors = nullptr)
Finds the eigenvalues and eigenvectors of a square, complex-valued matrix.
template<typename T>
auto dip::Sum(dip::uint n, ConstSampleIterator<T> input) -> T
Computes the sum of the values of a vector.
template<typename T>
auto dip::SumAbsSquare(dip::uint n, ConstSampleIterator<T> input) -> FloatType<T>
Computes the sum of the square of the values of a vector.
template<typename T>
auto dip::Product(dip::uint n, ConstSampleIterator<T> input) -> T
Computes the product of the values of a vector.
template<typename T>
auto dip::Norm(dip::uint n, ConstSampleIterator<T> input) -> FloatType<T>
Computes the norm of a vector.
template<typename T>
auto dip::SquareNorm(dip::uint n, ConstSampleIterator<T> input) -> FloatType<T>
Computes the square norm of a vector.
auto dip::Determinant(
Computes the determinant of a square, real-valued matrix.
auto dip::Determinant(
Computes the determinant of a square, complex-valued matrix.
template<typename T>
dip::uint n, ConstSampleIterator<T> input) -> T
Computes the determinant of a diagonal matrix.
template<typename T>
auto dip::Trace(dip::uint n, ConstSampleIterator<T> input) -> T
Computes the trace of a square matrix.
template<typename T>
auto dip::TraceDiagonal(dip::uint n, ConstSampleIterator<T> input) -> T
Computes the trace of a diagonal matrix.
dip::uint m, dip::uint n, ConstSampleIterator<dip::dfloat> input, SampleIterator<dip::dfloat> output, SampleIterator<dip::dfloat> U = nullptr, SampleIterator<dip::dfloat> V = nullptr)
Computes the “thin” singular value decomposition of a real-valued matrix
dip::uint m, dip::uint n, ConstSampleIterator<dip::dcomplex> input, SampleIterator<dip::dcomplex> output, SampleIterator<dip::dcomplex> U = nullptr, SampleIterator<dip::dcomplex> V = nullptr)
Computes the “thin” singular value decomposition of a complex-valued matrix
void dip::Inverse(
Computes the inverse of a square, real-valued matrix.
void dip::Inverse(
Computes the inverse of a square, complex-valued matrix.
void dip::PseudoInverse(dip::uint m, dip::uint n, ConstSampleIterator<dip::dfloat> input, SampleIterator<dip::dfloat> 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, ConstSampleIterator<dip::dcomplex> input, SampleIterator<dip::dcomplex> 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(
Computes the rank of a real-valued matrix.
auto dip::Rank(
Computes the rank of a complex-valued matrix.
void dip::Solve(
Solves a system of real-valued equations, using the Jacobi SVD decomposition.
ConstSampleIterator<dip::dfloat> data, SampleIterator<dip::dfloat> responsibilities, dip::uint size, dip::uint numberOfGaussians, dip::uint maxIter = 20, Option::Periodicity periodicity = Option::Periodicity::NOT_PERIODIC) -> std::vector<GaussianParameters>
Determines the parameters for a Gaussian Mixture Model.

## Operators

auto dip::operator+(
Combine two accumulators
auto dip::operator+(dip::VarianceAccumulator lhs, dip::VarianceAccumulator const& rhs) -> dip::VarianceAccumulator
Combine two accumulators
auto dip::operator+(
Combine two accumulators
auto dip::operator+(
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: int

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, ConstSampleIterator<dip::dfloat> input, SampleIterator<dip::dfloat> lambdas, SampleIterator<dip::dfloat> 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, ConstSampleIterator<dip::dfloat> input, SampleIterator<dip::dfloat> 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, ConstSampleIterator<dip::dfloat> input, SampleIterator<dip::dfloat> 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, ConstSampleIterator<dip::dfloat> input, SampleIterator<dip::dfloat> lambdas, SampleIterator<dip::dfloat> 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, ConstSampleIterator<dip::dfloat> input, SampleIterator<dip::dcomplex> lambdas, SampleIterator<dip::dcomplex> 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, ConstSampleIterator<dip::dcomplex> input, SampleIterator<dip::dcomplex> lambdas, SampleIterator<dip::dcomplex> 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, ConstSampleIterator<T> input)

Computes the sum of the values of a vector.

`input` is a pointer to `n` values.

### template<typename T> FloatType<T> dip::SumAbsSquare(dip::uint n, ConstSampleIterator<T> 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, ConstSampleIterator<T> input)

Computes the product of the values of a vector.

`input` is a pointer to `n` values.

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

Computes the norm of a vector.

`input` is a pointer to `n` values.

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

Computes the square norm of a vector.

`input` is a pointer to `n` values.

### dip::dfloat dip::Determinant(dip::uint n, ConstSampleIterator<dip::dfloat> 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, ConstSampleIterator<dip::dcomplex> 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, ConstSampleIterator<T> 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, ConstSampleIterator<T> 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, ConstSampleIterator<T> 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, ConstSampleIterator<dip::dfloat> input, SampleIterator<dip::dfloat> output, SampleIterator<dip::dfloat> U = nullptr, SampleIterator<dip::dfloat> 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, ConstSampleIterator<dip::dcomplex> input, SampleIterator<dip::dcomplex> output, SampleIterator<dip::dcomplex> U = nullptr, SampleIterator<dip::dcomplex> 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, ConstSampleIterator<dip::dfloat> input, SampleIterator<dip::dfloat> 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, ConstSampleIterator<dip::dcomplex> input, SampleIterator<dip::dcomplex> 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, ConstSampleIterator<dip::dfloat> input, SampleIterator<dip::dfloat> 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, ConstSampleIterator<dip::dcomplex> input, SampleIterator<dip::dcomplex> 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, ConstSampleIterator<dip::dfloat> 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, ConstSampleIterator<dip::dcomplex> 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, ConstSampleIterator<dip::dfloat> A, ConstSampleIterator<dip::dfloat> b, SampleIterator<dip::dfloat> output)

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

Solves , where `A` is an `m`x`n` 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(ConstSampleIterator<dip::dfloat> data, SampleIterator<dip::dfloat> responsibilities, dip::uint size, dip::uint numberOfGaussians, dip::uint maxIter = 20, 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 and the sigma and 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.