|
casa
$Rev:20696$
|
|
casa
$Rev:20696$
|
00001 from casac import casac 00002 from tw_utils import * 00003 import os 00004 from time import * 00005 import numpy as numarray 00006 import pylab 00007 import string 00008 00009 #from math import * 00010 #from Scientific.Functions import LeastSquares 00011 import tw_func 00012 00013 00014 class ImageTest: 00015 def __init__(self, imageName, write=True, 00016 resultDir='WEBPAGES/imageTest/', 00017 imDir='IMAGES/'): 00018 00019 import shutil 00020 self.imTool=casac.image() 00021 self.imTool.open(imageName) #open('tables/squint_corr.restored') 00022 00023 # fix up images that don't have CASA-canonical stokes and spec: 00024 # assume for now that direction is in 01 at least, 00025 mycs=self.imTool.coordsys() 00026 findstok=mycs.findcoordinate("stokes") 00027 if not findstok[0]: 00028 myImagename=imageName+".k" 00029 self.imTool.adddegaxes(stokes=True,outfile=myImagename,overwrite=True) 00030 mystokpix=self.imTool.summary()['ndim'] # ct from 0 00031 self.imTool.close() 00032 shutil.rmtree(imageName) 00033 shutil.move(myImagename,imageName) 00034 self.imTool.open(imageName) 00035 mycs.done() 00036 mycs=self.imTool.coordsys() 00037 else: 00038 mystokpix=findstok[1] 00039 00040 findspec=mycs.findcoordinate("spectral") 00041 if not findspec[0]: 00042 myImagename=imageName+".s" 00043 self.imTool.adddegaxes(spectral=True,outfile=myImagename,overwrite=True) 00044 myspecpix=self.imTool.summary()['ndim'] # ct from 0 00045 self.imTool.close() 00046 shutil.rmtree(imageName) 00047 shutil.move(myImagename,imageName) 00048 self.imTool.open(imageName) 00049 mycs.done() 00050 mycs=self.imTool.coordsys() 00051 else: 00052 myspecpix=findspec[1] 00053 00054 curr_order=[mystokpix,myspecpix] 00055 if curr_order != [2,3]: 00056 myImagename=imageName+".t" 00057 self.imTool.transpose(order="01%1i%1i" % (mystokpix,myspecpix),outfile=myImagename,overwrite=True) 00058 shutil.rmtree(imageName) 00059 shutil.move(myImagename,imageName) 00060 self.imTool.open(imageName) 00061 00062 self.rgTool=casac.regionmanager() 00063 self.qaTool=casac.quanta() 00064 self.getArr() #instead make explicit call to getArr() 00065 self.write=write 00066 self.imageName=imageName 00067 self.iterate=0 #for multiple cubeFit() tests 00068 00069 if self.write: 00070 self.resultDir=resultDir+strftime('/%Y_%m_%d/') 00071 if os.access(self.resultDir,os.F_OK) is False: 00072 print self.resultDir+' directory DNE, so am making one!' 00073 os.mkdir(self.resultDir) 00074 else: 00075 print self.resultDir+' directory exists; will add to it!' 00076 self.imDir=imDir 00077 if os.access(imDir,os.F_OK) is False: 00078 print imDir+' directory DNE, so am making one!' 00079 os.mkdir(imDir) 00080 else: 00081 print imDir+' directory exists; will add to it!' 00082 00083 t=localtime( time() ) 00084 self.fname='Regression-%s-%s-%s-%s-%s-%s.html'%(t[0],t[1],t[2],t[3],t[4],t[5]) 00085 self.html=self.resultDir+self.fname 00086 self.body1=[] 00087 self.body2=[] 00088 self.htmlPub=htmlPub(self.html,'Image tests') 00089 else: 00090 print 'stats-only mode; will not write to html file!' 00091 00092 def changeImage(self, imageName): 00093 self.imTool.open(imageName) #open('tables/squint_corr.restored') 00094 self.getArr() #instead make explicit call to getArr() 00095 self.imageName=imageName 00096 00097 def getArr(self): 00098 d=self.imTool.getchunk(axes=[],dropdeg=False) 00099 self.b=pylab.array(d) 00100 00101 00102 def simple_stats(self,sigma=10,plane=0): 00103 nchan=self.imTool.shape()[3] 00104 rmsmax=0 00105 rms1=0 00106 chan=0 00107 00108 for k in range(0,nchan/2+1): 00109 rms1=pylab.rms_flat(self.b[:,:,plane,k]) 00110 if(rms1 > rmsmax): 00111 rmsmax=rms1 00112 chan=k 00113 00114 rms1=rmsmax 00115 min1,max1=self.min_max(self.b[:,:,plane,chan]) 00116 self.show(self.b[:,:,plane,chan]) 00117 if self.write: 00118 header='Channel %d pol %d from image %s'%(chan,plane,self.imageName) 00119 body1=['The image generated with pylab:'] 00120 body2=['maximum: %f'%(max1),'minimum: %f'%(min1),'rms: %f'%(rms1)] 00121 #saveDir=self.imDir+self.fname[11:-5]+'-channel%d-pol%d.png'%(chan,plane) 00122 listnam=string.split(self.imTool.name(strippath=False), '/') 00123 imnam=listnam[len(listnam)-2]+'_'+listnam[len(listnam)-1] 00124 saveDir=self.imDir+imnam+'-channel%d-pol%d.png'%(chan,plane) 00125 pylab.savefig(saveDir) 00126 self.htmlPub.doBlk(body1, body2, saveDir,header) 00127 returnFlag= 1 00128 if(rms1 > 2*sigma): returnFlag=-1 00129 return rms1, max1, min1, returnFlag 00130 00131 def min_max(self,image): 00132 s=numarray.shape(image) 00133 min=99 00134 max=0 00135 for x in range(s[0]): 00136 for y in range(s[1]): 00137 if image[x,y]>max: max=image[x,y] 00138 if image[x,y]<min: min=image[x,y] 00139 return min,max 00140 00141 def p_min_max(self,n=5): #returns positions of n brightest pixels 00142 b=self.b 00143 s=numarray.shape(b) 00144 val=[] 00145 xp=[] 00146 yp=[] 00147 avoid=[] 00148 for z in range(n): 00149 v=0 00150 max=0 00151 x_p=0 00152 y_p=0 00153 for x in range(s[0]): 00154 for y in range(s[1]): 00155 if [x,y] not in avoid: 00156 v=b[x,y,0,0] 00157 if v>max: 00158 max=v 00159 x_p=x 00160 y_p=y 00161 val.append(max) 00162 xp.append(x_p) 00163 yp.append(y_p) 00164 avoid.append([x_p,y_p]) #ignores previous bright peaks 00165 return val,xp,yp 00166 00167 def done(self) : 00168 if self.write: 00169 self.htmlPub.doFooter() 00170 print 'webpage construction successful!' 00171 print 'images in '+os.path.abspath(self.imDir) 00172 print 'webpage at '+os.path.abspath(self.html) 00173 return '%s'%(os.path.abspath(self.html)) 00174 else: #return 0 if no writing of file is done 00175 return 'none' 00176 00177 def model(self,xx=1024,yy=1024): #make model array/image 00178 c=[] 00179 for x in range(xx): 00180 c.append([]) 00181 for y in range(yy): 00182 c[x].append([]) 00183 c[x][y].append(0.0) #note the floating point: very important! 00184 self.m=pylab.array(c) 00185 00186 def amodel(self,v,x,y): #alter selected 'pixels' of model 00187 self.m[x,y,0]=v 00188 00189 00190 def bmodel(self, XY=None, plane=0): 00191 shp=self.imTool.shape() 00192 result=[] 00193 blc=[int(0),int(0),int(plane),int(0)] 00194 trc=[int(shp[0]-1), int(shp[1]-1), int(plane), int(0)] 00195 reg=self.rgTool.box(blc=blc, trc=trc) 00196 dat=self.imTool.getchunk(blc=blc, trc=trc, dropdeg=True) 00197 self.show(dat) 00198 a={'return':{}} 00199 residual = 'framework.resid.tmp' 00200 00201 try: 00202 a = self.imTool.fitcomponents(region=reg, residual=residual) 00203 if (not a['converged']): 00204 a = self._retryFit(shp, plane, residual) 00205 except: 00206 a = self._retryFit(shp, plane, residual) 00207 resid = [] 00208 if (a['converged']): 00209 origName = self.imTool.name() 00210 self.imTool.open(residual) 00211 resid = self.imTool.getchunk() 00212 residshape = resid.shape 00213 resid = resid.reshape(residshape[0], residshape[1]) 00214 self.imTool.open(origName) 00215 else: 00216 resid = pylab.array([]) 00217 body2=[] 00218 fit_results = a['results'] 00219 tostr=lambda a: str(a[0]) if (type(a)==list) else str(a) 00220 if(fit_results.has_key('component0')): 00221 ra = tostr(self.qaTool.time(fit_results['component0']['shape']['direction']['m0'])) 00222 dec = tostr(self.qaTool.angle(fit_results['component0']['shape']['direction']['m1'])) 00223 bmaj= tostr(self.qaTool.angle(fit_results['component0']['shape']['majoraxis'], form='dig2')) 00224 bmin = tostr(self.qaTool.angle(fit_results['component0']['shape']['minoraxis'], form='dig2')) 00225 bpa = tostr(self.qaTool.angle(fit_results['component0']['shape']['positionangle'])) 00226 flux = str('%4.2f'%fit_results['component0']['flux']['value'][0])+fit_results['component0']['flux']['unit'] 00227 result.append([ra, dec, bmaj, bmin, bpa, flux]) 00228 ss='fit:\testimated center: %s %s\n'%(ra,dec)+'\tMajAxis : %s \tMinAxis: %s \tPosAng: %s'%(bmaj, bmin, bpa)+' flux= '+flux 00229 body2.append('<pre>%s</pre>'%ss) 00230 else: 00231 result.append(False) 00232 body2.append('<pre>Failed to converge in fitting</pre>') 00233 #write to web page 00234 if self.write: 00235 header='Image of plane%d of %s'%(plane,self.imageName) 00236 body1=['<pre>The image generated with pylab:</pre>'] 00237 #body2=['maximum: %f'%(max1),'minimum: %f'%(min1),'rms: %f'%(rms)] 00238 saveDir=self.imDir+self.fname[11:-5]+'-model%d.png'%(plane) 00239 pylab.savefig(saveDir) 00240 self.htmlPub.doBlk(body1, body2, saveDir,header) 00241 rms=0.0 00242 if(result[0] != False): 00243 self.show(resid) 00244 rms=pylab.rms_flat(resid) 00245 min1,max1=self.min_max(resid) 00246 print 'rms of residual image: %f'%(rms) 00247 #write to web page 00248 if self.write: 00249 header='Residual from plane%d of image %s'%(plane,self.imageName) 00250 body1=['<pre>The image generated with pylab:</pre>'] 00251 body2=['<pre>maximum: %f</pre>'%(max1),'<pre>minimum: %f</pre>'%(min1),'<pre>rms: %f</pre>'%(rms)] 00252 saveDir=self.imDir+self.fname[11:-5]+'-resid%d.png'%(plane) 00253 pylab.savefig(saveDir) 00254 self.htmlPub.doBlk(body1, body2, saveDir,header) 00255 return result, rms 00256 00257 def _retryFit(self, shp, plane, residual): 00258 ###lets limit the region and try again 00259 blc=[int(0.45*shp[0]),int(0.45*shp[1]) ,plane,0] 00260 trc=[int(0.55*shp[0]),int(0.55*shp[1]) ,plane,0] 00261 reg=self.rgTool.box(blc=blc, trc=trc) 00262 dat=self.imTool.getchunk(blc=blc, trc=trc, dropdeg=True) 00263 self.show(dat) 00264 a = {} 00265 try: 00266 a = self.imTool.fitcomponents(region=reg, residual=residual) 00267 except: 00268 a = {'converged': False, 'results':{}} 00269 return a 00270 00271 def bmodel_old(self,XY=None,plane=0): 00272 self.model() 00273 result=[] 00274 chi2=[] 00275 data=[] 00276 output=[] #to be returned: [x0,y0,FWHM_x,FWHM_y] 00277 body2=[] 00278 r,d=self.fitPeaks(XY,plane) 00279 for i in r: 00280 result.append(i[0]) 00281 chi2.append(i[1]) 00282 for i in range(len(d)): 00283 data.append([]) 00284 for j in d[i]: 00285 data[i].append(j[0]) 00286 #start computing models 00287 for i in range(len(result)): 00288 for point in data[i]: 00289 x=point[0] 00290 y=point[1] 00291 v=tw_func.gauss3d(result[i],[x,y]) 00292 self.amodel(v,x,y) 00293 print 'done building model' 00294 print 'stats for fits:' 00295 #rms=pylab.rms_flat(self.m[:,:,0]) 00296 #min1,max1=self.min_max(self.m[:,:,0]) 00297 for i in range(len(r)): 00298 self.show(self.m[:,:,0]) 00299 sigmaX=pylab.sqrt(1/(2*r[i][0][4])) 00300 sigmaY=pylab.sqrt(1/(2*r[i][0][2])) 00301 if(XY==None): 00302 XY=[[0,0]] 00303 ss='fit #%d:\testimated center: %s optimized center: [%.2f,%.2f]\n'%(i,XY[i],r[i][0][5],r[i][0][3])+'\tFWHM in x: %.3f pixels FWHM in y: %.3f\n\tchi2: %f'%(2.355*sigmaX,2.355*sigmaY,r[i][1]) 00304 print ss 00305 body2.append('<pre>%s</pre>'%ss) 00306 output.append([r[i][0][5],r[i][0][3],2.355*sigmaX,2.355*sigmaY]) 00307 #write to web page 00308 if self.write: 00309 header='Model of plane%d of image %s'%(plane,self.imageName) 00310 body1=['<pre>The image generated with pylab:</pre>'] 00311 #body2=['maximum: %f'%(max1),'minimum: %f'%(min1),'rms: %f'%(rms)] 00312 saveDir=self.imDir+self.fname[11:-5]+'-model%d.png'%(plane) 00313 pylab.savefig(saveDir) 00314 self.htmlPub.doBlk(body1, body2, saveDir,header) 00315 #return stuff 00316 return output 00317 00318 def subtract(self,plane=0): 00319 resid=[] 00320 b=self.b[:,:,plane,0] 00321 m=self.m[:,:,0] 00322 s=numarray.shape(b) 00323 for x in range(s[0]): 00324 resid.append([]) 00325 for y in range(s[1]): 00326 resid[x].append(0.0) #note the floating point: very important! 00327 self.resid=pylab.array(resid) 00328 for x in range(s[0]): 00329 for y in range(s[1]): 00330 self.resid[x,y]=b[x,y]-m[x,y] 00331 self.show(self.resid) 00332 rms=pylab.rms_flat(self.resid) 00333 min1,max1=self.min_max(self.resid) 00334 self.show(self.resid) 00335 print 'rms of residual image: %f'%(rms) 00336 #write to web page 00337 if self.write: 00338 header='Residual from plane%d of image %s'%(plane,self.imageName) 00339 body1=['<pre>The image generated with pylab:</pre>'] 00340 body2=['<pre>maximum: %f</pre>'%(max1),'<pre>minimum: %f</pre>'%(min1),'<pre>rms: %f</pre>'%(rms)] 00341 saveDir=self.imDir+self.fname[11:-5]+'-resid%d.png'%(plane) 00342 pylab.savefig(saveDir) 00343 self.htmlPub.doBlk(body1, body2, saveDir,header) 00344 #return rms 00345 return rms 00346 00347 def findPeaks(self,y,x,plane=0): #note inversion in x,y coord 00348 #use like findPeaks(x,y) 00349 #inversion should be self-consistent if take output as [x,y] 00350 r=50 #search 'radius' 00351 rms=pylab.rms_flat(self.b[:,:,plane,0]) 00352 thold=rms #limit of pixels to be considered 00353 00354 val=[] 00355 xp=[] 00356 yp=[] 00357 avoid=[] 00358 data=[] 00359 if(x==None): 00360 x=[len(self.b[:,0,plane,0])/2] 00361 y=[len(self.b[0,:,plane,0])/2] 00362 r=min(x[0], y[0])-2 00363 for i in range(len(x)): #len(x)==len(y)! 00364 v=0 00365 max=0 00366 x_p=0 00367 y_p=0 00368 for j in range(x[i]-r,x[i]+r): 00369 for k in range(y[i]-r,y[i]+r): 00370 if [j,k] not in avoid: 00371 v=self.b[j,k,plane,0] 00372 if v>max: 00373 max=v 00374 x_p=j 00375 y_p=k 00376 val.append(max) 00377 xp.append(x_p) 00378 yp.append(y_p) 00379 avoid.append([x_p,y_p]) #ignores previous bright pixels 00380 #generate data sets 00381 for i in range(len(x)): 00382 vv=val[i] 00383 xx=xp[i] 00384 yy=yp[i] 00385 while vv>thold: 00386 xx+=1 00387 vv=self.b[xx,yp[i],plane,0] 00388 vv=val[i] 00389 while vv>thold: 00390 yy+=1 00391 vv=self.b[xp[i],yy,plane,0] 00392 dx=xx-xp[i] 00393 dy=yy-yp[i] 00394 r=(dx+dy)/2 00395 if(r >12): 00396 r=12 00397 print 'PEAK pos, val and r ', xp[i], yp[i],val[i], r 00398 data.append([]) 00399 for j in range(xp[i]-r,xp[i]+r): 00400 for k in range(yp[i]-r,yp[i]+r): 00401 if (pylab.sqrt((xp[i]-j)*(xp[i]-j)+(yp[i]-k)*(yp[i]-k)) <=r) & (self.b[j,k,plane,0]>0): 00402 data[i].append(([j,k],self.b[j,k,plane,0])) 00403 return data,xp,yp 00404 00405 def fitPeaks(self,XY=None,plane=0): 00406 x=[] 00407 y=[] 00408 if(XY != None): 00409 for n in XY: 00410 x.append(n[0]) 00411 y.append(n[1]) 00412 else: 00413 x=None 00414 y=None 00415 data,x0,y0=self.findPeaks(x,y,plane) 00416 00417 tparam=[] 00418 result=[] 00419 for i in range(len(data)): #len(data) is number of fits to be made 00420 tparam.append([0.,1.,1.,x0[i],1.,y0[i]]) 00421 #tparam.append([1.,1.,1.,y0[i],1.,x0[i]]) 00422 # r=LeastSquares.leastSquaresFit(tw_func.gauss3d,tparam[i],data[i]) 00423 result.append(r) 00424 return result,data 00425 00426 def show(self,image): 00427 pylab.clf() 00428 zmin=-2.0*pylab.rms_flat(image) 00429 zmax=5.0*pylab.rms_flat(image) 00430 # pylab.imshow(image,interpolation='bilinear', origin='lower', cmap=pylab.cm.Greys, vmin=zmin, vmax=zmax) 00431 pylab.imshow(image,interpolation='bilinear', origin='lower', cmap=pylab.cm.RdYlBu_r, vmin=zmin, vmax=zmax) 00432 pylab.colorbar() 00433 def ishow(self,plane=0): #show current image 00434 pylab.clf() 00435 image=self.b[:,:,plane,0] 00436 zmin=-2.0*pylab.rms_flat(image) 00437 zmax=5.0*pylab.rms_flat(image) 00438 pylab.imshow(image,interpolation='bilinear', origin='lower', cmap=pylab.cm.Greys, vmin=zmin, vmax=zmax) 00439 pylab.colorbar() 00440 00441 def mshow(self): #show currentmodel 00442 pylab.clf() 00443 image=self.m[:,:,0] 00444 zmin=-2.0*pylab.rms_flat(image) 00445 zmax=5.0*pylab.rms_flat(image) 00446 pylab.imshow(image,interpolation='bilinear', origin='lower', cmap=pylab.cm.Greys, vmin=zmin, vmax=zmax) 00447 pylab.colorbar() 00448 00449 def rshow(self): #show current residual 00450 pylab.clf() 00451 image=self.resid 00452 zmin=-2.0*pylab.rms_flat(image) 00453 zmax=5.0*pylab.rms_flat(image) 00454 pylab.imshow(image,interpolation='bilinear', origin='lower', cmap=pylab.cm.Greys, vmin=zmin, vmax=zmax) 00455 pylab.colorbar() 00456 00457 def save(self,path='junk.png'): 00458 pylab.savefig(path) 00459 print 'image saved to %s'%(path) 00460 00461 def ishift(self,image,x,y,plane=None,cube=None): #shift the image by x,y 00462 #x,y > 0 shift to right, up 00463 #x,y < 0 shift to left, down 00464 if (plane!=None) & (cube!=None): #choose both plane and cube, if ever applicable 00465 image=image[:,:,plane,cube] 00466 elif plane!=None: #select plane of image 00467 image=image[:,:,plane,0] 00468 elif cube!=None: #select section of image cube 00469 image=image[:,:,0,cube] 00470 s=numarray.shape(image) 00471 c=[] 00472 for i in range(s[0]): 00473 c.append([]) 00474 for j in range(s[1]): 00475 c[i].append(0.0) 00476 c=pylab.array(c) 00477 for i in range(s[0]): 00478 for j in range(s[1]): 00479 try: 00480 dx=i-x 00481 dy=j-y 00482 if (dx>=0) and (dy>=0): 00483 c[i,j]=image[dx,dy] 00484 else: 00485 c[i,j]=0 00486 except IndexError: 00487 c[i,j]=0 00488 self.disp(c) 00489 return c 00490 00491 def ima(self,im1,im2): #multiply,add images 00492 s=numarray.shape(im1) #assume image shapes are same 00493 c=0 00494 for i in range(s[0]): 00495 for j in range(s[1]): 00496 c+=im1[i,j]*im2[i,j] 00497 return c 00498 00499 def mima(self,im1,im2,XY1,XY2,p1=None,p2=None,c1=None,c2=None): #c1,c2 for image cube slices only 00500 x,y,x2,y2=[],[],[],[] 00501 for n in XY1: 00502 x.append(n[0]) 00503 y.append(n[1]) 00504 for n in XY2: 00505 x2.append(n[0]) 00506 y2.append(n[1]) 00507 result=[] 00508 for n in range(len(x)): 00509 c=self.ishift(im1,x[n],y[n],p1,c1) 00510 d=self.ishift(im2,x2[n],y2[n],p2,c2) 00511 result.append(self.ima(c,d)) 00512 return result 00513 00514 def disp(self,image,plane=None,cube=None,xlabel='',ylabel=''): #cube arg applies to image cube of shape (:,:,1:,n) 00515 #cube arg specifies nth slice of cube 00516 if (plane!=None) & (cube!=None): #choose both plane and cube, if ever applicable 00517 image=image[:,:,plane,cube] 00518 elif plane!=None: #select plane of image 00519 image=image[:,:,plane,0] 00520 elif cube!=None: #select section of image cube 00521 image=image[:,:,0,cube] 00522 zmin=-2.0*pylab.rms_flat(image) 00523 zmax=5.0*pylab.rms_flat(image) 00524 00525 pylab.clf() 00526 pylab.xlabel(xlabel) 00527 pylab.ylabel(ylabel) 00528 pylab.imshow(image,interpolation='bilinear', origin='lower', cmap=pylab.cm.Greys, vmin=zmin, vmax=zmax) 00529 00530 #utilties for image cubes 00531 def drill(self,image,x0,y0): #'drills' into cube and pulls all pixels having coord x,y 00532 s=numarray.shape(image) 00533 x,y=[],[] 00534 for n in range(s[3]): 00535 x.append(n) 00536 y.append(image[x0,y0,0,n]) 00537 pylab.clf() 00538 return x,y 00539 00540 def cubePeaks(self,image,x,y): 00541 dataMax=[] 00542 dataMin=[] 00543 optimumXY=[] 00544 xp=[0,0] 00545 minX=0 00546 minY=999999 00547 maxX=0 00548 maxY=0 00549 mean=pylab.mean(y) 00550 #find max point 00551 for n in x: 00552 if y[n]>maxY: 00553 maxX=n 00554 maxY=y[n] 00555 xp[0]=xp[1]=maxX 00556 dataMax.append((xp[0],y[xp[0]])) 00557 #input data on each side of peak 00558 try: 00559 while y[xp[0]]>mean: #mean->tHold #ensures ~1 point before max 00560 xp[0]-=1 00561 dataMax.append((xp[0],y[xp[0]])) 00562 except IndexError: pass 00563 try: 00564 while y[xp[1]]>mean: #mean->tHold #ensures ~1 point beyond max 00565 xp[1]+=1 00566 dataMax.append((xp[1],y[xp[1]])) 00567 except IndexError: pass 00568 #find min point 00569 for n in x: 00570 if y[n]<minY: 00571 minX=n 00572 minY=y[n] 00573 xp[0]=xp[1]=minX 00574 dataMin.append((xp[0],y[xp[0]])) 00575 #input data on each side of peak 00576 try: #try statements guard against guassian peak going out bounds 00577 while y[xp[0]]<mean: 00578 xp[0]-=1 00579 dataMin.append((xp[0],y[xp[0]])) 00580 except IndexError: pass 00581 try: 00582 while y[xp[1]]<mean: 00583 xp[1]+=1 00584 dataMin.append((xp[1],y[xp[1]])) 00585 except IndexError: pass 00586 optimumXY.append([minX,minY]) 00587 optimumXY.append([maxX,maxY]) 00588 s='minX: %d minY: %f maxX: %d maxY: %f mean: %f'%(minX,minY,maxX,maxY,mean) 00589 if self.write: self.body2.append('<pre>%s</pre>'%s) 00590 print s 00591 return dataMin,dataMax,optimumXY 00592 def auto_fitCube2(self,verbose=0): 00593 ib=self.imTool.moments(moments=7, outfile='Some_momtest_7.im', overwrite=True) 00594 a=ib.statistics() 00595 x=int(a['maxpos'][0]) 00596 y=int(a['maxpos'][1]) 00597 ib.done() 00598 XY,fwhm=self.fitCube2(x,y) 00599 return XY,fwhm 00600 00601 def auto_fitCube(self,image,verbose=0): 00602 x,y=[],[] 00603 xlist,ylist=[],[] 00604 max_rms=0 00605 xp,yp=0,0 00606 diff=0 00607 max_diff=0 00608 s=numarray.shape(image) 00609 if s[0]>100: 00610 # xlist=range(.25*s[0],.75*s[0]) 00611 xlist=range(s[0]/4,3*s[0]/4) 00612 else: 00613 xlist=range(s[0]) 00614 if s[1]>100: 00615 # ylist=range(.25*s[1],.75*s[1],4) #later add the 4th pixel thing 00616 ylist=range(s[1]/4,3*s[1]/4,4) 00617 else: 00618 ylist=range(s[1]) 00619 for i in xlist: 00620 for j in ylist: 00621 x,y=self.drill(image,i,j) 00622 rms=pylab.rms_flat(y) 00623 diff=abs(pylab.max(y)-rms) 00624 if(diff > max_diff): 00625 max_diff,xp,yp=diff,i,j 00626 # if rms>max_rms: 00627 # max_rms,xp,yp=rms,i,j 00628 # if verbose: print i,j,rms 00629 print 'optimum fit to [%d,%d] rms=%.6f'%(xp,yp,max_rms) 00630 XY,fwhm=self.fitCube2(xp,yp) 00631 return XY,fwhm 00632 00633 def fitCube(self,image,x0,y0,fit=True): #fits all layers of cube at pixel x0,y0 00634 #fits min/max peaks w/ gaussian 00635 x,y=self.drill(image,x0,y0) 00636 result=[] 00637 if fit: 00638 maxFitted=0 #boolean flag 00639 minFitted=0 00640 s1,s2,t1,t2=[],[],[],[] 00641 s='fit stats (numbered in order of fits)' 00642 if self.write: self.body2.append('<pre>%s</pre>'%s) #assemble body2 for html 00643 print s 00644 dataMin,dataMax,optimumXY=self.cubePeaks(image,x,y) 00645 tparam0=[1.,-1.,1.,optimumXY[0][0]] 00646 tparam1=[1.,1.,1.,optimumXY[1][0]] 00647 try: #new min fit: reduce dataset size whenever overflow error occurs 00648 dataMin.sort() 00649 print 'tparam0', tparam0 00650 print dataMin 00651 dataMin.sort() 00652 # result.append(LeastSquares.leastSquaresFit(tw_func.gauss,tparam0,list(dataMin))) 00653 minFitted=1 00654 except: 00655 print '!\n overFlowError in fitting max peak\n will reduce dataset size' 00656 # print ' initial dataset size: %d'%len(dataMin) 00657 # shrink=0 00658 # while (shrink<25) & (minFitted==0): 00659 # dataMin=dataMin[:-1] 00660 # try: 00661 # result.append(leastSquaresFit(model=tw_func.gauss,parameters=tparam0,data=dataMin, stopping_limit=0.1)) 00662 # minFitted=1 00663 # shrink+=1 00664 # except: pass 00665 # print ' final dataset size: %d\n!'%len(dataMin) 00666 # except: 00667 # print 'error in finding a negative peak' 00668 try: #new max fit: reduce dataset size when things don't work right 00669 dataMax.sort() 00670 print tparam1 00671 # result.append(LeastSquares.leastSquaresFit(tw_func.gauss,tparam1,list(dataMax))) 00672 maxFitted=1 00673 except: 00674 print 'error in fitting a positive peak' 00675 # print '!\n overFlowError in fitting max peak\n will reduce dataset size' 00676 # print ' initial dataset size: %d'%len(dataMax) 00677 # shrink=0 00678 # while (shrink<25) & (maxFitted==0): 00679 # dataMax=dataMax[:-1] 00680 # try: 00681 # result.append(leastSquaresFit(tw_func.gauss,tparam1,dataMax)) 00682 # maxFitted=1 00683 # shrink+=1 00684 # except: pass 00685 # print ' final dataset size: %d\n!'%len(dataMax) 00686 #resume old 00687 for i in range(len(result)): 00688 print 'result ', result[i] 00689 sigma=pylab.sqrt(1/(2*result[i][0][2])) 00690 s='fit #%d:\toptimized coord: [%.6f,%.6f]\n\tFWHM: %f pixels\n\tchi2: %f'%(i,result[i][0][3],result[i][0][0]+result[i][0][1],2.355*sigma,result[i][1]) 00691 print s 00692 if self.write: self.body2.append('<pre>%s</pre>'%s) 00693 #draw nice guassian curves on chart 00694 if minFitted: 00695 tmax=0 00696 tmin=999999 00697 for j in dataMin: #->better gaussian chart 00698 s1.append(j[0]) 00699 for j in s1: 00700 if j>tmax: tmax=j 00701 if j<tmin: tmin=j 00702 s1=pylab.arange(tmin,tmax,0.1) 00703 for j in s1: 00704 t1.append(tw_func.gauss(result[0][0],j)) 00705 pylab.plot(s1,t1,'k.-.') #plot for min peak 00706 if maxFitted: 00707 tmax=0 00708 tmin=999999 00709 #print dataMax 00710 for j in dataMax: #->better gaussian chart 00711 s2.append(j[0]) 00712 for j in s2: 00713 if j>tmax: tmax=j 00714 if j<tmin: tmin=j 00715 #print tmin,tmax 00716 s2=pylab.arange(tmin,tmax,0.1) 00717 for j in s2: 00718 if minFitted: 00719 t2.append(tw_func.gauss(result[1][0],j)) 00720 else: 00721 t2.append(tw_func.gauss(result[0][0],j)) 00722 pylab.plot(s2,t2,'k.-.') #plot for max peak 00723 pylab.xlabel('Layer of image cube @ pixel [%d,%d]'%(x0,y0)) 00724 pylab.ylabel('Intensity') 00725 #pylab.plot(x,y,lw=0.4) 00726 pylab.plot(x,y,'b.-.',lw=0.4) 00727 #webpage write out 00728 if self.write: 00729 header='Image cube: %s'%(self.imageName) 00730 self.body1=['<pre>Plot of layers @ pixel [%d,%d]:</pre>'%(x0,y0)] 00731 saveDir=self.imDir+self.fname[11:-5]+'-cube-min_max%d.png'%self.iterate 00732 pylab.savefig(saveDir) 00733 self.htmlPub.doBlk(self.body1, self.body2, saveDir,header) 00734 self.iterate+=1 00735 #now return stuff: center coord, FWHM values 00736 XY=[] 00737 fwhm=[] 00738 if(len(result)==0): 00739 XY.append([-1., -1.]) 00740 XY.append([-1., -1.]) 00741 fwhm.append(-1) 00742 fwhm.append(-1) 00743 00744 for i in range(len(result)): 00745 sigma=pylab.sqrt(1/(2*result[i][0][2])) 00746 fwhm.append(2.355*sigma) 00747 XY.append([result[i][0][3],result[i][0][0]]) 00748 00749 return XY,fwhm 00750 00751 def fitCube2(self,x0,y0): 00752 result=[] 00753 box = str(x0) + ", " + str(y0) + ", " + str(x0) + ", " + str(y0) 00754 reg=self.rgTool.box(blc=[x0,y0], trc=[x0,y0]) 00755 model = "mymodel.im" 00756 resid = "myresid.im" 00757 myia = casac.image() 00758 for x in ([model, resid]): 00759 if (os.path.exists(x)): 00760 myia.open(x) 00761 myia.remove() 00762 myia.close() 00763 retval=self.imTool.fitprofile(box=box, ngauss=1, model=model, residual=resid) 00764 myia.open(model) 00765 theFit = myia.getchunk().flatten() 00766 myia.close() 00767 myia.open(resid) 00768 theResid = myia.getchunk().flatten() 00769 myia.done() 00770 for x in ([model, resid]): 00771 if (os.path.exists(x)): 00772 myia.open(x) 00773 myia.remove() 00774 myia.close() 00775 result.append([retval['gs']['amp'].flatten()[0], retval['gs']['center'].flatten()[0], retval['gs']['fwhm'].flatten()[0]]) 00776 result.append([retval['gs']['ampErr'].flatten()[0], retval['gs']['centerErr'].flatten()[0], retval['gs']['fwhmErr'].flatten()[0]]) 00777 s='fit at [%d,%d]\n\tFWHM: %f \n\peak: %f \t with errors: %f, %f '%(x0,y0, result[0][2], result[0][0], result[1][2], result[1][0]) 00778 print s 00779 if self.write: self.body2.append('<pre>%s</pre>'%s) 00780 data=self.imTool.getchunk(blc=[int(x0),int(y0)], trc=[int(x0),int(y0)], dropdeg=True) 00781 pylab.clf() 00782 pylab.plot(data,'r-') 00783 pylab.plot(theFit, 'bo') 00784 pylab.plot(theResid,'g.-.') 00785 pylab.xlabel('Layer of image cube @ pixel [%d,%d]'%(x0,y0)) 00786 pylab.ylabel('Intensity') 00787 pylab.title('Fit on Image '+self.imageName) 00788 if self.write: 00789 header='Image cube Fit: %s'%self.imageName 00790 self.body1=['<pre>Plot of layers @ pixel [%d,%d]:</pre>'%(x0,y0)] 00791 saveDir=self.imDir+self.fname[11:-5]+'-cube-min_max%d.png'%self.iterate 00792 pylab.savefig(saveDir) 00793 self.htmlPub.doBlk(self.body1, self.body2, saveDir,header) 00794 self.iterate+=1 00795 XY=[] 00796 fwhm=[] 00797 for i in range(len(result)): 00798 fwhm.append(result[i][2]) 00799 XY.append([result[i][0],result[i][1]]) 00800 00801 return XY, fwhm
1.8.0
