pyminuit icon indicating copy to clipboard operation
pyminuit copied to clipboard

Pyminuit Richardson lucy implementation

Open GoogleCodeExporter opened this issue 10 years ago • 4 comments

I am trying to use Pyminuit to iterate the richardson lucy(RL) algorithm as a 
possible method of image deconvolution. However I am getting a few errors 
through migrad(). I have tried two RL lambda functions, one using a summation 
method and the other a convolution method. Both can be found here, 
http://en.wikipedia.org/wiki/Richardson%E2%80%93Lucy_deconvolution. Here is my 
code snippet. 

    def Richardsonlucy_deconvolve_min(image_data='/Users/krishnamooroogen/Documents/physics/MSSL_projectfiles/VSOdata/AIA/2011_06_07/aia_lev1_211a_2011_06_07t06_27_36_63z_image_lev1_fits.fits' ,psf='/Users/krishnamooroogen/Documents/physics/MSSL_projectfiles/psf/psf_211_AIA.fits'):
    #opens image data and psf
    #creates inverse psf
    #sets lambda function (ricahrdson lucy)
    #sets start value
    #iterates gradient

        image_data = pyfits.open(image_data)[0].data
        psf = pyfits.open(psf)[0].data
        psf_hat = psf[::-1]

        #RL = lambda latent: latent*(convolute((image_data/(convolute(latent,psf))),psf_hat))
        RL = lambda latent: latent*(((image_data/(latent*psf).sum())*psf).sum())

        m = minuit.Minuit(RL,latent = numpy.array([numpy.array([1.0]*len(image_data)),numpy.array([1.0]*len(image_data))]*int(len(image_data)/2)))
        m.tol = 1000000000
        m.printMode = 1
        m.migrad()
        m.hesse()
        latent = m.values['latent']

The convolute function is here, 

    def convolute(input,input2):

    fourier_space_convolution = (numpy.fft.rfft2(input,input.shape))*(numpy.fft.rfft2(input2,input2.shape))
    convert_realspace_convolution = numpy.fft.fftshift(numpy.fft.irfft2(fourier_space_convolution,fourier_space_convolution.shape))
    return convert_realspace_convolution

When using the summation method of RL I get the following error, 

    --> 359     Richardsonlucy_deconvolve_min()

    /Users/krishnamooroogen/Documents/Modules/MSSL_solar/solar_usr.pyc in Richardsonlucy_deconvolve_min(image_data, psf)
    276         m.tol = 1000000000
    277         m.printMode = 1
    --> 278         m.migrad()
    279         m.hesse()
    280         latent = m.values['latent']

    /Users/krishnamooroogen/Documents/Modules/MSSL_solar/solar_usr.pyc in <lambda>(latent)
    271 
    272         #RL = lambda latent: latent*  (convolute((image_data/(convolute(latent,psf))),psf_hat))
    --> 273         RL = lambda latent: latent*(((image_data/(latent*psf).sum())*psf).sum())
    274 
    275         m = minuit.Minuit(RL,latent = numpy.array([numpy.array([1.0]*len(image_data)),numpy.array([1.0]*len(image_data))]*int(len(image_data)/2)))

    TypeError: only length-1 arrays can be converted to Python scalars 

The error when using the convolute method, 

    AttributeError: 'float' object has no attribute 'shape' 


However outside of minuit, both the convolute function and the RL summation 
method for the latent image estimation work fine. Can someone shed some light, 
is this a Minuit specific problem? Can it be resolved? 

Original issue reported on code.google.com by [email protected] on 21 Nov 2013 at 1:00

GoogleCodeExporter avatar Mar 15 '15 22:03 GoogleCodeExporter