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

Gael Varoquaux gael.varoquaux@normalesup....
Mon Oct 3 02:50:06 CDT 2011

It seems to me that you have what is known as a registration problem. It
has many local minima and is a nasty optimization problem. I am not an
expert in registration, but I suggest that you try a multi-scale approach
if your angles and translations are large: downsample by a large factor
the two images (target and data), fit these, and use the rotation and
translation parameters learned to initialize less downsampled fits.
Iterate a few times for different values of downsampling. This helps
avoiding the local minima.

All the registration strategies are very dependant on the kind of data
that you have. Registration is more of an art than a science.


On Fri, Sep 30, 2011 at 03:18:33PM +0000, Niemi, Sami wrote:
> 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)

> _______________________________________________
> SciPy-User mailing list
> SciPy-User@scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-user

    Gael Varoquaux
    Research Fellow, INSERM
    Associate researcher, INRIA
    Laboratoire de Neuro-Imagerie Assistee par Ordinateur
    NeuroSpin/CEA Saclay , Bat 145, 91191 Gif-sur-Yvette France
    Phone:  ++ 33-1-69-08-78-35
    Mobile: ++ 33-6-28-25-64-62

More information about the SciPy-User mailing list