PSF fourier transforms and noise

[darkownt]

I posted a similar inquiry... probably in the wrong area.

Imagine the following:  you have a first specific PSF1 and a second also very specific (and different) PSF2.

Take Image0 and convolve it with PSF1 to get Image1 and convolve it with image0 to get Image2.

Add a certain moderate amount of noise of the same type and magnitude to both image1 and image2.


If you were given Image1 or Image2 without knowing which one it was or even that it had undergone a convolution by a PSF you would likely apply normal "blind" noisereduction.

If you are however given Image1 or Image2 with the knowledge of the PSF1 or PSF2 having been used to convolve an earlier unknown image (image0) would you be able to take into account the form of PSF1 or PSF2 when applying noise reduction?  i.e. if the PSF attenuates certain frequecies or impose slope/gradient limits cannot this be taken into account to when reducing noise i.e. certain sizes or types of features are now more easily identified as noise ?   

Further can the noise be "better" removed from an image which has certain restrictions and conditions imposed by being a convolved version of an unknown, than noise being removed from a "complete" unknown image (which could simply be noise itself or look much like the noise which was added).

Any thoughts... especially if they can involve using the FFT capabilities of PixInsight would be greatly appreciated.

DO

[Carlos Milovic]

Do you just want to remove the noise, or recover the ideal image ("Image0")?

[darkownt]

If I can do one I would be able to do the other.   

Really I want to know if some information (true signal has undergone a convolution with a PSF) is better than NO information (true signal unknown) and how to use that information to advantage.

 

[Carlos Milovic]

If I can do one I would be able to do the other.

