Arcsinh like stretching in Affinity Photo

[windjammer]

Hi all
So, on a dreary Sunday I settled down to get my brain around arcsinh stretching and have a play to see what its all about.  I'm using Affinity Photo and got a blank - it doesn't do it.  So here is a close (?) work around for arcsinh and a preset setup for a live filter layer tool that does the business. I've got a headache now but wanted to put this out before the urge passes.

Comments welcome

Simon

Arcsinh stretching is used to stretch the image histogram while preserving the colour in the image. Generally the idea is to define a function f(x) so that:

x' = f(Z*x) where x' and x are the new and old pixel values for R,G, B and L, luminance, and Z is a parameter to control the strength of the effect (a larger Z gives a bigger effect).  

The function f(x) in the Affinity world has a value in the range 0,1, as do all the R,G,B and L values across the image.

The colour preserving bit of the arcsinh recipe comes from swapping R,G,B for L in the argument to f(x),
so that:

R' = (R/L)*f(Z*L), G' = (G/L)*f(Z*L) and B' = (B/L)*f(Z*L), and in the arcsinh method, 

f(Z*L) = arcsinh(Z*L)/arcsinh(Z), with Z in range 1 to 1000s.  The first pic shows the family of curves f(Z*L) for a range of Z values, 2 to 500.

Now, the arcsinh function, asinh(x) = ln(x+sqrt(x*x+1)), where ln is the natural log function.  Unfortunately Affinity Photo's version of pixelmaths does not supply a ln function so it cannot be computed in AP.

So, are there any curves that are similar in shape to the arcsinh family that could be used instead? The second pic shows the family of curves for x^p, where p is in range 0.17 to 0.8, which have similar shapes to arcsinh with Z in range 2 to 500.  So we can supply a new f(x):

f(L) = L^p, and R' = (R/L)*L^p, G'=(G/L)*L^p and B'=(B/L)*L^p

1. To add this into AP, open an image in AP.  Click Layer->Live filter layer->Colours->Procedural Texture.

2. Click on the '+' symbol on the left hand side, half way down the Procedural Texture (PT) window.

3. Paste the following string into the equation window, making sure the R channel box is selected:
var lum=clamp(0.2126*R + 0.7152*G + 0.0722*B-b, 0, 1);(R-b)*pow(lum,a)/lum

4. Click on the '+' symbol again to create a new equation window, with the G channel box selected. Paste in the same string with G substituted for R:
var lum=clamp(0.2126*R + 0.7152*G + 0.0722*B-b, 0, 1);(G-b)*pow(lum,a)/lum

5. And '+' again to create the Blue channel equation:
var lum=clamp(0.2126*R + 0.7152*G + 0.0722*B-b, 0, 1);(B-b)*pow(lum,a)/lum

6. Now create user control variable a:  At the bottom left of the PT window, click on the 0,1 symbol.  A symbol 'a' appears and a title box to the right.  Type in a useful name such as 'power to raise Lum'.

7.Now create user control variable b:  Click on the 0,1 symbol again.  A symbol 'b' appears and a title box to the right.  Type in a useful name such as 'black point subtract'.

8. Go to the top right of the PT window and click on the small pull down menu to the left of the Preset name window.  Select 'create preset' and type in useful name, eg 'arcsinh'.  The now created preset macro can be used again and in the future to do arcsinh-like stretching. If the Procedural Texture window does not behave as expected, the chances are there is an error in the equation strings, so check the syntax carefully.

9.  Dismiss the PT window.  To invoke, select the layer you want to process.  Layer->Live filter layer->Colours->Procedural Texture.  Click on the Preset pulldown, and select the arcsinh preset, (or whatever you called it).

10. To use the preset, first move the 'a' control to the far right.  This sets 'a', the power to raise L, to value of 1: L is thus unchanged.  Next move the control 'b' from the left (b=0) to the right so as to set a black point - watch the histogram window to locate the histogram near the left hand side of the graph.  

From the family of curves you can see that very strong stretching occurs at the smallest values of luminance, so setting the black point enables you to decide how much of the crud at the bottom you are going to stretch.  After that, move control 'a' left (p ->0) or right (p->1) to get desired stretching.  Smaller values of p give greater stretching.  Values between 0.4 and 0.8 estimated from the bar slider seem most useful.

What do the equations do:
var lum=clamp(0.2126*R + 0.7152*G + 0.0722*B-b, 0, 1);(R-b)*pow(lum,a)/lum

