Experimental 3x3 polar decomposition helpers used for color transform analysis. More...

#include <complex.h>
#include <math.h>
#include <stdbool.h>
#include <stdio.h>
#include <stdlib.h>

Include dependency graph for polar_decomposition.h:

This graph shows which files directly or indirectly include this file:

Macros
#define	TYPE double

#define	ABS(n) fabs(n)

#define	SQRT(n) sqrt(n)

Functions
void	swap_rows (double *m, size_t rows, size_t cols, size_t row0, size_t row1)
	Swap two rows of a row-major dense matrix in place.

void	swap_cols (double *m, size_t rows, size_t cols, size_t col0, size_t col1)
	Swap two columns of a row-major dense matrix in place.

void	abs_matrix (double *matrix_abs, size_t num_el)
	Replace each element of a vector with its absolute value.

void	matrix_multiply (const double A, const double B, double *RES, int rows_A, int cols_A, int cols_B)
	Multiply two dense row-major matrices.

double	max_val_matrix (double *matrice, size_t size)
	Return the maximum entry of a flat matrix buffer.

static void	normalize_array (double *vector, int size)
	Normalize a vector to unit Euclidean norm.

void	stampa_matrice (double *m, size_t rows, size_t cols)
	Print a row-major matrix for manual debugging.

void	compute_null_space (double *nullspace, const double a, const double b, const double c)
	Compute one non-zero vector in the null space of a symmetric 2x2 matrix.

void	copy (double dest, double source, size_t num_el)
	Copy a flat buffer.

void	orthonormalize_v_with_qr (double v0, double v1, const double v00, const double v10, const double v01, const double v11)
	Build two orthonormal candidate eigenvectors from symbolic QR factors.

void	multiply_il_v (double RES, const double IL01, const double IL02, const double IL03, const double IL12, const double IL13, const double v)
	Multiply a vector by the lower-triangular inverse factor used in the LDL solve.

void	multiply_id_v (double RES, const double ID00, const double ID11, const double ID, const double *v)
	Multiply a vector by the block-diagonal inverse factor used in the LDL solve.

void	multiply_v_il (double RES, const double v, const double IL01, const double IL02, const double IL03, const double IL12, const double IL13)
	Multiply a vector by the transpose-side lower-triangular factor used in the LDL solve.

void	orthonormalize_v_with_qr_single (double v0, double v1)
	Re-orthonormalize two 4D vectors after repeated inverse iteration steps.

void	multiply_minus_v_d (double RES, double v, double *D)
	Apply the negative D factor of the LDL factorization to a vector.

double	vector_scalar (const double a, const double b, const size_t num_el)
	Compute the scalar product of two vectors.

void	polar_decomposition (double A[3][3], double Q[3][3], double H[3][3])
	Compute the polar decomposition of a 3x3 matrix.

Detailed Description

Experimental 3x3 polar decomposition helpers used for color transform analysis.

This implementation follows the 3x3 polar decomposition strategy described in "An algorithm to compute the polar decomposition of a 3 x 3 matrix" by Nicholas J. Higham and Vanni Noferini, 2016.

Let A be a non-singular 3×3 matrice, like the ones used in channel mixer or in cameras input profiles. Such matrices define transforms between RGB and XYZ spaces depending on the vector base transform. The vector base is the 3 RGB primaries, defined from/to XYZ reference (absolute) space. Converting between color spaces is then only a change of coordinates for the pixels color vector, depending on how the primaries rotate and rescale in XYZ.

RGB spaces conversions are therefore linear maps from old RGB to XYZ to new RGB. Geometrically, linear maps can be interpreted as a combination of scalings (homothety), rotations and shear mapping (transvection).

But they also have an interesting property :

For any 3×3 invertible matrice A describing a linear map, the general linear map can be decomposed as a single 3D rotation around a particular 3D vector.

That is, there is a factorization of A = Q * H, where Q is the matrice of rotation around a axis of vector H.

This is interesting for us, on the GUI side. 3×3 matrices (9 params) are not intuitive to users, and the visual result of a single coefficient change is hard to predict. This method allows us to reduce 9 input parameters to :

6 : 3 angles of rotation, and the 3D coordinates of the (non-unit) rotation axis vector,
7 : 3 angles of rotation, the 3D coordinates of the unit rotation axis vector, and a scaling factor for this vector.

Usually, this is achieved by using HSL spaces, which suck because they work only for bounded signals in [ 0 ; 1 ]. Also, they are not colorspaces, not connected to either physics or psychology, so they are bad. Anyone saying otherwise is fake news.

The present method generalizes the HSL approach to XYZ, LMS and weird spaces, with none of the drawbacks of the the cheapo lazy-ass maths-disabled HSL bullshit. It's great. You should try it some time. Simply the best.

Reference paper: https://www.researchgate.net/publication/296638898_An_algorithm_to_compute_the_polar_decomposition_of_a_3_3_matrix/link/56e29d8c08ae03f02790a388/download

Reference implementation: https://github.com/higham/polar-decomp-3by3

