Gradient HDR compression

[georg.viehoever]

Hi,

I have been working on  a process for enhancing faint gradients, while suppressing strong gradients. It should be ideal for enhancing faint structures, especially in imkages with high dynamic range such as M42, M31, NGC7000, ... . It is an algorithm from the family of gradient domain algorithms, that I have already been playing with a while ago with the MergeMosaic script http://pixinsight.com/forum/index.php?topic=2124.0.

The GradientHdrCompression process works by computing the gradients in an image, and then manipulating the gradients. The manipulated gradients are then used tor reconstruct the image. The process allows 3 kinds of manipulations:
- scale gradients by pow(gradient, exp). If exp<1.0 this enhances faint gradients
- clip large gradients. The effect is similar to the pow-approach
- as an experiment, I also permit to clip low gradients to 0.0. I thought this might be useful for supressing noise, but my first experiments are inconclusive.
The remaining two parameters allow to influence the behaviour of the linear solver that is below the process. The default settings should be ok for many uses.

Attached is a screenshot of a results I obtained using HDR luminance data of M42 kindly provided by Harry. The left hand side shows the original, the right hand side the processed version. No processing beyond this and STF. The processes image  shows details of a kind I have not seen before - all kinds of waves, streams, structures, clounds.  No ringing (Gibbs effect). Stars are almost perfectly suppressed, only the delicate structures within the nebula are visible. HDR Wavelets, wavelets etc. all enhance a different kind of detail. For estethic images, it is probably necessary to combine the obtained image with some other data.
 
I must stress the whole implementation is still experimental (for example, it is not yet clear how to use the knobs for best effect, what images work and which don’t, and it is too slow/memory consuming for images beyond 1000x1000). If you want to test the process:
- I did my work on a virtual 1GB Fedora 14 Linux machine. For anything beyond 800x600, memory becomes a problem. On larger memory system you should be able to do larger images - if you have the time.
- compute times range from seconds to hours. There is some output at stdout that allows to monitor progress. Canceling the process does not work yet.
- I strongly recommend to try small snippets (100x100 or 200x200) first. For them, runtime is a few seconds.
- I did not yet try it on color images. For these, it works on luminance.

The next  message contains the module, sources and some more technical comments.

Georg

[georg.viehoever]

Here are the sources:

- The module has been tested and compiled on Fedora14-x64 with PI1.7.0 only.
- There is no reason why it should not work on the other platforms supported by PI, including Windows. In addtion to PCL you need gmm+ (http://download.gna.org/getfem/html/homepage/gmm.html), a library for vector/matrix operations (header only), On Fedora, it can easily be installed using "yum install gmm".
- The core of the operations happens in gradients.h where you also find additional information. Most other files are boilerplate code required by PI.
- I am not really happy with the speed of the linear solver. If you are more educated in sparse matrix operations than I am, and know how to quickly solve Poisson Equation Type PDE, you are welcome to join the project. I am currently using the Conjugate Gradients implementation of gmm+, and it is simply too slow for certain types of algorithms that I have in mind.
- Consider the module to be early beta. It most likely has many quirks.

If you need the module compiled for Fedora14-x64: Let me know by private message. It is too large to be attached here,

Enjoy,
Georg

[Carlos Milovic]

Hi Georg
I've been working in the same problem for a very long time. I wrote the whole structure for a Gradient Domain HDR Compression algoirthm, but my problem always was how to solve Poisson's equation. I tested several approaches, but neither of them worked quite well. One of the fastest was a Cosine transform based solver, with a Wiener-like deconvolution. Also I tried multigrid methods, but they did not work well (most probably, my fault). I wanted to implement the optimized multigrid solver you used for a previous script, but the code is way beyond my capabilities. I was hoping Juan would spare some time with me solving this... but it has been a very long wait.
If you don't mind, I'm gonna "steal" your module, and add the processing flow I designed, and integrate it to yours. I'll upload the new code ASAP. Then, I'll try to help you with the solver optimization...

[Harry page]

Hi

You are all very clever blokes  8)

Harry

[georg.viehoever]

Hi Carlos,

interesting to hear what you did. Regarding the solver:
- one of my problems is that for the kind of work I did with this module, I needed a solver that was fairly accurate. The multigrid solver I found for MergeMosaic was accurate enough for 8 bit pictures, but not the float kind I needed here.
- I understand that a Cosine transform based solver would only work for rectangular grids. I may be wrong about it ... its way beyond my mathematics. For the kind of processes I am currently thinking about, I would also need non-rectangular boundaries.
- Feel free to steal what you need 