1. var lum=clamp(0.2126*R + 0.7152*G + 0.0722*B-b, 0, 1);

This computes the luminance value at a pixel from the R,G,B values, according to a recipe off the internet. The value of 'b' set by the slider subracts the black value, so luminance is zero at the chosen black value.  The clamp function just constrains the value of lum between 0 and 1.

2. (R-b)*pow(lum,a)/lum

This computes the black corrected R,G,B value to luminance ratio, before multiplying with luminance to the power set by 'a'.  So we get for each channel R' = (R/L)*L^p, G'=(G/L)*L^p and B'=(B/L)*L^p as discussed earlier.

ends

[vlaiv]

1 hour ago, windjammer said:

Arcsinh stretching is used to stretch the image histogram while preserving the colour in the image. Generally the idea is to define a function f(x) so that:

x' = f(Z*x) where x' and x are the new and old pixel values for R,G, B and L, luminance, and Z is a parameter to control the strength of the effect (a larger Z gives a bigger effect).  

The function f(x) in the Affinity world has a value in the range 0,1, as do all the R,G,B and L values across the image.

The colour preserving bit of the arcsinh recipe comes from swapping R,G,B for L in the argument to f(x),
so that:

R' = (R/L)*f(Z*L), G' = (G/L)*f(Z*L) and B' = (B/L)*f(Z*L), and in the arcsinh method, 

f(Z*L) = arcsinh(Z*L)/arcsinh(Z), with Z in range 1 to 1000s.  The first pic shows the family of curves f(Z*L) for a range of Z values, 2 to 500.

Why is this "method" special?

It is very clear that it is simply "repacked" simple RGB ratio method.

If R' = (R/L) * f(L)

