Alternative Image Stretching for Mathematically bent Pixinsight Users: from Linear to non-linear - quicker, easier,.... and better? [Deep Sky] Processing techniques · David Payne · ... · 77 · 11271 · 14

How do I get it:

The attached file, I am hoping, can be saved to where you normally keep your Process Icons (or wherever).  It defines a Pixelmath function that can stretch your images and can be loaded under the Process-Process Icons-Load Process Icons.  

GenHypStretchv1.xpsm

Alternatively, you can load the following into Pixelmath via copy and paste into the expression editor....  Two sections, between the, not include the ***s....

Into RGB/K:
*****
iif(b==0,EC=1,EC=0);
iif(b>0,Ds=D*b,Ds=D);
iif(b>0,bs=b,bs=1);
iif(EC==1,q0=exp(-Ds*SP),q0=(1+Ds*SP)^(-1/bs));
iif(EC==1,qWP=2-exp(-Ds*(HP-SP)),qWP=2-(1+Ds*(HP-SP))^(-1/bs));
iif(EC==1,q1=2-2*exp(-Ds*(HP-SP))+exp(-Ds*(2*HP-SP-1)),q1=2-2*(1+Ds*(HP-SP))^(-1/bs)+(1+Ds*(2*HP-SP-1))^(-1/bs));
iif($T<SP,EC*exp(-Ds*(SP-$T))+(1-EC)*(1+Ds*(SP-$T))^(-1/bs)-q0,iif($T>HP,2-EC*(2*exp(-Ds*(HP-SP))+exp(-Ds*(2*HP-$T-SP)))+(1-EC)*(2*(1+Ds*(HP-SP))^(-1/bs)+(1+Ds*(2*HP-$T-SP))^(-1/bs))-q0,2-EC*exp(-Ds*($T-SP))-(1-EC)*(1+Ds*($T-SP))^(-1/bs)-q0))/(q1-q0);
****
& into symbols:
****
D = 1.0;
b =0.5;
SP =0.00;
HP =1.00;
Rnorm;
q0;
qWP;
q1;
Ds;
bs;
EC;
*****

So what does it do?
You may be suprised at what these little equations can do for stretching...

I consider image stretching my make or break process in creating an image - unfortunately it is generally hit-or miss for me.   In an effort to potentially improve upon and parameterize image stretching, I have put together a Pixelmath sequence to stretch images in a natural, but controllable manner, that in my initial testing, is showing to improve upon the multiple use of curves or histogram stretch and is easier to duplicate than "hand drawn" curves, or my previously favourite method - Arcsinh stretch + Curves.  

The basis of this stretching are general hyperbolic decline decline/growth functions used in empirical forecasting of natural time series processes such as nuclear decay (exponential), chemical reaction rates (concentrations, catalyst site availability), oil,gas, and water well production decline (Arps (Arps, not Arp) equations), and even epidemiology. 

I believe that the general hyperbolic stretch routines are suitable to for stretching (even in extremes) due to the properties that: they are piecewise continuous given natural look to the result, their first derivative is always positive (will always maintain the rank order of pixel intensity - avoiding multiple stretch artifacts such as backwards stretches) and piecewise continuous (will avoid spike and backwards stretch situation that can arise with curves), are normalized to run between 0 and 1 and will never "clip" data.  In its basic form, the equations provide a stretch to the lower intensity pixels (increasing the contrast), while decreasing the relative contrast at high intensity pixels.   Note that this is always the trade-off with stretching, as whatever stretching function is applied from 0 to 1 intensity level, the average slope must equal 1. 

In base form, two parameters control the stretch: "D" and "b", with the maximum stretch intensity (dynamic range extension) occuring at 0 intensity (darkest part of image) and continuously reducing to the right hand side of the historgram.   In this way, the left side of the histrogram has its contrast between pixels increased while the right hand side has it increased.   D is used to control the amount of stretch that is applied in a give application.  The equations will accept any value of D>0 (with 0 meaning no stretch/change at all) and 10 being the maximum that I find you can routinely use while controlling the results and <1 used to tweak results and final adjustments.    Shown here, is the stretching function, with b=1 for various values of D - the degree of stretching:



If you have used curves before to stretch your image, you would likely have created such a curve within the Curves Process.   Particularly with large stretches, using handdrawn points joined by spline can lead to significant artifacts.   Not so with these equations, they are precisely and accurately defined everywhere and you will get an exact stretch within the digital precision being used.  

The other parameter "b" can be thought of as the "stretch intensity", or how focused the actual stretch is around a single point.   For bigger b, the stretch will be greater focused around a single intensity, while a lower b will spread the stretch around.   Mathematically, a b=0 represents a pure exponential stretch, while 0<b<1 represents a hyperbolic stretch, b=1 is a harmonic stretch, and b>1 is a highly intense, super-hyperbolic stretch.  I find keeping b<2 useful, unless the image has very poor dynamic stretch.   Here for a constant D of 2, is the effect of changing b:



As b is increase, the contrast around 0 pixel intensity is increased the most through range expansion - this is at the expense of constrast for high intensity pixels, where stretching reduces contrast.   Think of b and D as the trade-off between the two.  D is the blunt instrument, controlling how much stretch is used, and b is the surgical instrument - controllling how focussed the stretch is at, in this case, 0 intensity.  Early in the stretch, you will be inclined to to use a large factor for D and b. 


Two additional parameters are defined:  SP and HP.  SP which defaults to 0 specifies the pixel intensity at which the maximum stretch is applied.   Here is the effect of SP on the stretch function for various values of b and SP.



 The b factor will control how focused the stretch will be at that point, and the level of b can be determined by examining how spread out the histogram is.   If it is narrower, then use a larger value of b.   The way to use SP and b is either use the target image to find the intensity level where you would like to increase contrast the most, and set SP to this intensity level - It should also be the place on the histogram where the slope is greatest.  If you want narrowly increase the contrast at this point, use a higher b.   If it is a more general level of contrast increase, set b to a lower level.      Finally, the histogram peak (mode) location along with the image itself will indicate how much stretch to apply via D.



HP which defaults to 1, is a method to reserve dynamic range for the stars, or the higher end of the intensity charge.  Here some charts for various values for b and HP so you can see what the effect is.  If you want to increase contrast to the brighter parts the 



Stars can easily be overstretched, so it is not uncommon to use a HP value as low as 0.6 or 0.65 to protect high intensity stars with halos.   Ideally HP should be set between the values of the brightest nebula and the intensity of the dimmest halos in danger of bloating.  A level where there is already a large amount of contrast is best.  Alternatively, if stars are already are clipped at 1, or are close to it, using a value of HP>1, can create some room to stretch these stars by controlling their brightness (again no clipping by the stretch routine will occur).   Other ways of dealing with bright stars exist, such as removing the stars and stretching independently.   Another is to invert the image (so that bright stars appear dark), then stretch this low area using this Pixelmath routine,  trying to keep nebulosity intact, and then re-invert the image.

That's it, I could talk endlessly on how to do more complicated things with this stretch, but I was interested if anyone wanted to try it out.   I, myself, have just created this function, and will be trying it out (other than testing) on  some real image processing shortly.   If you would like to try it, I suggest trying it on some of your own unstretched images and play with it.   Please provide feedback though - and show me some images.

One last problem:

I would love to have these equations actually draw the curves (as you see above) like is done in the Curves Transformation Process in Pixinsight - show the before and after histograms (I now use Histogram Transformation for this), and actually show a preview image - much like Histogram.   However, I never learned Javascript (although I do know C++), and don't have the skills to put this into a Pixinsight script or process.   If I could, I would and I believe that would make this a great tool.   If you know how, and have some time, let me know.

I hope this is useful to you.

Dave

Dave,
That looks very interesting!  I’m going to have to play with it .  I agree with you that it would be nice to put into a tool where you could adjust the parameters and see the stretch curve before applying the transformation to the data.  (I’m in the same boat as you.  Long ago I wrote a lot of C++ code but my programming skills have since “aged out” so I’d have trouble implementing this sort of thing as a PI tool). Thanks for sharing your thoughts on this!

John

