Colours to matrix and polynomials

What -color-matrix or -function Polynomial will make one image look like another?

When we have two same-sized images with corresponding pixels, so there is a colour transformation but no geometric transformation, we can calculate the parameters for "-color-matrix" or "-function Polynomial" that approximate the transformation.

The method assumes the transformation is global rather than local, eg the entire image has been shifted towards red, rather than one part shifted towards red and another part towards blue.

Scripts on this page assume that the version of ImageMagick in %IM7DEV% has been built with various process modules. See Process modules.

References

The method

IM's "-color-matrix" and "-function Polynomial" operations have mathematical definitions: each channel of each output pixel is calculated from a formula that has the channels of the corresponding input pixel and a list of numbers. We use these operations by supplying the list of numbers, so IM can calculate the output pixels.

This method operates the reverse process: given the input and output pixels, it calculates the list of numbers.

How does it do this? There is one equation per channel per pixel, with known channel values and unknown numbers in the list, but we know the same numbers are used in all the equations. So we have a bunch of simultaneous equations that can be solved for all the unknowns, which gives us the list of numbers. We often have more equations than unknowns, so we find the unknowns that minimizes the sum-of-squared-errors.

Why do we want this? We might have a set of input images that we want to normalise to a common standard. For example, time-lapse photos of a building site that we want to make into a movie. Or aerial photos that we want to join together. Or photos of an object taken under different lighting conditions. The same matrix or polynomial is then applied to any image for which we want the same transformation.

Before giving the method details, we will play with IM's colour matrix and polynomial to see how they work.

The following descriptions assume the three colour channels are RGB, but they could be Lab or YIQ or any other three-channel colorspace.

IM's colour matrix

IM's operation "-color-matrix" calculates new pixels as the sum of a DC offset plus the colour components each multiplied by a factor. Hence, it provides for cross-feed between channels.

For RGB with offsets, we use a 6x6 matrix. The six columns and rows are the five channels (R, G, B, K, A) and an offset. "K" is the black channel of CMYK, and we don't use this. The top three rows specify the red, green and blue outputs. The bottom three rows are constant, and we don't usually want to change colours according to alpha, or to change alpha according to colours, so we have only 12 important values:

When we know the values of a..m, each output R',G',B' is then calculated from the corresponding input R,G,B:

If we use only the diagonal elements a,f,k and the offsets d,h,m, setting other values to zero, we have the equivalent of the gain and bias method, where a,f,k are the gain multipliers and d,h,m are the bias offsets.

How do we get the 12 numbers a..m from two images? If the images each have N pixels, then we have N input values for each of R, G and B, and N output values for each of R', G' and B'. This gives us 3*N simultaneous equations that can be solved by Gauss-Jordan elimination.

For example, suppose an input pixel is RGB=(10%,20%,30%) and the corresponding output pixel is R'G'B'=(21%,22%,14%). This gives us three equations:

For this pair of pixels, the terms are (10,20,30,1), and the results are (21,22,14).

Four pairs of pixels would supply 12 simultaneous equations, so we could solve for the 12 unknowns a..m. If we have more than four pixels, the problem is over-determined and there may be no exact solution. But we can find a solution that minimizes the squared errors.

The output in each channel depends on the input in every channel, so we can't represent the transformation by simple "in versus out" curves.

Constraint: no cross-feed

We may wish to prevent the solution from cross-feeding between channels. The solution then has only six important values, gain a,f,k and bias d,h,m:

The equations are fully independent; they have no elements in common. Hence they are solved as three independent problems.

This is simple linear regression. The output in each channel depends on the input of that channel only, so we can represent the transformation by three simple curves (which are all straight lines):

Constraint: gain only

Constraining to just the gain can be useful. The solution then has only three important values, a,d,f:

Again, the equations are fully independent, so they are solved as three independent problems.

IM's polynomial

IM's operation -function polynomial is closely related to -colour-matrix. We can apply a different polynomial to each channel, and that is how we will use it. There is no cross-channel mixing. Each output channel value is the sum of the input raised to a number of integer powers, each multipled by a coefficient.

If u is the input value normalised to typically between 0.0 to 1.0, and the polynomial is degree n, then the output u' is...

Given a number of inputs u and corresponding outputs u', the module calculates the coefficients a₀ ... a_n. The number of terms, and number of coefficients, is (n+1). We have one set of coefficients per channel, for a total of 3*(n+1) coefficients.