(I'm omitting that parameter Z as it is really unimportant for the argument)

that can be rearranged as

R' = R * ( f(L) /L )

Similarly if we do the same for other two components we get:

G' = G * ( f(L) / L )

B' = B * ( f(L) / L )

f(L) / L is simply some g(L), and can be considered constant on pixel level (single luminance value per pixel)

We are simply multiplying every component with a constant (on pixel level) which preserves R:G:B ratio of that pixel (this is desirable quality).

Choice of f and consequently g is irrelevant, and we can use power law (gamma) or any other histogram transformation on luminance we desire and still preserve color in the same manner.

 

 

[windjammer]

So, as you point out, a colour preserving transform, which is (as I understand it) one of the plus points of the arcsinh transform, and that, we agree, is how it does it - and is of course entirely general for any constant across the channels.   Colour gets lost in endless cycles of manually driven curves.  (I'm not quite sure why that happens, but it does).  

The  other advantage of arcsinh is (I'm told) that it gets to a stretched result in less iterations than curves cycles, and being algorithmic is (in principle - pesky sliders reappear) reproducible, which is my interest.    Who really thinks that after  grinding round a dozen curves iterations, shifting points around on the curve, you would get the same result if you started from scratch? Similar, yes, but not the same.

But the point of the post was to find an alternative to arcsinh, which seems to be very popular for the above reasons, that was within the computational capabilities of Affinity Photo and would give similar results.   So a simple function along the lines of x^p which has similar properties to arcsinh, and integrates easily into a non-destructive live layer with only 2 controls: black point and stretch.  What's not to like ?

Simon

[vlaiv]

9 hours ago, windjammer said:

Colour gets lost in endless cycles of manually driven curves.  (I'm not quite sure why that happens, but it does).  

Because important step of RGB ratios is not performed

9 hours ago, windjammer said:

The  other advantage of arcsinh is (I'm told) that it gets to a stretched result in less iterations than curves cycles, and being algorithmic is (in principle - pesky sliders reappear) reproducible, which is my interest.    Who really thinks that after  grinding round a dozen curves iterations, shifting points around on the curve, you would get the same result if you started from scratch? Similar, yes, but not the same.

Any function that looks like this:

or to put it into words - operates on domain [0-1] and transforms it in [0-1] range and is monotonically increasing, will give you results that you need.

It does not have to be arcsinh - it can be any other function. It can be "free formed" - that is what curves do, and no, you don't have to do multiple rounds - you can only do a single round to reproduce exact result, or you can do analytical function.

What is important is to calculate "ratio" or coefficient c for each pixel that is in form - transformed_lum / lum, or to put in words - "how much we amplified lum for that pixel", and then we use that amplify each component of RGB.

As a final step - I would advocate exact same thing be done not in RGB space, but instead in XYZ color space, and then result transformed into suitable working color space (usually sRGB, but could be wide gamut like Adobe RGB or other if need be).

[ollypenrice]

My problem with all of this is that I only perform a global stretch up to the point at which the background has reached a desirable level. After that, I stretch manually but only above the background level. I can see no point in continuing to stretch the background (where signal is weak and the noise floor low) once it's reached its final level.

Olly

Edited December 6, 2022 by ollypenrice

[Martin Meredith]

An exceedingly simple function I use for stretching is

   f(x) = a * x / (x + b)

for suitable choice of a and b. I've not seen it used elsewhere in an image context (*).

A quick comment on the x^p family: it gets noisy near x=0 due to the gradient going to infinity. There is a modified gamma function that handles this by linearising the low intensity part of the range, although it adds a little to the computational complexity. I believe this also addresses the point Olly is making to some extent.

Martin

 

(*) This function was inspired by a 1985 approximation to the input-output function of the turtle cochlea hair cell that I picked up in the course of my doctoral research decades ago. The auditory system has to solve a similar problem to AP: compressing an enormous range of values into more meaningful equally-discriminable units.

[vlaiv]

13 minutes ago, Martin Meredith said:

An exceedingly simple function I use for stretching is

   f(x) = a * x / (x + b)

for suitable choice of a and b. I've not seen it used elsewhere in an image context (*).

It is actually basis for DDP - digital development process, one of first stretching functions developed to mimic analog film:

http://www.astrosurf.com/buil/us/iris/iris5.htm

where

a would be coefficient * level

and b = level

 

[windjammer]

>>after that, I stretch manually but only above the background level.

yes, the control described here has a black point set so that you can decide how much of the crud at the bottom you want to play with.  

>>f(x) = a * x / (x + b)

Thats interesting, I will give it a try!  Do you have a feeling for what a and b correspond to, or just take them as magic numbers ?

Simon

[Martin Meredith]

After some experimentation I've been using 

b = (1.02 - y) / 50

where y varies from 0 to 1 to control the degree of stretch (1 being most)

and a = (1 + b) to ensure the output goes from 0-1 too

These values give a pretty aggressive stretch for my (EEVA) purposes, but I suspect if the goal is to match arcsinh other values will be needed. No doubt some optimisation routine could be used to find the best match, or perhaps even some progress could be made analytically.

 

Edited December 6, 2022 by Martin Meredith
typo

[windjammer]

Hi Martin

That is a *very* fierce stretch - see first pic.  After playing around with the parameters I went back to your basic f(x) = a * x / (x + b).  A bit of a struggle with normalisation I ended up with f(x) =x* (1+b)/(x+b), b:(0.1).  See 2nd pic.

This looks very workable - looks as good as arcsinh or x^p, but very non linear in its variation with b.  So I think a tool would come in 2 versions.  Strong XXX, where b is scaled 0 to 0.15 and Auntie where b is scaled 0.2 to 1.

I think I am still waiting for the magic function where the whole stretch bore is done and dusted in one click.  Maybe AI is the answer...

Simon

[Martin Meredith]

Nice work Simon -- I'll take a look at 'Auntie' myself. Yes, fierce is the word. It really brings out faint galaxy features/tidal tails etc when used interactively in an EEVA-observing context but you wouldn't want it in a prize-winning AP situation, and most likely does crazy things to colour too.

[vlaiv]

27 minutes ago, Martin Meredith said:

and most likely does crazy things to colour too.

As long as you apply RGB ratio in the same way - it will preserve physical color.

[iantaylor2uk]

Not sure if there is any value using an S-curve. These have the properties that the gradient is zero when x=0 and x=1 which may or not be beneficial. 

I've just published a paper on using such curves to predict friction (the paper is open access if you want to read it at: https://www.sciencedirect.com/science/article/pii/S0301679X2200408X)

Based on the work in that paper you could try a curve such as: (1+(1/(cx))^k)^(-a) where k c, and a are constants (with c being 10 or 100 for example). 

[windjammer]

Hello Ian

Thats an interesting function - see pic.  There is a possibility of a set of curves in there that deal with the foot of the histogram in a more natural way than just cutting it off at the knees.  And also a set of curves that naturally move the histogram to the left after a stretch.  It needs a play to narrow the parameter space a bit.

Simon

[vlaiv]

@windjammer

What are you hoping to achieve with any given function?

Maybe we can define a set of requirements and then simply find function that behaves like that?

Here is an example, that might be useful for automatic stretching.

Background level and standard deviation is calculated.

"First part" of our stretch function will position background level at some predefined (or maybe parameter) value - like %8 or 0.08 (if we observe 0-1 range) - function itself will be linear in 0 - 0.08+3*stddev range - having that slope. From that value on we "connect" it to power low with such properties that:

1. power function at 0.08 + three standard deviations matches in value to linear first part

2. gradient at 0.08 + three standard deviations matches gradient of linear first part

3. power part reaches 1 at saturation value (again some value measured from the image to be saturation value)

[windjammer]

V - Still thinking abut this.  I have some formalism but doesn't go very far beyond stating the obvious:

Say G(x) is the transformed histogram H(x), G(x) = H(g(x)), g(x) the stretch function we want.

The derivative G' = (dH/dg).(dg/dx)

We specify an Xb, the bottom of the histogram below which we choose some linear function.  

We want G(Xb) = H(Xb) and G'(Xb) = H'(Xb) = (dH/dg).(dg/dx)

Which doesn't offer many more clues on the g(x).  I might play with adhoc methods to join up a linear bottom end with an interesting g(x).  Perhaps just join up with something like quadratic interpolation through a couple of points - that would satisfy your G=H and G'=H' at Xb criterion.

If it looks promising try and regularise it afterwards.

An Xt can be defined as you suggest, the top of the histogram at which we want the pixel to be saturated,  so that g(Xt) = 1. 

And all the fun is in between.

Simon

[vlaiv]

47 minutes ago, windjammer said:

V - Still thinking abut this.  I have some formalism but doesn't go very far beyond stating the obvious:

Here is what I had in mind for automatic stretching. I'll work with example as it might be easier to understand (I might lack proper terminology to explain things precisely).

Let's say that we have 16 bit data to start with. We need to determine three numbers in order to derive histogram transform:

1. Background level. Say that we identify background level to be around 900ADU

2. Standard deviation of pixel values around that background level (we identify pixels that we think are background and we compute standard deviation for sample) - say it is 20ADU

3. Maximum brightness value that won't saturate in the image - let's say we have nebulosity peak value at say 17000ADU. Stars can saturate - we don't care about those as they don't contain too much information, but nebula core does.

Now we have couple of equations to solve.

First is our background level - which we can set to say 8% of 0-1 range.

So we know that 900ADU maps to 0.08

Next, we want background noise to be virtually invisible. We ensure this by saying that +/- 3 standard deviations around that background value (3 here can be parameter) will fall into range that we can barely detect. We can barely detect say 0.5% difference in output value - or 0.005 out of 0-1 range.

900 + 60 = 0.08 + 3 * 0.005

900 - 60 = 0.08 - 3 * 0.005

We have our two points for linear part.

Now we have to connect that to a power law.

First we define general shape as base^(x*a + b)

Where X is our ADU value.

We now need several equations to form our system of equations to solve.

First condition is that for 17000 we want to have output value of say 0.96 (again - this value can be tweaked - but that is brightest part of target).

so

0.96 = base^(17000*a + b)

Second equation is that we want to "connect" above linear function to our power function at some point - let's say that is after those 3 standard deviations.

0.095 = base^(960*a + b)

Third requirement is that we have gradient of both linear function and our power law be same at that point so we have "smooth" joint.

Gradient is dy / dx so 6 * 0.005 / 6 * 20 = 0.005 / 20 = 0.00025

Now we find first derivative of our power law that will be

base^(x*a+b) * ln(base) * a

And that must be equal to 0.00025 for X = 960 so we plug in numbers for third equation:

base^(960*a+b) * ln(base) * a = 0.00025

Now, if we can solve this set of equations for base, a and b - we have our histogram transform that will satisfy our initial conditions.

[Martin Meredith]

The modified gamma function is worth considering perhaps 

        s = g / (a0 * (g - 1) + a0 ** (1 - g))
        d = (1 / (a0 ** g * (g - 1) + 1)) - 1
        y[x < a0] = x[x < a0] * s
        y[x >= a0] =  (1 + d) * (x[x >= a0] ** g) - d
 

where g is the exponent and a0 is the x-value below which the regularised part comes in (s is the slope).

This is Python code so ** is ^

I use this in my EEVA software to stretch the RGB channels semi-automatically (user supplies g only) so that the user doesn't have to mess around with the individual channels and can focus on stretching luminosity.