Changeset 20285

Timestamp:

Oct 20, 2008, 11:29:18 PM (18 years ago)

Author:

Sebastian Jester

Message:

Working version of Gauss-fitting functions

File:

• branches/sj_ippTests_branch_20080929/ippTests/compIPPphoto.py (modified) (6 diffs)

branches/sj_ippTests_branch_20080929/ippTests/compIPPphoto.py

@@ -759 +759 @@
     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):
@@ -770 +769 @@
             inradius = (((x-x0)**2 + (y-y0)**2) > blankRadius**2)
         inradius = reshape(inradius,(inradius.shape[0]*inradius.shape[1],))
-        resid[inradius==1] = 0
+        resid[inradius] = 0
+    # win.showArray(reshape(resid,model.shape))
     return resid
     
@@ -794 +794 @@
         theta -= 180
     bestp=(scale,x0,xsig,y0,ysig,theta)
-    resid=reshape(gaussResid(bestp,yx,data,blankrad),data.shape)
+    resid=reshape(gaussResid(bestp,yx,data,blankrad,blankInner),data.shape)
     mod=data-resid
     return bestp,resid,mod
@@ -812 +812 @@
     
 def doGaussPsfPhotometry(infile,initialGuessList=[(2000.,1011.,2.,1176.,2.,0)],\
-                             blankrad=15,defHalfwinsize=15,HDU=0):
+                             blankrad=9,defHalfwinsize=15,HDU=0):
     """Fit Gaussians to image in infile's HDU near the
     positions and with fluxes and sigmas given in initialGuessList"""
+    # XXX Todo: Make frame with residuals. Make blank array, then insert residuals.    
     import pyfits
     from numpy import nan
+    from numarray import ones,zeros,Float32
+    from scipy import ogrid
     f_handle=pyfits.open(infile)
     img = f_handle[HDU].data
+    resid_img = 1e-7*ones(img.shape,Float32)
+    print resid_img.shape,type(resid_img),img.shape,type(img)
     fitlist = []
     for initialGuessTuple in initialGuessList:
@@ -829 +834 @@
             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]
+        # ds9-based coords start at 1,1, while numpy's start at 0,0,
+        # so we need to add 1 to all internally used coords to get the
+        # fitted coordinates to match with what we see in ds9
+        yx=ogrid[1+int(y0)-halfwinsize:1+int(y0)+halfwinsize+1,1+int(x0)-halfwinsize:1+int(x0)+halfwinsize+1]
         try:
             par,res,bestfit=fitGauss(yx,subim,initialGuessTuple,blankrad,blankInner=False)
         except:
             par=(nan,nan,nan,nan,nan,nan)
+            res = subim
         (scale,x,xsig,y,ysig,theta)=par
+        win.showArray(res)
+        resid_img[int(y0)-halfwinsize:int(y0)+halfwinsize+1,int(x0)-halfwinsize:int(x0)+halfwinsize+1] = res
         print initialGuessTuple, scale,x,xsig,y,ysig,theta
         fitlist.append(par)
     f_handle.close()
-    return par
+    outname = 'gaussfit_resid_'+infile
+    primHDU = pyfits.PrimaryHDU(resid_img)
+    hdulist = pyfits.HDUList([primHDU])
+    hdulist.writeto(outname,clobber=True)
+    return fitlist
         
 def doEfbComparison(simfile='simtest.1.00a.fits',efbfile='efbsim.1.00a.fits',simtab='simtest.1.00a.cmf'):
     import pyfits
-    from numpy import array
+    from numpy import array, rec
+    import re
     tab = pyfits.open(simtab)
     tdata = tab[1].data
@@ -861 +874 @@
 
     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'])
+    saverows = 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)
-
-    
+    fitsend = re.compile('(\.[sc]mf|\.fits?)$')
+    saveHDU.writeto('gaussfit_'+fitsend.sub('',efbfile)+'___'+simfile,clobber=True)
+    
+
+