Macro Definition Documentation

◆ ABS

#define ABS ( n ) fabs(n)

◆ SQRT

#define SQRT ( n ) sqrt(n)

◆ TYPE

#define TYPE double

Function Documentation

◆ abs_matrix()

void abs_matrix	(	double *	matrix_abs,
		size_t	num_el
	)

Replace each element of a vector with its absolute value.

Parameters

[in,out]	matrix_abs	Flat array to process in place.
[in]	num_el	Number of elements to process.

References ABS, and i.

◆ compute_null_space()

void compute_null_space	(	double *	nullspace,
		const double	a,
		const double	b,
		const double	c
	)

Compute one non-zero vector in the null space of a symmetric 2x2 matrix.

The coefficients correspond to the matrix [ a b ; b c ], which is the shape produced by the reduced system in the polar decomposition code.

Parameters

[out]	nullspace	Output vector of length 2.
[in]	a	Matrix coefficient (0, 0).
[in]	b	Matrix coefficient (0, 1) and (1, 0).
[in]	c	Matrix coefficient (1, 1).

Referenced by polar_decomposition().

◆ copy()

void copy	(	double *	dest,
		double *	source,
		size_t	num_el
	)

Copy a flat buffer.

Copy a flat buffer used to initialize test matrices.

Parameters

[out]	dest	Destination buffer.
[in]	source	Source buffer.
[in]	num_el	Number of elements to copy.

Referenced by _masks_gui_menu_item_forward_event(), _set_help_string(), _update_tree_view_nodes_font(), dt_cairo_sharpen_surface_rgb24(), dt_collection_load_filmroll(), dt_ioppr_extract_multi_instances_list(), dt_ioppr_merge_multi_instance_iop_order_list(), gui_init(), main(), and polar_decomposition().

◆ matrix_multiply()

void matrix_multiply	(	const double *	A,
		const double *	B,
		double *	RES,
		int	rows_A,
		int	cols_A,
		int	cols_B
	)

Multiply two dense row-major matrices.

The function assumes the matrix shapes are compatible and does not perform dimension validation. It is used as a small utility inside the decomposition code where the shapes are known statically at the call site.

Parameters

[in]	A	Left matrix stored row-major with shape `rows_A` x `cols_A`.
[in]	B	Right matrix stored row-major with shape `cols_A` x `cols_B`.
[out]	RES	Product matrix stored row-major with shape `rows_A` x `cols_B`.
[in]	rows_A	Number of rows in `A`.
[in]	cols_A	Number of columns in `A` and rows in `B`.
[in]	cols_B	Number of columns in `B`.

References A, B, i, k, RES, and TYPE.

Referenced by main(), and polar_decomposition().

◆ max_val_matrix()

double max_val_matrix	(	double *	matrice,
		size_t	size
	)

Return the maximum entry of a flat matrix buffer.

Parameters

[in]	matrice	Flat matrix storage.
[in]	size	Number of elements in `matrice`.

Returns: Largest value found in the buffer.

References i, max, size, and TYPE.

◆ multiply_id_v()

void multiply_id_v	(	double *	RES,
		const double	ID00,
		const double	ID11,
		const double *	ID,
		const double *	v
	)

Multiply a vector by the block-diagonal inverse factor used in the LDL solve.

Parameters

[out]	RES	Output vector of length 4.
[in]	ID00	Reciprocal of the first diagonal coefficient.
[in]	ID11	Reciprocal of the second diagonal coefficient.
[in]	ID	Pointer to the trailing 2x2 block stored row-major.
[in]	v	Input vector of length 4.

References RES, and v.

◆ multiply_il_v()

void multiply_il_v	(	double *	RES,
		const double	IL01,
		const double	IL02,
		const double	IL03,
		const double	IL12,
		const double	IL13,
		const double *	v
	)

Multiply a vector by the lower-triangular inverse factor used in the LDL solve.

Parameters

[out]	RES	Output vector of length 4.
[in]	IL01	Lower-triangular coefficient.
[in]	IL02	Lower-triangular coefficient.
[in]	IL03	Lower-triangular coefficient.
[in]	IL12	Lower-triangular coefficient.
[in]	IL13	Lower-triangular coefficient.
[in]	v	Input vector of length 4.

References RES, and v.

◆ multiply_minus_v_d()

void multiply_minus_v_d	(	double *	RES,
		double *	v,
		double *	D
	)

Apply the negative D factor of the LDL factorization to a vector.

Parameters

[out]	RES	Output vector of length 4.
[in]	v	Input vector of length 4.
[in]	D	Pointer to the 4x4 factor stored row-major.

References RES, and v.

◆ multiply_v_il()

void multiply_v_il	(	double *	RES,
		const double *	v,
		const double	IL01,
		const double	IL02,
		const double	IL03,
		const double	IL12,
		const double	IL13
	)

Multiply a vector by the transpose-side lower-triangular factor used in the LDL solve.

Parameters

