|
casa
$Rev:20696$
|
|
casa
$Rev:20696$
|
00001 from taskinit import * 00002 import numpy 00003 import sys 00004 import re 00005 00006 # Writes out regions above threshold to regionfile+'.box' 00007 00008 def boxit(imagename, regionfile, threshold, maskname, chanrange, polrange, minsize, diag, boxstretch, overwrite): 00009 00010 casalog.origin('boxit') 00011 00012 if not(regionfile): 00013 regionfile = imagename + '.box' 00014 if not regionfile.endswith('.box'): 00015 regionfile = regionfile + '.box' 00016 00017 if not(overwrite): 00018 if(os.path.exists(regionfile)): 00019 casalog.post('file "' + regionfile + '" already exists.', 'WARN') 00020 return 00021 if(maskname and os.path.exists(maskname)): 00022 casalog.post('output mask "' + maskname + '" already exists.', 'WARN') 00023 return 00024 00025 # If no units, assume mJy for consistency with auto/clean tasks. 00026 # But convert to Jy, because that's what units the images are in. 00027 threshold = qa.getvalue(qa.convert(qa.quantity(threshold,'mJy'),'Jy'))[0] 00028 casalog.post("Setting threshold to " + str(threshold) + "Jy", "INFO") 00029 00030 newIsland = numpy.zeros(1, dtype=[('box','4i4'),('npix','i4')]) 00031 # Find all pixels above the threshold 00032 ia.open(imagename) 00033 # CAS-2059 escape characters in image name that will confuse the lattice expression processor 00034 escaped_imagename = imagename 00035 for escapeme in ['-', '+', '*', '/' ]: 00036 escaped_imagename = re.sub("[" + escapeme + "]", "\\" + escapeme, escaped_imagename) 00037 mask = escaped_imagename+'>'+str(threshold) 00038 fullmask = ia.getregion(mask=mask, getmask=True) 00039 if not(fullmask.max()): 00040 casalog.post('Maximum flux in image is below threshold.', 'WARN') 00041 return 00042 writemask = bool(maskname) 00043 csys = ia.coordsys() 00044 00045 shape = fullmask.shape 00046 00047 nx = int(shape[0]) 00048 ny = int(shape[1]) 00049 n2 = n3 = 1 00050 if len(shape)==3: 00051 n2 = int(shape[2]) 00052 if len(shape)==4: 00053 n2 = int(shape[2]) 00054 n3 = int(shape[3]) 00055 #print 'nx:', nx, 'ny:', ny, 'n2:', n2, 'n3:', n3 00056 00057 #casa generated images always 4d and in order of [ra, dec, stokes, freq] 00058 #other images can be in the order of [ra, dec, freq, stokes] 00059 nms = csys.names() 00060 #axes=[] 00061 #for i in xrange(len(nms)): 00062 # if nms[i]=='Right Ascension': 00063 # axes[0]=i 00064 # if nms[i]=='Declination': 00065 # axes[1]=i 00066 # if nms[i]=='Stokes' 00067 # axes[2]=i 00068 # if nms[i]=='Frequency': 00069 # axes[3]=i 00070 00071 chmax=n3 00072 pomax=n2 00073 chmin=0 00074 pomin=0 00075 if len(nms)==4 and nms[3]=='Stokes': 00076 chmax=n2 00077 pomax=n3 00078 00079 if chanrange: 00080 try: 00081 if str.count(chanrange, '~') == 1: 00082 ch1=int(str.split(chanrange, '~')[0]) 00083 ch2=int(str.split(chanrange, '~')[1]) 00084 else: 00085 ch1=int(chanrange) 00086 ch2=ch1 00087 if ch1 > ch2 or ch1 < 0 or ch2 > chmax: 00088 casalog('invalid channel range', "SEVERE") 00089 return 00090 if ch1>chmin: 00091 chmin = ch1 00092 if ch2<chmax: 00093 chmax = ch2 00094 except: 00095 casalog('bad format for chanrange', "SEVERE") 00096 return 00097 #print 'chmin', chmin, 'chmax', chmax 00098 00099 if polrange: 00100 try: 00101 if str.count(polrange, '~') == 1: 00102 po1=int(str.split(polrange, '~')[0]) 00103 po2=int(str.split(polrange, '~')[1]) 00104 else: 00105 po1=int(polrange) 00106 po2=po1 00107 if po1 > po2 or po1 < 0 or po2 > pomax: 00108 casalog( 'invalid stokes range', "SEVERE") 00109 return 00110 if po1>pomin: 00111 pomin = po1 00112 if po2<pomax: 00113 pomax = po2 00114 except: 00115 casalog( 'bad format for polrange', "SEVERE") 00116 return 00117 #print 'pomin', pomin, 'pomax', pomax 00118 00119 n1=pomin 00120 n2=pomax 00121 n3=chmax 00122 n4=chmin 00123 if len(nms)==4 and nms[3]=='Stokes': 00124 n1=chmin 00125 n2=chmax 00126 n3=pomax 00127 n4=pomin 00128 #print 'n1:', n1, 'n2:', n2, 'n4:', n4, 'n3:', n3 00129 00130 f = open(regionfile, 'w') 00131 totregions = 0 00132 outputmask = [] 00133 if writemask: 00134 outputmask = ia.getchunk() 00135 outputmask.fill(False) 00136 00137 for i3 in xrange(n4, n3): 00138 for i2 in xrange(n1, n2): 00139 regions = {} 00140 boxRecord = {} 00141 if len(shape)==2: 00142 mask = fullmask 00143 if len(shape)==4: 00144 mask = fullmask[:,:,i2,i3].reshape(nx, ny) 00145 islands = [] 00146 pos = numpy.unravel_index(mask.argmax(), mask.shape) 00147 #print pos 00148 while(mask[pos]): 00149 # found pixel in new island 00150 island = {} 00151 island['box'] = [pos[0], pos[1], pos[0], pos[1]] 00152 island['npix'] = 1 00153 find_nearby_island_pixels(island, mask, pos, shape[0], shape[1], diag) 00154 islands.append(island) 00155 pos = numpy.unravel_index(mask.argmax(), mask.shape) 00156 for record in islands: 00157 box = record['box'] 00158 # check size of box 00159 if record['npix'] < minsize: 00160 continue 00161 totregions += 1 00162 # need pixel corners, not pixel centers 00163 box[0] -= boxstretch 00164 box[1] -= boxstretch 00165 box[2] += boxstretch 00166 box[3] += boxstretch 00167 # in case we used boxstretch < 0 and a one-pixel sized box: 00168 if box[0] > box[2]: 00169 box[0] += boxstretch 00170 box[2] -= boxstretch 00171 if box[1] > box[3]: 00172 box[1] += boxstretch 00173 box[3] -= boxstretch 00174 if writemask: 00175 # avoid pixels in boxes that have been stretched beyond the image limits 00176 for ii in range(max(0, box[0]), min(box[2], outputmask.shape[0] - 1)): 00177 for jj in range(max(0, box[1]), min(box[3], outputmask.shape[1] - 1)): 00178 outputmask[ii][jj][i2][i3] = True 00179 blccoord = [box[0]-0.5, box[1]-0.5, i2, i3] 00180 trccoord = [box[2]+0.5, box[3]+0.5, i2, i3] 00181 # note that the toworld() calls are likely very expensive for many boxes. But then 00182 # again, the box-finding algorithm itself seems pretty inefficient, but resource 00183 # constraints only permit a band aid fix at this time. 00184 blc = ia.toworld(blccoord, 'm')['measure'] 00185 trc = ia.toworld(trccoord, 'm')['measure'] 00186 # RA/Dec reference frame 00187 outstring = "worldbox " + blc['direction']['refer'] 00188 # RA blc/trc 00189 outstring += " [" + quantity_to_string(blc["direction"]["m0"], "rad") + ", " 00190 outstring += quantity_to_string(trc["direction"]["m0"], "rad") + "]" 00191 # Dec blc/trc 00192 outstring += " [" + quantity_to_string(blc["direction"]["m1"], "rad") + ", " 00193 outstring += quantity_to_string(trc["direction"]["m1"], "rad") + "]" 00194 # frequency blc/trc 00195 freq = blc["spectral"]['frequency']['refer'] + " " 00196 freq += quantity_to_string(blc["spectral"]["frequency"]["m0"], "Hz", False) 00197 outstring += " ['" + freq + "', '" + freq + "']" 00198 # Stokes blc/trc 00199 outstring += " ['" + blc["stokes"] + "', '" + trc["stokes"] + "']" 00200 # add the mask flag 00201 outstring += " " + str(1) 00202 f.write(outstring + "\n") 00203 casalog.post("Wrote " + str(totregions) + " regions to file " + regionfile, 'INFO') 00204 if writemask: 00205 myia = iatool() 00206 myia.fromimage(infile=imagename, outfile=maskname, overwrite=True) 00207 myia.done() 00208 myia.open(maskname) 00209 myia.putchunk(outputmask) 00210 myia.done() 00211 f.close() 00212 return True 00213 00214 def quantity_to_string(quantity, unit=None, quotes=True): 00215 if unit != None: 00216 quantity = qa.convert(quantity, unit) 00217 string = str(quantity['value']) + quantity['unit'] 00218 if quotes: 00219 string = "'" + string + "'" 00220 return string 00221 00222 def find_nearby_island_pixels(island, mask, pos, xmax, ymax, diag): 00223 # blank this position so we don't deal with it again 00224 mask[pos] = False 00225 xref = pos[0] 00226 yref = pos[1] 00227 for x in range(max(xref-1, 0), min(xref+2, xmax)): 00228 for y in range(max(yref-1, 0), min(yref+2, ymax)): 00229 if x == xref and y == yref: 00230 # same position as reference 00231 continue 00232 if ( (not diag) and (x-xref != 0) and (y-yref != 0)): 00233 # diagonal pixels only used if diag is true 00234 continue 00235 if mask[x][y]: 00236 # found another pixel in this island 00237 island['box'][0] = min(island['box'][0],x) 00238 island['box'][1] = min(island['box'][1],y) 00239 island['box'][2] = max(island['box'][2],x) 00240 island['box'][3] = max(island['box'][3],y) 00241 island['npix'] += 1 00242 # look for island pixels next to this one 00243 find_nearby_island_pixels(island, mask, (x,y), xmax, ymax, diag) 00244 return 00245 00246
1.8.0
