comp.lang.idl-pvwave archive: archive

Home » Public Forums » archive » simple deconvolution

Show: Today's Messages :: Show Polls :: Message Navigator
E-mail to friend

Re: simple deconvolution [message #75287 is a reply to message #75210]

Wed, 23 February 2011 08:14

pgrigis
Messages: 436
Registered: September 2007

Senior Member

On Feb 22, 10:00 am, chris <rog...@googlemail.com> wrote:
> Hi folks,
> I want to implement an image deconvolution into a larger package. The
> following code performs either the Iterative Wiener (by A.W.
> Stevenson) or the Richardson-Lucy deconvolution, but both go wrong for
> the recovery of both smoothed images and smoothed images plus noise .
> I'm a little bit confused about that. Maybe somebody could help me?
> The implemented CONVOLVE comes from the Astrolib. I'm using IDL 8 and
> the code is not optimised as you can see :)
>
> function cr_deconv,im,psf,method,small=small
> sz1 = size(im,/dimensions)
> sz2 = size(psf,/dimensions)
> small=~n_elements(small)?1e-5:small
> if total(sz1 eq sz2) ne 0 then begin
> p=fltarr(sz1)
> p[(sz1[0]/2)-(sz2[0]/2) ,(sz1[1]/2)-(sz2[1]/2)]=psf
> endif
> p/=total(psf)
> p[where(p lt small)]=small
> if method eq 'iwiener' then begin
> psf_fft=fft(p)
> psf_fft[where(abs(psf_fft) lt small)]=small
> snr=mean(median(im,3))/stddev(im-median(im,3)) : snr
> pc=psf_fft*conj(p)
> pc[where(abs(pc) lt small)]=small
> filter=pc
> filter/=(filter+1./snr)
> filter[where(abs(filter) lt small)]=small
> res=abs(fft(filter*fft(im)/psf_fft,/inverse))
> for i=0l,iter-1l do begin
> res+=abs(fft((fft(convolve(i eq 0?im:res,p)-im)/psf_fft)*$
> (pc/(pc+(1./snr))),/inverse))
> snr=mean(median(res,3))/stddev(res-median(res,3))
> endfor
> else begin
> corr_kernel=rot(p,180)
> for i=0l,iter-1l do $
> res=(i eq 0?im:res)*convolve(im/convolve(i eq 0?
> im:res,p),corr_kernel)
> endelse
> return,res
> end
>
> Thanks in advance
>
> CR

My understanding is that the Richardson-Lucy algorithm
works as follows.

Given an Image IM and a point-spread function PSF.

Initialization:
O=IM

Loop:
IHAT=CONV(PSF,O)
O=O*CORR(IM/IHAT,PSF)

After somewhere between 10 to 50 iterations, O is going to
be an approximation to the the deconvolved version of IM.

Here CONV and CORR are the usual convolution and correlation
functions. Some care need to be taken with normalization, but
this is the skeleton of the algorithm.

I do not see that your algorithm is performing this operation,
or is it? Also you may want to implement the convolutions and
correlations manually yourself using FFT - this way you have
more control over what is happening.

Ciao,
Paolo

Report message to a moderator

[Message index]

		simple deconvolution By: rogass on Tue, 22 February 2011 07:00
		Re: simple deconvolution By: rogass on Wed, 23 February 2011 12:32
		Re: simple deconvolution By: pgrigis on Wed, 23 February 2011 08:14
		Re: simple deconvolution By: penteado on Wed, 23 February 2011 07:45
		Re: simple deconvolution By: David Fanning on Wed, 23 February 2011 07:30
		Re: simple deconvolution By: Jeremy Bailin on Wed, 23 February 2011 06:53
		Re: simple deconvolution By: rogass on Tue, 22 February 2011 23:40
		Re: simple deconvolution By: James[2] on Tue, 22 February 2011 15:49
		Re: simple deconvolution By: David Fanning on Tue, 22 February 2011 15:24
		Re: simple deconvolution By: James[2] on Tue, 22 February 2011 15:15
		Re: simple deconvolution By: wlandsman on Tue, 22 February 2011 15:03
		Re: simple deconvolution By: David Fanning on Tue, 22 February 2011 09:10
		Re: simple deconvolution By: rogass on Tue, 22 February 2011 08:49
		Re: simple deconvolution By: penteado on Tue, 22 February 2011 08:35

Previous Topic:	Re: Pixelwise temporal trend, problem with REFORM
Next Topic:	IDL on win7 crashing

-=] Back to Top [=-

[ Syndicate this forum (XML) ] [

] [

]

Current Time: Wed Oct 08 19:20:42 PDT 2025

Total time taken to generate the page: 0.00234 seconds