[out]	RES	Output vector of length 4.
[in]	v	Input vector of length 4.
[in]	IL01	Lower-triangular coefficient.
[in]	IL02	Lower-triangular coefficient.
[in]	IL03	Lower-triangular coefficient.
[in]	IL12	Lower-triangular coefficient.
[in]	IL13	Lower-triangular coefficient.

References RES, and v.

◆ normalize_array()

static void normalize_array	(	double *	vector,
		int	size
	)

inlinestatic

Normalize a vector to unit Euclidean norm.

Parameters

[in,out]	vector	Vector to normalize in place.
[in]	size	Number of elements in `vector`.

References size, SQRT, and TYPE.

Referenced by orthonormalize_v_with_qr(), orthonormalize_v_with_qr_single(), and polar_decomposition().

◆ orthonormalize_v_with_qr()

void orthonormalize_v_with_qr	(	double *	v0,
		double *	v1,
		const double	v00,
		const double	v10,
		const double	v01,
		const double	v11
	)

Build two orthonormal candidate eigenvectors from symbolic QR factors.

Parameters

[out]	v0	First output vector of length 4.
[out]	v1	Second output vector of length 4.
[in]	v00	First symbolic coefficient.
[in]	v10	Second symbolic coefficient.
[in]	v01	Third symbolic coefficient.
[in]	v11	Fourth symbolic coefficient.

References normalize_array().

◆ orthonormalize_v_with_qr_single()

void orthonormalize_v_with_qr_single	(	double *	v0,
		double *	v1
	)

Re-orthonormalize two 4D vectors after repeated inverse iteration steps.

Parameters

[in,out]	v0	First vector of length 4.
[in,out]	v1	Second vector of length 4.

Todo:: : find and equivalent in C

References ABS, d, f, factor, g, normalize_array(), and TYPE.

◆ polar_decomposition()

void polar_decomposition	(	double	A[3][3],
		double	Q[3][3],
		double	H[3][3]
	)

Compute the polar decomposition of a 3x3 matrix.

The function factorizes the input matrix as A = Q * H, where Q is orthogonal and H is symmetric positive semidefinite. The implementation first normalizes the input matrix, then builds a reduced eigenproblem whose dominant eigenvector is converted into the quaternion-like parametrization used to reconstruct the orthogonal factor.

Parameters

[in,out]	A	Input 3x3 matrix, temporarily normalized in place and restored before returning.
[out]	Q	Orthogonal factor.
[out]	H	Symmetric factor.

References A, ABS, B, compute_null_space(), copy(), d, H, i, k, L, m, matrix_multiply(), normalize_array(), p, r, S, SQRT, swap_cols(), swap_rows(), t, TYPE, v, vector_scalar(), and x.

Referenced by main().

◆ stampa_matrice()

void stampa_matrice	(	double *	m,
		size_t	rows,
		size_t	cols
	)

Print a row-major matrix for manual debugging.

Parameters

[in]	m	Matrix data stored row-major.
[in]	rows	Number of rows.
[in]	cols	Number of columns.

References i, and m.

◆ swap_cols()

void swap_cols	(	double *	m,
		size_t	rows,
		size_t	cols,
		size_t	col0,
		size_t	col1
	)

Swap two columns of a row-major dense matrix in place.

Parameters

[in,out]	m	Matrix data stored row-major.
[in]	rows	Number of rows in `m`.
[in]	cols	Number of columns in `m`.
[in]	col0	Index of the first column to swap.
[in]	col1	Index of the second column to swap.

References i, m, and TYPE.

Referenced by polar_decomposition().

◆ swap_rows()

void swap_rows	(	double *	m,
		size_t	rows,
		size_t	cols,
		size_t	row0,
		size_t	row1
	)

Swap two rows of a row-major dense matrix in place.

Parameters

[in,out]	m	Matrix data stored row-major.
[in]	rows	Number of rows in `m`.
[in]	cols	Number of columns in `m`.
[in]	row0	Index of the first row to swap.
[in]	row1	Index of the second row to swap.

References i, m, row0, row1, and TYPE.

Referenced by polar_decomposition().

◆ vector_scalar()

double vector_scalar	(	const double *	a,
		const double *	b,
		const size_t	num_el
	)

Compute the scalar product of two vectors.

Parameters

[in]	a	First vector.
[in]	b	Second vector.
[in]	num_el	Number of elements in both vectors.

Returns: Dot product of a and b.

References TYPE.

Referenced by polar_decomposition().

Macros

Functions

Detailed Description

Macro Definition Documentation

◆ ABS

◆ SQRT

◆ TYPE

Function Documentation

◆ abs_matrix()

◆ compute_null_space()

◆ copy()

◆ matrix_multiply()

◆ max_val_matrix()

◆ multiply_id_v()

◆ multiply_il_v()

◆ multiply_minus_v_d()

◆ multiply_v_il()

◆ normalize_array()

◆ orthonormalize_v_with_qr()

◆ orthonormalize_v_with_qr_single()

◆ polar_decomposition()

◆ stampa_matrice()

◆ swap_cols()

◆ swap_rows()

◆ vector_scalar()