monnier@physics2 (John D. Monnier) wrote:
>
>
> I am in need of routine to fit a 2-D Gaussian to a small set of image data.
> I would like the model gaussian fit to be a function of 6 parameters;
> Height, Major Axis sigma, Minor Axis sigma, Centroid X, Centroid Y,
> and Major Axis Angle. As you can see, this is a 6-parameter, non-linear
> model so the fitting procedure is a bit complicated.
I wrote the following using the NR_POWELL routine. I also wrote a
program using NR_DFPMIN (uses derivatives) and I'll post it later
as soon as I put some comments on it.
I also wrote a routine for simplex minimization. Maybe later on, I'll
send that too. It needs a bit more programming refinement.
Stefano
<--<>--<>--<>--<>--<>--<>--<>--<>--<>--<>--<>--<>-- <>--<>--<>--<>--<>-->
Stefano Polizzi |
Dipartimento di Chimica Fisica | tel: +39-41-5298618
Universita' di Venezia | FAX: +39-41-5298594
Dorsoduro 2945 | e-mail:polizzi@unive.it
I-30123 VENEZIA (Italy) |
FUNCTION gaussian,par
common dati,image,x,y,f0
buffer=size(image)
col=buffer(1) & row=buffer(2)
if (col lt row)then maxr=row else maxr=col
;make an square array for your x and for y
xx=replicate(1.,maxr)#(findgen(maxr)-par(1))
yy=(findgen(maxr)-par(2))# replicate(1.,maxr)
;rotate it
x=(cos(par(5))*xx-sin(par(5))*yy)/par(3)
y=(sin(par(5))*xx+cos(par(5))*yy)/par(4)
f0=par(0)*exp(-(x^2+y^2)) ;calculate z (gaussian)
f0=f0(0:col-1,0:row-1) ;cut it for your rectangular image
f=total((f0-image)^2) ;sum of residuals square (you can weigh it, if you need)
print,f
return,f
END
pro powell
common dati,image,x,y,f0
col=61 &row=69 ;your image dimension
filename='bcampio.bin' ;your image
image=findgen(col,row)
openr,unit,filename,/get_lun
readu,unit,image ;read them
free_lun,unit
Gtol=1.0e-2 ;convergence parameter
;initial parameters (height,sigma-x,sigma-y,xzero-position,yzero-position,angle)
par=[max(image),30.,33.,4.,33.,-50*!DtoR]
buffer=size(par)
n=buffer(1)
xi=fltarr(n,n)
for i=0,n-1 do begin
xi(i,i) =1.0 ;initial direction vectors
endfor
;the following to see your initial guess
;f=gaussian(par)
;window,/free,xsize=col,ysize=row
;tvscl,image
;window,/free,xsize=col,ysize=row
;tvscl,f0
;window,/free,xsize=500,ysize=350
;surface,image,zrange=[0,max(image)] & surface,f0,zrange=[0,max(image)],color=12,/noerase
;stop
NR_POWELL,par,Xi,Gtol,Fmin,'gaussian',/double,iter=niter ;minimize
print,niter,par ,par(5)*!RaDeg,fmin
window,/free,xsize=col,ysize=row
tvscl,f0
window,/free,xsize=500,ysize=350
surface,image,zrange=[0,max(image)] & surface,f0,zrange=[0,max(image)],color=12,/noerase
;stop
END
|
comp.lang.idl-pvwave archive
Members
Search
Help
Login
Home