Georg

[Carlos Milovic]

- one of my problems is that for the kind of work I did with this module, I needed a solver that was fairly accurate. The multigrid solver I found for MergeMosaic was accurate enough for 8 bit pictures, but not the float kind I needed here.

Same here... that's why I wanted to translate the code into the PCL, and did not just use it directly. But, as I said, is beyond my capabilities

- I understand that a Cosine transform based solver would only work for rectangular grids. I may be wrong about it ... its way beyond my mathematics. For the kind of processes I am currently thinking about, I would also need non-rectangular boundaries.

The whole point in the cosine transform, is that it satisfies the Neumann boundary conditions. Other boundary contitions may just use the FFT...
Is there any reason not to use rectangular boundaries? You may expand the "domain", inserting virtual data, or unmodified one...

[georg.viehoever]

Hi Carlos,

just have a look at the code in gradients.cpp/h, and see if you get this translated into DCT/FFT.
- the process at hand uses Von Neumann boundary condtions.
- I had the intention to provide something like a clone stamp tool. This would require Dirichlet boundary conditions.
- I am not sure: If I want to do something like MergeMosaic as a module, the background (value 0) of the merged images is a "dont care" region. The foreground will be a set of tilted rectangles. Dont know how this influences the optimal solver.

Georg

[Carlos Milovic]

If I am not wrong, dirichlet can be implemented with the sine transform... anyway, we have to go through the FFT to calculate it.

I though that the solver processed the whole image, in mergemosaic... I'll see it carefully (not in the short term, sorry), and see if I have any idea.

[Juan Conejero]

Hi Georg,

Fantastic work. This is an exciting new project that may give raise to extremely interesting tools in PixInsight. I'm afraid I can't help you too much in the mathematical part but count on me for the rest, including compilation/adaptation to all platforms and revision of PCL based code. A first suggestion would be integrating GMM header files (very nice library by the way, I didn't know it) with the module's source code. In this way you'd have a completely independent module. I have read the license terms and there is no problem in doing that, especially given that you are releasing your module as an open source product.

Keep it up!

[Juan Conejero]

By the way, I suppose you already know this but in case you don't:

http://www.cs.huji.ac.il/~danix/hdr/

[Carlos Milovic]

Georg, (et al)

Here is the code of my module. It is highly commented, since I was trying to just solve the poisson problem. The GDHDRC part (described in the paper Juan mentioned in that link) is already implemented, and so far looks to work fine (the mapping function seems ok... I have not tested the color correction).

As you will see, I split the code, taking off the solver part. I had implemented two options, that I'll describe briefly:

- Constrained Least Squares: this is a very naive implementation, and I hoped to make it work ... it is close, but not enough. Basically, I was hoping that a deconvolution would yield the result. Maybe with some more changes to the code it may do the job, at least as an approximation. Right now, it is generating artificial gradients that I'm not quite sure yet how to remove.

- Multigrid Algorithm: Once I made it work pretty close to the result... I modified some of the code, and now it doesn't I'm going to look for the original source that I tried to "translate" into the PCL.

[georg.viehoever]

Hi Carlos,

maybe you find this one useful, at least for the vonNeumann type. It appears to be fast, accurate, and pretty straight forward for someone trained in FFT: http://www.cs.cmu.edu/~jmccann/single-fft.pdf

Georg

[Carlos Milovic]

Thanks!
BTW, I read your code. Why do you leave PCL so soon? I think that the only problematic part is the poisson solving... so you may just go to gmm for that procedure. My approach is to use a header file with a "SolvePoisson( Image laplacian, etc etc)" function, that does this job You did not make things easy for jumping in

[georg.viehoever]

...BTW, I read your code. Why do you leave PCL so soon? ...

Well, every module has a history, and this one started as a command line tool, without all the PI boilerplate code (which -quite frankly- can be tedious to change). The easiest way to get it running in PI was to pass a PI image into the gmm structure - it has no significant performance penalty now, since most of the time is spent in the solver.

Georg

[Carlos Milovic]

I just read the paper. Looks quite straightforward. I think that I may implement it, with a little help of my friends

So, to begin with... The DCTII is the real part of the FFT, right? I'll have to dig a bit in the PCL's 1D FFT, and see how to implement it.


BTW, I liked the numerical accuracy of this solver I would choose that over the double speed of multigrid methods.