A polynomial of degree zero adds a constant value (a bias-only operation). A polynomial of degree one multiplies, then adds a constant value (a linear polynomial; a gain and bias operation). Polynomials of degree two are called quadratic; of degree three are called cubic. Higher degrees describe more complex transfer curves. A polynomial degree n can be used to transform an image of n pixels to any other image.

By applying the same operations to a gradient, we can get a transfer curve for each channel.

So, we can create a polynomial for each channel, and tweak a photo so its colour patches roughly match those of a reference. Higher degrees of polynomials will match the patches more closely, until perfection when the degree equals the number of patches. But this comes at a price: with higher degrees, the curve becomes increasingly eratic at colours other than the patches. In practice, degrees above three do not seem useful, and degree two is often sufficient.

Example photos

The process module cols2mat

When we have two same-sized images, "-process cols2mat" calculates the 12 numbers a..m and hence the 6x6 colour matrix, or the polynomial of required degree, that best transforms the first image to the second.

Option	Description
Short form	Long form
`m string`	`method string`	Method for calculating the matrix, one of: `Cross` include cross-channel multipliers (12 terms); `NoCross` exclude cross-channel multipliers (6 terms); `NoCrossPoly` polynomial without cross-channel (3*`degreePoly`+3 terms); `GainOnly` include only this-channel multipliers (3 terms). Default = `Cross`.
`d integer`	`degreePoly integer`	For method `NoCrossPoly`, degree of polynomial. For example: degree 3 gives v' = av³ + bv² + c*v + d. Default: `2`.
`w number`	`weightLast number`	Weight for last line of image. For example: more than 1.0 (eg 10, 100) to give greater weight to last line, between 0.0 and 1.0 to give less weight. Default: `1.0`.
`wa`	`weightAlpha`	Multiplies weight by product of the pixel alphas.
`x`	`noTrans`	Don't replace images with transformation.
`f string`	`file string`	Write text data (the colour matrix) to `stderr` or `stdout`. Default = `stderr`.
`v`	`verbose`	Write some text output to stderr.
	`version`	Write version information to stdout.

The module needs two or three input images. The first two inputs must be the same size as each other. It replaces all the inputs with a single output image. It calculates the colour matrix or polynomial from the first two inputs. If there are only two inputs, the output is the first transformed by the colour matrix or polynomial. If there are three inputs then the output is the third transformed by the colour matrix or polynomial.

The noTrans option will leave the image list unchanged. For applications that don't need a transformed image, it saves some time.

The inputs must have three colour chanels, representing RGB, L*a*b*, YIQ or whatever. For RGB images, I use it with sRGB colorspace, but I expect it will work with any profiled 3-channel colorspace.

Typically the first two inputs are small, for example 6x4 pixels. But they can be any size, provided they are the same size and the pixels correspond. It takes about four seconds to process a pair of 35 MP images.

For the NoCrossPoly method, the text output is three lists of polynomial coefficients. For the other methods, the text output is a single list of numbers in the 6x6 matrix. (A more general polynomial with cross-channel terms is possible, but IM has no operation for this.)

By default the calculation gives equal weight to all the input pixels, but a different weight may be applied to the last row on input pixels, and the weight may be multiplied by the product of the alphas.

Extracting the matrix

The calculated 6x6 color matrix is sent as text to stderr or stdout, in a line that starts with "c2matrix=". The 36 numbers are separated by commas, with no spaces.

Extracting the polynomial

Testing the matrix

We test a round trip: now we know the matrix that most closely transforms imgA to imgB, we can apply the matrix to imgA and the result should be close to imgB.

We also test a case where the inputs are equal. This should create the identity matrix, which has one in the diagonal elements and zero elsewhere.

Limited arithmetic precision has prevented the cross-channel multipliers from being exactly zero. "method NoCross" will calculate the solution with zeros in those positions.

For some cases, the module will not find a unique solution, usually because the problem is under-constrained (there are fewer simultaneous equations than unknowns). For example:

No output line starts with "c2matrix=". Scripts should check for this condition.

(In fact, two solutions are possible: the blue output is the red input multiplied by one, or the blue output is set to an offset of one. Each of these solutions has an infinite number of variations.)

A single solution is found: the blue output is set to an offset of one, with no channel multipliers.

Testing the polynomial

We test a round trip: now we know the polynomials that most closely transforms c2mp_ph1_mat.png to c2mp_ph4_mat.png, we can apply them to c2mp_ph1_mat.png and the result should be close to c2mp_ph4_mat.png.

We also test a case where the inputs are equal. This should create polynomials with one in the v¹ elements and zero elsewhere.

