[SciPy-User] Fitting a model to an image to recover the position

Niemi, Sami sniemi@email.unc....
Fri Sep 30 10:18:33 CDT 2011


Hello,

I am trying to solve a problem where I need to fit a model to an image to recover the position (of the model in that image). The real-life application is more complicated (fitting sparse field spectroscopic data to an SDSS r-band image, if you are interested in) than the simple example I give below, but it is analogous. The most significant difference being that in the real-life application I need to allow rotations (so I need to find a position x and y and rotation r that minimizes for example the chi**2 value) and that the difference between the model and image might be larger than the small random error applied in the simple example (and that there is less information in one of the directions because of the finite slit width).

The simple example I show below works for the real-life problem, but it's far from being effective or elegant. I would like to use some in-built minimization methods like scipy.optimize.leastsq to solve this problem but I haven't found a way to get it work on my problem. Any ideas how I could do this better?


Cheers,
Sami



Simple Example:

import numpy as np

def findPosition(image, model, xguess, yguess, length, width, xrange=20, yrange=20):
    '''                                                                                                                                                                           
    Finds the x and y position of a model in an image that minimizes chi**2 by looping                                                                                            
    a grid around the initial guess given by xguess and yguess.                                                                                                                   
                                                                                                                                                                                  
    This method is far from optimal, so the question is how to do this                                                                                                            
    with the scipy.optimize.leastsq or some other built-in min. algorithm.                                                                                                        
    '''
    out = []
    for x in range(-xrange, xrange):
        for y in range(-yrange, yrange):
            obs = image[yguess+y:yguess+y+length, xguess+x:xguess+x+width].ravel()
            chi2 = np.sum((obs - model)**2 / model)
            out.append([chi2, x, y])
    out = np.array(out)

    return out, np.argmin(out[:,0])


#create an image of 100 times 100 pixels of random data as an example, this represents the imaging data                                                                           
image = np.random.rand(10000).reshape(100,100) + 1000

#take a slice and add a small error, this represents the model that had been recovered somewhere else                                                                             
xpos, ypos = 80, 20
length, width = 21, 4
model = (image[ypos:ypos+length, xpos:xpos+width].copy() + np.random.rand(1)[0]*0.1).ravel()

#initial guess of the position (this is close, the correct one is 80, 20)                                                                                                         
xguess, yguess = 75, 33
#find the position using the idiotic algorithm                                                                                                                                    
out, armin = findPosition(image, model, xguess, yguess, length, width)

#print out the recovered shift and position                                                                                                                                       
print 'xshift yshift chi**2'
print out[armin,1], out[armin, 2], out[armin, 0]
print 'xfinal yfinal'
print xguess+out[armin,1], yguess+out[armin, 2], 'which should be %i, %i' % (xpos, ypos)

 


More information about the SciPy-User mailing list