Well, I ask because when dealing with inverse problems (recover the original, ideal image) one usually deal with both issues at the same time, in a regularized deconvolution. Often this is done with an iterative approach, where you are using either a dynamic equation, or minimizing a functional (data fidelity plus regularization terms).
A very simple way to deal with the problem is the Tikhonov regularization ( http://en.wikipedia.org/wiki/Tikhonov_regularization ). This is also the same principle behind SVD regularization and the Wiener filter. In a nutshell, since noise lives in high frequency components, you amplify the low frequency components of the image until a cut off frequency that is related to the noise. This way, you enhance the image, without noise amplification. In the Wiener filter the transition is smooth and you may achieve some noise reduction.
Of course there are much better ways to deal with the problem. Since the Tikhonov regularization searches for an optimal solution that is closer to the ideal image in a least squares sense, loss of definition in the edges and other artifacts are generated (ringing, due the high degree continuity that it imposes). Total variation based regularizations, for example, assume that the image is piecewise constant, thus achieving sharp edges, at the expense of generating stairlike gradients, or discontinuities in smooth transitions. The GRAYCstoration is a complex example of this kind of regularization, used for noise reduction purposes (also, in the page of the authors of that algorithm, you may see other uses, like inpainting).
There is also a more complex regularization, that we are working on, that is called Total Generalized Variation. In this algorithm, the image is assumed to be piecewise smooth, avoiding the stairlike artifacts and most of the ringing.

So, if you want to get into this field, I recommend you to do some research on inverse problems, and regularization methods. Also, you may look into representation in other spaces (like wavelets), that may be more effective to isolate noise than the Fourier transform. Usually the DFT is a great tool to simplify the convolution processes, but is not that great to reduce noise o model it, unless it is a clear periodical signal.

[darkownt]

I have come to similar conclusions regarding the use of FFT.  General non periodic noise is not easy to deal with although some of my experiments with gating based on the PSF do reduce noise some artifacts creep in.  I am thinking of working in the spatial domain...

As for research... this forum is my first step in research... and the concept, the idea, that I am trying to squeeze... is that there is a fundamental difference between processing noise in a vacum (i.e. absent any further information) and processing noise on an image when you know the true image has undergone a particular convolution.  To my mind, although I do not know the techniques, it seems that the particular imformation (the PSF) should be usable SOMEhow.  I'm no genius someone else must have realized this... I am wondering HOW they used it... i.e. how CAN I use it?   

Clearly an image after going through a gaussian filter (rad >1) would not naturally have salt and pepper noise.. or noise that has spacial slopes greater than that possible after filtering with the specific gaussian, therefore a smart enough algorithm should be able to somehow understand this.

I like the idea of wavelets... can they be defined/based on a specific form such as a PSF profile? taking the coefficient (strength) of the PSF at each radius (circularly symmetric) how can we use wavelets?  Can I decompose an image arbitrarily in PixInsight in some basis determined by components of a PSF?

[Carlos Milovic]

Instead of assuming that the image has undergo a certain PSF as prior information, just for noise reduction information, it is more useful to think on other kind a prior informations, or "guesses". For example, in the Compressed Sense field, the image is assumed to be sparse (only a few pixels with values different to zero). Then, they minimize the L0-norm of the image, or the L1-norm.

The thing is that if you don't have a cutoff frequency, you may end deleting true signal at high frequencies that have a low amplitude in the Fourier domain, due the PSF. The problem is how to determine what is signal, and what is noise, and I'm afraid that this kind of analysis in the Fourier domain is quite complex (you may see the optimized or adaptive notch filter, wich tries to atenuate certain frequencies until they match the texture of the image, minimizing variation). Wavelet domain techniques may be easier to perform, following the line of the soft and hard thresholds. The scaling function may be whatever convolution you want, and certainly you may use the PSF. I'm not sure if this is the best approach, but indeed it should "work". I recommend you the book: Sparse Image and Signal Processing (Wavelets, Curvelets, Morphological Diversity) from Jean-Luc Starck, Fionn Murtagh and Jalal M. Fadili.

[darkownt]

Hello

Thanks for your reply.  I have been using the scaling function (kernal) option in the atrouswavlet function (script) and have looked at the results when it matches a PSF previously applied to an image and I am having a bit of a difficult time trying to interpret the significance of the results.

I did some reading, and it seems that the scaling function and the wavlet function are related by W(x) = 2*S(2x) - S(x), or in other words the wavelet function is a combination of a narrow tall version of the scaling function minus a normal height twice wider version of the scaling function.  I am wondering if I should set my scaling function kernal to my PSF 2x larger?

What the actual significance of applying this is still something I need to consider...

cheers

[Carlos Milovic]

You may better understand the wavelet layers if you take a look at the gaussian and laplacian pyramids. In the ATrous Transform, instead of scaling down each layer, the filter is sort of scaled up.

[darkownt]

OK so I am using the

image.aTrousWaveletTransform

function to extract layers using an arbitrary scaling function (a grey iiimage of the PSF)

I am trying to interpret the results,  is there anything other than the layers that can be gotten from the function?

For example can I get access to the coefficients of the wavlets used to generate each layer rather than each layer itself?

[Carlos Milovic]

The ATrousWavelets Transform, or Starlets, is a redundant algorithm. This means, that each layer containing the wavelet coefficients has the same dimensions as the original image. So, if you want the coefficients, just read the pixel values (and also, this is why the inverse starlet transform is so easy: just add all the coefficients for a given pixel, at each layer plus the remanent).

[darkownt]

Hi Carlos:

When I do a wavelet transform of an image using the ATrouswavelettransform script, I get a number of layers which are rescaled versions of the extracted layers.   I understand each layer is a superposition of wavelets of the particular scale (each wavelet being related to the scaling function as described in my previous email).  We know this since each of these layers so decomposed can be added to acheive the entire image.

What I mean by coefficients are the factors (along with positions) which determine the magnitude and position of each wavelet used for each said superposition.

Each layer is made up of wavelet functions of the same size, having specific "centers" and magnitudes.  More simply, each pixel (position) can be understood to have a factor associated with the magnitude of the wavelet (in RGB space) of the particular size for the layer centered on that pixel.  Since the wavelets are larger than a single pixel, the superposition of all of them (i.e. the final rendered individual layer) does not just consist of the coefficients (factors).

I just had a thought.  If I generate the wavelets manually as a deconvolution kernal and use it to deconvolve each appropriate wavelet layer I should be able to in theory recover the coefficients...

I think I may have a solution to my problem...

[Carlos Milovic]

If you are using the preview option, or the script, the layers are rescaled in the sense of using the whole normalized range (i.e., minimum value is assigned to zero, and the larger to one). True wavelet layers, as you can get when using floating point images and the PCL or PJSR classes, have negative and positive values, usually in the [-1,1] range. Each layer is directly the wavelet coefficients. No combination. The only layer that is a linear combination of layers is the residual one. In the ATWT, all layers have the same dimensions as the original image. What changes in size is the scaling function, the convolution operation. That's why I suggested to you to take a look at the Laplacian pyramid. The ATWT builds its layers in the same way, but upscaling the function rather than downsampling each layer.

[darkownt]

Good.  Thats good to know.

I have also been thinking about the form of the scaling function I am using.  Since the Wavelet function W is related to the Scaling function

W(x) proportional to S(2x) - 1/2 S(x)

I am thinking of using a scaling function which is related to the PSF P(x) by the following:

S(x) =   P(2x) + 1/2 P(x) + 1/4 P(x/2) + 1/8 P(x/4) + ...

This leads to a wavelet function which is proportional to the PSF but with a tiny magnitude very spread out "left over" portion.

I was trying to build an image S from a PSF P by multiplying with pixel math (no problems 0.5, 0.25, etc.) and scaling with the resampling tool and cropping out the center. Unfortunately resampling causes the center to drift substantially, is there a resampling method that is smooth and accurate to create something like P(x/8)  (i.e. scaling of 800%) from P(x)?

For now I will try integer resampling but that is not smooth... 

[EDIT:  Corrected the equations...]