The library infrastructure » Numeric algorithms and constants module

#include "diplib.h"

struct dip::GaussianParameters

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

enum class dip::Option::DecompositionMethod: uint8

Select the algorithm to use with dip::SymmetricEigenDecomposition2 and dip::SymmetricEigenDecomposition3.

enum class dip::Option::Periodicity: uint8

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

template<typename T>

dip::dfloat dip::maximum_gauss_truncation( ) 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.

dip::uint dip::gcd( dip::uint a, dip::uint b) constexpr

Compute the greatest common denominator of two positive integers.

template<typename T, <SFINAE> = 0>

T dip::div_ceil( T lhs, T rhs) constexpr

Integer division, unsigned, return ceil.

template<typename T, <SFINAE> = 0>

T dip::div_floor( T lhs, T rhs) constexpr

Integer division, unsigned, return floor.

template<typename T, typename <SFINAE> = T>

T dip::div_round( T lhs, T rhs) constexpr

Integer division, return rounded.

dip::uint dip::modulo( dip::uint value, dip::uint period) constexpr

Integer modulo, result is always positive, as opposed to % operator.

dip::sint dip::modulo( dip::sint value, dip::sint period) constexpr

Integer modulo, result is always positive, as opposed to % operator.

template<typename T, typename <SFINAE>>

dip::sint dip::floor_cast( T v) constexpr

Fast floor operation, without checks, returning a dip::sint.

template<typename T, typename <SFINAE>>

dip::sint dip::ceil_cast( T v) constexpr

Fast ceil operation, without checks, returning a dip::sint.

template<typename T, typename <SFINAE>>

dip::sint dip::round_cast( T v)

Fast round operation, without checks, returning a dip::sint.

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.

template<typename T>

T dip::abs( T value) constexpr

constexpr version of std::abs. Prefer std::abs outside of constexpr functions.

template<typename T>

T const& dip::clamp( T const& v, T const& lo, T const& hi) constexpr

Clamps a value between a min and max value (a.k.a. clip, saturate, etc.).

dip::dfloat dip::pow10( dip::sint power) constexpr

Computes integer powers of 10, assuming the power is relatively small.

bool dip::ApproximatelyEquals( dip::dfloat lhs, dip::dfloat rhs, dip::dfloat tolerance = 1e-6) constexpr

Approximate floating-point equality: abs(lhs-rhs)/lhs <= tolerance.

dip::uint dip::LengthUnicode( dip::String const& string)

Counts the length of a (UTF-8 encoded) Unicode string.

dip::dfloat dip::BesselJ0( dip::dfloat x)

Computes the Bessel function J of the order 0 (with around 7 digits of precision).

dip::dfloat dip::BesselJ1( dip::dfloat x)

Computes the Bessel function J of the order 1 (with around 7 digits of precision).

dip::dfloat dip::BesselJN( dip::dfloat x, dip::uint n)

Computes the Bessel function J of the order n (with around 7 digits of precision).

dip::dfloat dip::BesselY0( dip::dfloat x)

Computes the Bessel function Y of the order 0 (with around 7 digits of precision).

dip::dfloat dip::BesselY1( dip::dfloat x)

Computes the Bessel function Y of the order 1 (with around 7 digits of precision).

dip::dfloat dip::BesselYN( dip::dfloat x, dip::uint n)

Computes the Bessel function Y of the order n (with around 7 digits of precision).

dip::dfloat dip::Sinc( dip::dfloat x)

Computes the sinc function.

dip::dfloat dip::Phi( dip::dfloat x)

Computes phi, the integral of the PDF of a Normal distribution with unit variance and zero mean from minus infinity to x.

dip::dfloat dip::Phi( dip::dfloat x, dip::dfloat m, dip::dfloat s)

Computes phi, the integral of the PDF of a Normal distribution with standard deviation s and mean m from minus infinity to x.

dip::dfloat dip::HypersphereSurface( dip::uint n, dip::dfloat r) constexpr

Computes the surface area of an n-dimensional hypersphere with radius r.

dip::dfloat dip::HypersphereVolume( dip::uint n, dip::dfloat r) 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.

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.

If n is 2 or 3, prefer to use dip::SymmetricEigenDecomposition2 or dip::SymmetricEigenDecomposition3.

void dip::SymmetricEigenDecomposition2( dip::ConstSampleIterator input, dip::SampleIterator lambdas, dip::SampleIterator vectors = nullptr, dip::Option::DecompositionMethod method = Option::DecompositionMethod::PRECISE)

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, dip::Option::DecompositionMethod method = Option::DecompositionMethod::PRECISE)

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.

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 , 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.

dip::uint dip::RankFromPercentile( dip::dfloat percentile, dip::uint n) constexpr

Computes the rank (index into array) for a given percentile and an array of length n.

The rank is symmetric (i.e. if the 5^th percentile translates to rank 14, then the 95^th percentile translates to rank n-1-14).

percentile is clamped to the range [0, 100], no error is produced for a percentile outside the valid range.

dip::dfloat const dip::pi = 3.14159265358979323846264338327950288 constexpr

The constant π.

dip::dfloat const dip::nan = std::numeric_limits::quiet_NaN() constexpr

A NaN value.

dip::dfloat const dip::infinity = std::numeric_limits::infinity() constexpr

Infinity.