Thanks for the feedback.   On scripting, Pixinsight assumes you will be using one of the canned processes.   If you want to do something different, the scripting functionality/library is at a very low level.   For example, there is no "xy plot" facility in the library - you have to design it completely from scratch.   Way back in the early 90s, I was creating Unix/C programs with a graphical method of creating buttons, plots, slider etc - visual basic and .Net can graphically create dialogue boxes now too.   Not to knock Pixinsight too hard, because I love the program, but you really have to be a die-hard programmer to create a dialogue with widgets from their very low level Java script library.    Pixelmath would be the bee's knees if it a) could plot (simple x/y) the function you are entering (like the curves transformation plot) and b) have preview functionality (that you could also view on a histogram)… Then again why isn't there always a histogram tied to the current view!

Anyways, please try the Pixelmath equations/input - I am still coming up with ways to apply them (star protection, saturation adjustment, etc.) - and how to best locate optimal parameters.    If people try them and like them, it would go a long way to pushing the Pixinsight folks to make changes.

Cheers

Pulling @Roberto Sartori in the discussion, as this sounds like it's right up his alley!

@Salvatore Iovene thanks for bringing me in That's an interesting idea, I think that all the building blocks needed to make this a script are already available. 2D Plots and Image previews are used in several scripts so extracting the code and see how to use these components is quite easy, the rest is about putting the bricks together to build the wall. Let me brainstorm the idea:

Interface:
- you need an image preview control to preview the stretch
- you need a 2d plot component control to preview the stretching function
- you need few sliders that control  the stretch formula parameters
- you need a reset tool button to restore the curve default  values
- you need a cancel/execute button to finalize the operation on the working image or to cancel

When the script starts it could create a hidden copy of the active image, this will work as a preview image shown the preview control, and all default parameters must be set.

You have to build your "refresh" function that could work like this:

- instantiate PixelMath
- construct the stretch formula to be written into PixelMath using the current parameters
- set the formula in the PixelMath instance
- apply PixelMath to the preview image
- show the preview image in the preview control
- generate an array of values calculating the stretch formula at regular steps between 0 and 1
- plot the array of values as a 2D curve in the 2d plot component

Every time you change a slider (or reset them all)  you should update the corresponding formula parameters and re-execute the "refresh" function.

When you click the "execute" button you just need to create a PixelMath instance, set the formula and apply it to the working image.

On exit don't forget to delete the hidden preview image.

I hope this helps as a brief guide on how to do it at a very high level. For any specific questions regarding the implementation please visit PIxInsight's dedicated development forum section and you'll receive all support needed to make this idea a real script to be shared with the community!!!

Robyx

Benny Colyn:
I believe the original SubframeSelector script (before the SubframeSelector *Process*) displayed some graphs using SVG in the build-browser widget. Might be an option if you want code it.

Thanks Benny,  I will try and hunt it down.

Hi Dave

Really interesting approach you have here.  I am looking at putting together a script that would do what you want.  So I am trying to understand your equations but I notice that your calculated variable qWP does not get used in the final result anywhere - unless I have missed something.  Is that right?  Also, of course Rnorm isn't used either - I guess that is a hang over from a previous version.

All the best

Mike

I guess it is bad form to post to your own thread, but I thought I would give one more suggestion, I will try and keep it non-Pixinsight focused, and actually show some images!

