Index: /branches/sj_ippTests_branch_20080929/ippTests/compIPPphoto.py
===================================================================
--- /branches/sj_ippTests_branch_20080929/ippTests/compIPPphoto.py	(revision 20278)
+++ /branches/sj_ippTests_branch_20080929/ippTests/compIPPphoto.py	(revision 20279)
@@ -751,3 +751,124 @@
         box()
 
-    
+def gaussResid(p,yx,data,blankRadius=0,blankInner=True):
+    """return residuals after fitting two-D gaussian with parameters p
+    = scale,x0,xsig, y0, ysig, theta, blanking pixels inside or
+    outside blankRadius, depending on whether blankInner = True or
+    False"""
+    from numpy import reshape
+    scale,x0,xsig,y0,ysig,theta = p
+    model=scale*gauss2D(yx,x0,xsig,y0,ysig,theta)
+    win.showArray(data-model)
+    resid = reshape(data,(data.shape[0]*data.shape[1],))-reshape(model,(model.shape[0]*model.shape[1],))
+    if (blankRadius > 0):
+	# set residuals within or outside blankRadius to 0, i.e. set model = data there
+	x=yx[1]
+	y=yx[0]
+        if blankInner:
+            inradius = (((x-x0)**2 + (y-y0)**2) <= blankRadius**2)
+        else:
+            inradius = (((x-x0)**2 + (y-y0)**2) > blankRadius**2)
+        inradius = reshape(inradius,(inradius.shape[0]*inradius.shape[1],))
+	resid[inradius==1] = 0
+    return resid
+    
+def fitGauss(yx,data,p0,blankrad,blankInner):
+    """fit gaussian to data with starting parameters in p0"""
+    from scipy.optimize import leastsq
+    from numpy import reshape
+    # need to add weighting by errors, if we know them
+    plsq = leastsq(gaussResid,p0,args=(yx,data,blankrad,blankInner))
+    bestp = plsq[0]
+    # print plsq[1]
+    scale,x0,xsig,y0,ysig,theta = bestp
+    # remove 360 degree ambivalence
+    theta -= int(theta/180)*180
+    # make sure xsig is major axis
+    if xsig < ysig:
+	tmp=xsig
+	xsig=ysig
+	ysig=tmp
+	theta += 90
+    # make sure theta is in [-90,90]
+    if theta > 90:
+	theta -= 180
+    bestp=(scale,x0,xsig,y0,ysig,theta)
+    resid=reshape(gaussResid(bestp,yx,data,blankrad),data.shape)
+    mod=data-resid
+    return bestp,resid,mod
+
+def gauss2D(coordgrid,x0,xsig,y0,ysig,theta=0):
+    """compute a normalized gaussian at positions given in 2d-array
+    coords, centered at x0, y0, with xsig, ysig for theta=0;
+    coordinates in coordgrid given by ogrid"""
+    from numpy import cos,sin,pi,exp
+    cth = cos(theta/180.*pi)
+    sth = sin(theta/180.*pi)
+    
+    x=coordgrid[1]-x0
+    y=coordgrid[0]-y0
+    res = exp(-0.5*(((x*cth+y*sth)/xsig)**2+((-x*sth+y*cth)/ysig)**2))
+    return res/(2*pi*xsig*ysig)    
+    
+def doGaussPsfPhotometry(infile,initialGuessList=[(2000.,1011.,2.,1176.,2.,0)],\
+                             blankrad=15,defHalfwinsize=15,HDU=0):
+    """Fit Gaussians to image in infile's HDU near the
+    positions and with fluxes and sigmas given in initialGuessList"""
+    import pyfits
+    from numpy import nan
+    f_handle=pyfits.open(infile)
+    img = f_handle[HDU].data
+    fitlist = []
+    for initialGuessTuple in initialGuessList:
+        scale0,x0,xsig0,y0,ysig0,theta0 = initialGuessTuple
+        # Allow for variable blankrad leading to variable halfwinsize
+        # (no information beyond blankrad will be used), but avoid going to 0
+        if blankrad > 0:
+            halfwinsize = blankrad
+        else:
+            halfwinsize = defHalfwinsize
+        subim = img[int(y0)-halfwinsize:int(y0)+halfwinsize+1,int(x0)-halfwinsize:int(x0)+halfwinsize+1]
+        # ds9-based coords start at 1,1, while these start at 0,0, so
+        # we need to add 1 to all internally used coords to get them
+        # to match with what we see in ds9
+        yx=scipy.ogrid[1+y0-halfwinsize:1+y0+halfwinsize+1,1+x0-halfwinsize:1+x0+halfwinsize+1]
+        try:
+            par,res,bestfit=fitGauss(yx,subim,initialGuessTuple,blankrad,blankInner=False)
+        except:
+            par=(nan,nan,nan,nan,nan,nan)
+        (scale,x,xsig,y,ysig,theta)=par
+        print initialGuessTuple, scale,x,xsig,y,ysig,theta
+        fitlist.append(par)
+    f_handle.close()
+    return par
+        
+def doEfbComparison(simfile='simtest.1.00a.fits',efbfile='efbsim.1.00a.fits',simtab='simtest.1.00a.cmf'):
+    import pyfits
+    from numpy import array
+    tab = pyfits.open(simtab)
+    tdata = tab[1].data
+    x = array(tdata.field('x_psf'))
+    y = array(tdata.field('y_psf'))
+    counts = 10**(-0.4*array(tdata.field('psf_inst_mag')))
+    initlist = []
+    inparlist = []
+    nrows = len(tdata.field('IPP_IDET'))
+    for i in range(nrows):
+        initlist.append((counts[i],x[i],2.,y[i],2.,0.))
+        inparlist.append((counts[i],x[i],y[i]))
+    efbparlist = doGaussPsfPhotometry(efbfile,initlist)
+    simparlist = doGaussPsfPhotometry(simfile,initlist)
+
+    savelist = [ inparlist[i] + efbparlist[i] + simparlist[i] for i in range(nrows)]
+    saverows = numpy.rec.array(savelist,names=['counts_in','x_in','y_in',\
+                                                   'efbfit_counts','efbfit_x','efbfit_fwhm_x',\
+                                                   'efbfit_y','efbfit_fwhm_y','efbfit_theta',\
+                                                   'simfit_counts','simfit_x','simfit_fwhm_x',\
+                                                   'simfit_y','simfit_fwhm_y','simfit_theta'])
+    for col in ['efbfit_fwhm_x','efbfit_fwhm_y','simfit_fwhm_x','simfit_fwhm_y']:
+        saverows.field(col)[0:] *= 2.35
+    
+    saveHDU = tabHDUfromRecArray(saverows)
+    saveHDU.writeto('gaussfit_'+efbfile)
+
+