No output line starts with "PolyRed=" etc. Scripts should check for this condition.

A single solution is found: the blue output is set to an offset of one, with no channel multipliers.

Weighting by alpha

The weightAlpha option calculates the solution of the simultaneous equations weighted by the product of the alphas of the input and output pixels. So if either pixel is entirely transparent, that pixel will be disregarded for the purpose of calculating the matrix or polynomial.

(Note that weightAlpha doesn't assume that pixel colours are calculated from alpha, nor that alpha is calculated. That would treat alpha like the colour channels, and I haven't included that in the module. I might in the future.)

What colour matrix makes the grass of toes.png look like the grass of toes_x_holed.png?

Oops. Most of toes_x_holed.png is (transparent) black, so the matrix is heavily influenced by this. We can eliminate the transparent pixels from the calculation by using the weightAlpha option:

Even when both inputs are in the range 0 to 100%, the transformed output can be outside that range. If the image is saved to an integer format, values will be automatically clamped (aka clipped). If HDRI is used, pixels can be brought within gamut by "-clamp" or "-auto-level".

As in the Wolf reference, the module could be run for a pair of opaque inputs, then the difference between the output and second input is negated and used as the opacity for the first input of a re-run, and iterate until stability is reached. Thus outliers would count less towards a solution.

Future

IM has no built-in operator to process these, so that would also need to be written.

We could "regard alpha". This could calculate and use alpha in the same way as the colour channels.

We might have an option to save the colour matrix as a 6x6 image. Then another module can apply it to images. And another can do maths: concatenation.

All images on this page were created by the commands shown, using:

%IMG7%magick -version

Version: ImageMagick 7.1.0-42 Q16-HDRI x64 396d87c:20220709 https://imagemagick.org
Copyright: (C) 1999 ImageMagick Studio LLC
License: https://imagemagick.org/script/license.php
Features: Cipher DPC HDRI OpenCL 
Delegates (built-in): bzlib cairo freetype gslib heic jng jp2 jpeg jxl lcms lqr lzma openexr pangocairo png ps raqm raw rsvg tiff webp xml zip zlib
Compiler: Visual Studio 2022 (193231332)

%IM7DEV%magick -version

Version: ImageMagick 7.1.0-20 Q32-HDRI x86_64 2021-12-29 https://imagemagick.org
Copyright: (C) 1999-2021 ImageMagick Studio LLC
License: https://imagemagick.org/script/license.php
Features: Cipher DPC HDRI Modules OpenMP(4.5) 
Delegates (built-in): bzlib cairo fontconfig fpx freetype jbig jng jpeg lcms ltdl lzma pangocairo png raqm rsvg tiff webp wmf x xml zip zlib
Compiler: gcc (11.2)

Source file for this web page is col2mp.h1. To re-create this web page, execute "procH1 col2mp".

This page, including the images, is my copyright. Anyone is permitted to use or adapt any of the code, scripts or images for any purpose, including commercial use.

Anyone is permitted to re-publish this page, but only for non-commercial use.

Anyone is permitted to link to this page, including for commercial use.

Identity matrix. %IMG7%magick ^ toes.png ^ -color-matrix ^ 1,0,0,0,0,0,^ 0,1,0,0,0,0,^ 0,0,1,0,0,0,^ 0,0,0,1,0,0,^ 0,0,0,0,1,0,^ 0,0,0,0,0,1 ^ c2mp_cm1.png
Swap red and green channels. %IMG7%magick ^ toes.png ^ -color-matrix ^ 0,1,0,0,0,0,^ 1,0,0,0,0,0,^ 0,0,1,0,0,0,^ 0,0,0,1,0,0,^ 0,0,0,0,1,0,^ 0,0,0,0,0,1 ^ c2mp_cm2.png
%IMG7%magick ^ toes.png ^ -color-matrix ^ 1,0.25,0,0,0,-0.25,^ 0,1,0,0,0,0,^ 0,0,1,0,0,0,^ 0,0,0,1,0,0,^ 0,0,0,0,1,0,^ 0,0,0,0,0,1 ^ c2mp_cm3.png
%IMG7%magick ^ toes.png ^ -color-matrix ^ 1.1,0.25,0,0,0,-0.1,^ 0.1,1,0.1,0,0,-0.2,^ -0.1,0,1,0,0,0,^ 0,0,0,1,0,0,^ 0,0,0,0,1,0,^ 0,0,0,0,0,1 ^ c2mp_cm4.png

set imgA=toes.png
set imgB=toes_x.jpg