Once pixel values are normalized to a blackpoint (or minimum point), then any transformation can be applied to the data, to increase the contrast in one part of the image in favour of the other.   The only real condition placed on this is that the transformation or function applied must take the data as an argument, and the result must always increase as the data increases.   For example, a sqrt relationship can emphasis dim objects using the tranformation y=sqrt(x), or in Pixelmath:
sqrt(($T-BP)/(1-BP))
which reduces to sqrt($T), if $T is linearly normalized with the black point first.
where BP is the black-point.    Or another (I haven't actually tried this one)
sin($T*3.14159/2)
or even another I haven't tried is
$T*sin($T*3.14159/2)
which also always increases with $T, or:
ln($T+sqrt($T^2+1))
which happens to be the formula for arcsinh stretching (although Pixinsight modifies this with a parameter, b to control the stretch level....
ln(b*$T+sqrt((b*$T)^2+1))/ln(b+sqrt(b^2+1))

Any of these transforms can be used in series, including in series with the general hyperbolic stretch shown here.  For Pixinsight this could amount to creating a pixelmath process icon, labelled with whatever you want.   Just remember some processes work best on linear data (such as Deconvolution), and once you apply one of these non-linear transforms, the data will be non-linear.   As an aside, if you for photometry that your camera gives non-linear results, if one could find the function that descibes how non-linear it is, then applying the inverse function would make the data linear - allowing the use of a non-linear camera to be used for photometry!).

As a example, I took the lum channel of the Ring Nebula (M57), shown here in linear form:


When stretched via the ScreenTransferFunction and HistogramTransferFunction method, both the halo around M57 and the galaxy IC1296 are revealed:

but, this overstretches the ring nebula itself blows out the stars somewhat.   But under a gentler HT stretch, that works for the ring and stars, hides the galaxy the M57's Halo:


So somewhere along these lines are the best one can do with histogram stretch alone.
Instead, I first put my unstretched image through a pixelmath transformation with the formula:
$T^P
where in this case, P was set to 0.25.   (any P<1 will emphase shadows over bright areas).   The result was the greyish image (because the background got stretched too):

so I judiciously chose a BP for blackpoint normalization, and then stretch using the generalized hyperbolic transform a couple of times, trying to pick out specific areas to increase constrast and protect the stars.   The final result of the Pixelmath stretches is:

So the result, I believe has preserved the stars and the ring itself and is the better than juse using histogram stretch alone.   I would, at this point deal with the dynamic range further with HDMRT and LHE in Pixinsight.   
Now, one could say that they could do as good or even better with CurvesTransformation - and I can't argue that - especially if done by someone with a more artistic eye than myself is doing it.   However, doing it with pixelmath allows one with a bit different form of creativity to come up with reproduceable transforms.

If you have read this far,  thanks.   I will leave it alone again - its just that this has been a lot of fun and I hope useful to you.

Clear Skies!
Dave

I just finished trying out this new GHS script a few days ago and I am very pleased with David and Mike’s results.  I’ll continue to use this script routinely in my image processing workflow.  As I use it more, I’m certain that I’ll learn some further advantages from the script that I haven’t yet tried . . . or even imagined.

I used it for all the linear-to-nonlinear stretching on this image.  There were 3 separate sub images on which I used it to perform the stretch (RGB, Lum and HOO).  I used only the first 3 of the 5 script parameters (D, b and SP).  

I immediately noticed the improved protection against star bloat that I could achieve with the “b” parameter.  Being able to see the image histogram while testing (and being able to easily undo) various values for the "SP" parameter, gave me much better visibility. understanding and control of where I was placing the point of maximum contrast in the image to minimize stretching the background noise and maximize the contrast on the subject of interest.  Having these "levers” and having the fundamental workings of the parameters explained by David in the pedagogical manner provided in his primer was very helpful.  It made me feel as though I actually knew what I was doing to my data and knew that I was giving it the best chance to shine.

I also used the GHS script for several noise reduction masks in processing this image.  Using the GHS script, I felt I was better able to put the mask strength where I needed it for better NR results.

Overall I’m very satisfied with the GHS script and I think David and Mike have a winner here.

A tutorial is coming soon.  Will post in the next couple of days with a link.

If possible could you please include an initial setup with a Mac? How and where to load the files to make it operational would be greatly appreciated. Kudos for this script!

You just install it like any other PI script. Here's how I did it:

1. Download the release from the GitHub link.
2. Unzip it if your Mac didn't by default
3. Copy the folder GHS-1.0 to /Applications/PixInsight/src/scripts (this step is not necessary, but I did it to ensure I didn't accidentally delete the script from my Downloads folder)
4. Fire up PixInsight
5. Go to Script -> Feature Scripts
6. Click the Add button and go to that GHS-1.0 directory and click Open
7. The script will be added
8. Click Done

That's it. You'll now see the GHS under Script -> Utilities.

Thanks JB…it worked!! 

i've never had to learn anything about computers past the copy/paste function and reading the internet so all this manipulation of files, creating folders and such I've had to finally learn! It takes a while at age 67! I really do appreciate and am grateful for help from people like yourself who are so willing to lend a hand to an old dolt as myself. Thank you.

Here is a half hour video tutorial on getting started with the GHS script.   It was taken from a Zoom meeting where I was showing some fellow members of my local centre (Victoria) of RASC how to get started.  I have put the video on Youtube.  Please excuse the videography and elementary editing.

https://youtu.be/ke5B8XQBwVc

I hope this helps.
Thanks,
Dave