[SciPy-user] [noob inside]star detection

[Gael Varoquaux]

On Thu, May 17, 2007 at 02:44:19PM +0200, Jean-Baptiste BUTET wrote:
> I have some ideas how to begin this... but I think there are function
> in scipy that are dedicated to this. (scipy or numpy I use both)

> 1) find local minima
> 2) fit a gaussian/PSF on it
> 3) compute fwhm

> 1) or 2) are hot for me.

I think there is no need for fit when dealing with gaussian data: all you
need to do is find the moments of the distribution.

I will be lazy and steal some code I have already used to explain this,
you can find it on:
http://scipy.org/Cookbook/FittingData

OK, this code is 1D, and you want 2D code. I'll try not to be to lazy and
guive you some 2D code:

+++++++++++++++++++
from pylab import *
from numpy import *

# Create the gaussian data
gaussian = lambda x, y: 3*exp(-((30-x)**2+(20+y)**2)/20.)

Xin, Yin = mgrid[-50:51, -50:51]
data = gaussian(Xin, Yin)

matshow(data, cmap=cm.gist_earth_r)

# Now pretend we do not know how this data was created, and find the
# parameters.
total = sum(data)
X, Y = indices(data.shape)
x = sum(X*data)/total
y = sum(Y*data)/total
width_x = sqrt(abs(sum((X-x)**2*data)/total))
width_x = sqrt(abs(sum((Y-y)**2*data)/total))

max = data.max()

fit = lambda u,v : max*exp(-((u-x)**2+(v-y)**2)/(2*width**2))

matshow(fit(X, Y), cmap=cm.gist_earth_r)

show()
+++++++++++++++++++

I hope this code is understandable to you and that you elaborate on it.
If not, ask questions here.

-- 
    Gaël