|NRAO Home > CASA > CASA Toolkit Reference Manual||
Create and manipulate regions of interest
Overview of Regionmanager functionality
- Create simple world-coordinate based regions with functions
- frombcs (range along one axis)
- wrange (range along one axis)
Also related is function setcoordinates
- Convert pixel regions to world regions with
- There are some general utility functions, generally only of interest if you are
writing scripts. These are
- There are some functions relating to interactive creation of regions with the
Viewer. These are generally only of interest if you are writing scripts, and
- There are some functions relating to communications. These are generally only
of interest if you are writing scripts and are
Regions and the Regionmanager
When working with an image, one is generally interested in some part of that image on which astrophysical analysis is performed. This region-of-interest (or more generically and simply, the ‘region’) might be the whole image, or some subset of it.
Regions come in a few varieties. There are simple regular shapes (box, ellipsoid), simple irregular shapes (polygon), as well as compound regions made by combining other regions (simple or compound). For example unions or intersections of regions. In addition, the simple regions can be pixel-coordinate based or world-coordinate based. However, a compound regions must always comprise world-coordinate based regions.
It is the task of the Regionmanager (tool) to manage your regions; it creates, lists and manipulates them. Apart from a Regionmanager, the only other way to create a region-of-interest is with the viewier tool(type viewer within the casapy environment). This allows you to interactively make, with the cursor and an image display, a box or polygonal region. The region so made may be collected by the Regionmanager (in fact the complete process can be initiated from the Regionmanager).
The Regionmanager has a command line interface, but there are plans to have it interact directly with the CASA viewertool. in the future. Currently the only way to interact is to save the regions created with the viewer to a file or as a region in the image , and use the Regionmanager function fromfiletorecord or fromtabletorecord, repsectively, to bring the regions in the CLI.
It is probably fair to say that for the simplest regions, such as a pixel box, there is little to choose between making a region with the GUI or the command line interface. However, for the world regions, the GUI is significantly better; it hides the details of handling the coordinate information.
Simple Pixel and World regions
Pixel regions are specified purely in terms of pixel coordinates. Pixel coordinates are considered to run from 1 at the bottom left corner (blc) of an image to the image shape at the top-right corner (trc).
For example, a pixel region might be box that runs from a bottom blc of [20,20] to a trc of [64,90] in a 2D image. Such a pixel region is NOT very portable though. If you are interested in a region about an astronomical source, pixel coordinates are only useful for a particular image. In some other image, the source may well have different pixel coordinates. However, pixel regions are easy to make and use, and they are valuable.
So far, we have been talking about absolute pixel coordinates. However, pixel regions can also be specified in relative coordinates. This is controlled through the argument absrel which may take on one of the values ’abs’ (absolute coordinates) ’relref’ (relative to the reference pixel of the image) or ’relcen’ (relative to the center of the image). You can use the summary function of the image module to see the reference pixel. The sense of the offset is rel = abs - ref.
You may also define pixel regions in terms of fractional coordinates that are in the range [0,1] (blc to trc) - look for argument frac which can be T or F.
From the command line, let’s make a simple pixel box region with
- ia.open(’myimage’) - csys = ia.coordsys() - r1 = rg.box( blc=[10,20], trc=[30,40] ) - ia.statistics(region=r1); Number points = 441 Flux density = -1.728895e-03 Jy Sum = -9.522969e-02 Mean = -2.159403e-04 Variance = 1.518106e-05 Sigma = 3.896288e-03 Rms = 3.897854e-03
Minimum value is -1.208747e-02 at [13, 32] Maximum value is 1.287862e-02 at [21, 37]
You have now created a region tool called r1 which describes a 2D box from [10,20] to [30,40] in absolute pixel coordinates. The region is made ‘const’ as well so you can’t overwrite it by mistake. We then apply it to an image and evaluate some statistics in the specified region.
World regions are more portable, because the region is specified in terms of the world coordinates (e.g. RA and DEC, or Frequency etc). You can create a region with one image, and then apply it to another. For example, you might make an RA/DEC world box (i.e. blc and trc specified as RA and DEC) from an optical image. You can then apply it to your radio image of the same source. The software will look for the axes of the region (in this case RA and DEC), and then attempt to match them with the image to which the region is being applied. For each matching axis, the region will be applied.
When you make a world region, you must supply a coordinate system from the appropriate image to the Regionmanager. The coordinate system comes from the coordsys function of the Image tool. This associates a coordinate with an image pixel axis.
Now let’s do an example with a world box. It is better to make a world box with the GUI interface, but the command line is still manageable. First, let’s summarise the header of our image (the one from the previous example).
- ia.open( ’myimage’ ) - ia.summary() Image name : myimage Image mask : Absent Image units : JY/BEAM
Direction system : J2000
Name Proj Shape Tile Coord value at pixel Coord incr Units ————————————————————————— Right Ascension SIN 155 155 17:42:29.303 90.00 -1.000000e+00 arcsec Declination SIN 178 178 -220.127.116.110 90.00 1.000000e+00 arcsec T
And now we will make a world box with the blc at the reference pixel, and the trc offset somewhat. We use the quanta module for creating associations of values and units; in particular, the quantity function.
- blc = ”17:42:29.303 -18.104.22.1680”; - trc = ”17:42:28.303 -22.214.171.1240”; - csys = ia.coordsys() - r2 = rg.wbox(blc=blc, trc=trc, pixelaxes=[0,1], csys=csys.torecord() ) - ia.boundingbox(r2); [blc=[90 90], trc=[103 98], inc=[1 1], bbShape=[14 9], regionShape=[14 9], imageShape=[155 178] ]
As well as passing in the blc and trc quantities, we have also told the Regionmanager which axes these values pertain to (they are numbered in the order in which they list in the summary function), and we also supplied it with the coordinate system of the image. Note that with the GUI interface, both the axes and coordinates handling is hidden. If we wished, we could have left out the pixelaxes argument, and it would have defaulted to [0,1]. We then find the bounding box of the region when it’s applied to the image. You can see from the summary listing that it is correct.
We generally think about our images as being in some order; RA/DEC/Frequency (for the X/Y/Z axes) or whatever. In the application of world regions, this order is unimportant. If you created a world region from an RA/DEC/Frequency image, and applied it to a Frequency/DEC/RA image, that image order change would be accounted for.
Note that for pixel regions however, because the region is not tagged with a coordinate system, the order of the image is relevant. So if you made a pixel box from an RA/DEC image and applied it to a DEC/RA image, the RA range would be applied to the DEC axis and vice versa. The pixel regions are present for simple usage only.
Simple pixel regions as world regions
It is possible to create a world region that is specified in pixel coordinates, and maintains the coordinate axis information. This can be useful because compound regions must be world-coordinate based (see below). We can do this because we have defined in the Regionmanager, “new” world units, called “pix” (which are then known to the quanta module) and these are recognized in the world regions. Similarly, we have also defined a unit called “frac”. This allows you to specify world coordinates in units of the fraction of the image shape. For example, the blc of a 2D image has a “frac” coordinate of [0,0] and a trc of [1,1]. But be careful, such regions are still not portable to another image, because they are still effectively pixel-coordinate based (although masquerading as world regions).
- ia.maketestimage() - csys = ia.coordsys() - - r1 = rg.wbox(blc=”10pix,10pix”, trc=”20pix,20pix”, pixelaxes=[0,1], csys=csys.torecord()) - localstats = ia.statistics(region=r1, list=False, verbose=False) - print ’Region 1: number points = ’, localstats[’npts’] Region 1: number points = 121 - - r2 = rg.wbox(blc=”30pix 30pix”, trc=”40pix 40pix”, pixelaxes=[0,1], csys=csys.torecord()) - localstats=ia.statistics(region=r2, list=False, verbose=False) - print ’Region 2: number points = ’, localstats[’npts’] Region 2: number points = 121
- regions= - regions[’r1’]=[r1] - regions[’r2’]=[r2] - r3 = rg.makeunion(regions) - localstats=ia.statistics(region=r3, list=False, verbose=False) - print ’Region 3: number points = ’, localstats[’npts’] Region 3: number points = 242
In this example we make a union from two pixel boxes masquerading as world regions.
It is often desirable to make a region which combines others. These are called compound regions. For example, give me the union of 3 regions and evaluate the statistics in that union. Or intersect this region with that one and extract the image data from that region. Compound regions are fully recursive; you can nest them as deeply as you like.
You must be aware that compound regions can only be made from world-coordinated based regions. You can convert a pixel region to a world region with the pixeltoworldregion function (at some point this function will move from the Regionmanager tool to the Image tool.
- ia.open(’myimage’) - csys ia.coordsys(); - - blc = ”17:42:29.303, -126.96.36.1990”; - trc = ”17:42:28.303, -188.8.131.520”; - r1 = rg.wbox(blc=blc,trc=trc,pixelaxes=[0,1],csys=csys.torecord()); - ia.boundingbox(r1); [blc=[90 90] , trc=[103 98] , inc=[1 1], bbShape=[14 9] , regionShape=[14 9] , imageShape=[155 178] ] - - blc = ”17:42:27.303, -184.108.40.2060” - trc = ”17:42:26.303, -220.127.116.110” - r2 = rg.wbox(blc=blc,trc=trc,pixelaxes=[0,1],csys=csys.torecord()); - ia.boundingbox(r2); [blc=[116 90] , trc=[129 98] , inc=[1 1], bbShape=[14 9] , regionShape=[14 9] , imageShape=[155 178] ] - - regions= - regions[’r1’]=[r1] - regions[’r2’]=[r2] - r3 = rg.makeunion(regions) - bb = ia.boundingbox(r3) - bb [blc=[90 90] , trc=[129 98] , regionShape=[40 9] , imageShape=[155 178] ] - - s1 = ia.statistics(region=r1, list=False) - s2 = ia.statistics(region=r2, list=F) - s3 = ia.statistics(region=r3, list=F) - - np = (bb[’trc’]-bb[’blc’]+1)*(bb[’trc’]-bb[’blc’]+1) - print ’No. pts in r1, r2, r3, bb(r3) = ’, s1[’npts’], s2[’npts’], s3[’npts’], np No. pts in r1, r2, r3, bb(r3) = 126, 126, 252, 360 - - - r4 = rg.intersection(r1, r2) - ia.boundingbox(r4) - r4
We make two world boxes and find the bounding box of each. Then we make the union of the two and find the bounding box which we store in a record called bb. We then find some statistics over the union (supressing the logger output and storing the results in records). You can see that the number of points found in the region reflects the union, not the bounding box of the union.
Finally, we make an intersection of our two regions. Because regions r1 and r2 don’t actually intersect, so the region returned is empty.
One general philosophy behind the regions software is that if you apply a region to an image, but don’t explicitly specify all axes of the image in that region, then for the unspecified axes, the full image is taken. For example, imagine that from an RA/DEC optical image you made an RA/DEC region. You then applied it to an RA/DEC/Frequency radio spectral-line cube. For the frequency axis, which was not specified in the region, all frequency pixels would be selected. Another philosophy is that if, on application, a region extends beyond an image, the overhang is silently discarded.
- ia.open(’myimage’) - ia.boundingbox(rg.box()) [blc=[1 1] , trc=[155 178] , inc=[1 1] , bbShape=[155 178] , regionShape=[155 178] , imageShape=[155 178] ] - - r1 = rg.box(trc=) - ia.boundingbox(r1) [blc=[1 1] , trc=[20 178] , inc=[1 1] , bbShape=[20 178] , regionShape=[20 178] , imageShape=[155 178] ] - - r2 = rg.box(trc=[9000,20]) - ia.boundingbox(r2) [blc=[1 1] , trc=[155 20] , inc=[1 1] , bbShape=[155 20] , regionShape=[155 20] , imageShape=[155 178] ]
In the first example, none of the axes are specified in the region, so it defaults to the full image on all axes. In the second example only the first axis of the trc is specified. In the third example, the trc overhang is silently discarded.
Some comment about the combination of masks and regions-of-interest is useful here. Consider a simple polygonal region. This region-of-interest is defined by a bounding box, the polygonal vertices, and a mask called a region mask. The region mask specifies whether a pixel within the bounding box is inside or outside the polygon. For a simple box region-of-interest, there is obviously no need for a region mask.
Now imagine that you wish to recover the mask of an image from a polygonal region-of-interest. Now, necessarily, the mask is returned to you in regular Boolean array. Thus, the array shape reflects the bounding-box of the polygonal region. If the actual pixel mask that you apply is all good, then the retrieved mask would be good inside of the polygonal region and bad outside of it. If the actual pixel mask had some bad values in it as well, the retrieved mask would be bad outside of the polygonal region. Inside the polygonal region it would be bad if the pixel mask was bad.
More simply put, the mask that you recover is just a logical “and” of the pixel mask and the region mask; if the pixel mask is T and the region mask is T then the retrieved mask is T (good), else it is F (bad).
Many of the functions of a Regionmanager tool take vectors (numeric or strings) as their arguments. You can assume that
- For numeric vectors you can enter as an actual numeric vector (e.g. [’1.2’,’2.5’,’3’]), a vector of strings (e.g. ”1.2 2.5 3”) or even a string with white space and/or comma delimiters (e.g. ’1.2 2.5 3’).
- For string vectors you can enter as an actual vector of strings (e.g. ”1.2 2.5 3”) or a string with white space or comma delimiters (e.g. ’1.2 2.5 3’).
When specifying blc and trc vectors for some regions (usually boxes) it can often be that you wish to default one axis but not others.
For example, you may like to specify a value for axis 2 but leave axis 1 at its defaults. This is trivial with the GUI interface (you just don’t fill in the ones you are not interested in), but is a little harder with the command line interface.
There is a function called def which provides a magic value which can be used with pixel boxes for this purpose. For example,
THIS EXAMPLE IS NOT VALID YET
- ia.fromshape(’x’, [20,30]) - r = rg.box(trc=[rg.dflt(),10]) - ia.boundingbox(region=r) [blc=[1 1] , trc=[20 10] , regionShape=[20 10] , imageShape=[20 30] ]
You can see that trc for axis 1 has defaulted to the shape of the image.
With world boxes, we have defined a magic unit called ‘default’ instead. Thus
- ia.fromshape(’x’, [20,30]) - csys = ia.coordsys() - rg.setcoordinates(csys.torecord()) - r = rg.wbox(trc=”0dflt 10pix”) - ia.boundingbox(region=r) [blc=[1 1] , trc=[20 10] , regionShape=[20 10] , imageShape=[20 30] ]
Although some error checking is done when the region is created (in CASA), much of it does not occur until the region is applied to an image. In particular, when creating a compound region, it is not checked that the compound region contains only world-coordinate based regions. It is only when you apply the region to an image that these checks are made. Any errors that occur will cause exceptions to be generated. Thus, when writing scripts, you should always use try/except blocks as in this example.
- def getStats(imageobject, blc, trc, csys): try : r1 = rg.wbox(blc=blc, trc=trc, csys=csys.torecord()) imageobject.statistics(r1) except Exception, instance: print ’ERROR: ”, instance return False
This fairly silly function detects if the region creation failed (function given invalid arguments) and displays the exception if it did. If all is well statistics are then evaluated from the region, and True is returned.
How a region is stored in CASA
In CASA, a region is itself a tool. This means it can be transmitted about with Glish, saved to Tables and restored from Tables. A region is actually made with a generic container tool called an Itemcontainer. The regions that you make have no intelligence regarding their greater purpose in life (although of course an Itemcontainer does have some functions so they can be manipulated). All of the knowledge about regions lies with the Regionmanager. It should not be necessary for you to know anything about the insides of regions. You just need to know how to make them and how to use them.
THIS EXAMPLE IS NOT VALID ANYMORE
- r1 = rg.box([10,20], [30,40]) - r1 [gui=¡function¿, done=¡function¿, get=¡function¿, has_item=¡function¿, makeconst=¡function¿, length=¡function¿, names=¡function¿, set=¡function¿, add=¡function¿, fromrecord=¡function¿, torecord=¡function¿, type=¡function¿] - - r1.torecord() [name=LCSlicer, isRegion=3, blc=[10 20] , trc=[30 40] , inc=, fracblc=[F F] , fractrc=[F F] , fracinc=, arblc=[1 1] , artrc=[1 1] , oneRel=T, comment=]
First we make a simple pixel-box region. We then try and print it out. You can see that this isn’t very enlightening. What you see are the Itemcontainer’s functions. To really see inside the region, you can use the torecord function of Itemcontainer. This converts the region to a record which you can view if you really must (the record is of no functional use when dealing with regions).
It’s perhaps also useful to know that once you have made a region you can’t break it. You can delete it, maybe overwrite it, but you can’t mistakenly mess up its contents (even with the Itemcontainer functions with which it was created). It’s user proof - there’s a challenge for you !
|regionmanager||Construct a regionmanager|
|absreltype||Convert region type value to a string|
|box||Create a pixel box region|
|frombcs||Create a world coordinate region based on box-chan-stokes input|
|complement||Create the complement of a world region|
|concatenation||Concatenate world regions along a new axis|
|deletefromtable||Delete regions from a Table|
|difference||Create the difference of two world regions|
|done||Destroy this regionmanager|
|selectedchannels||Get an array of zero-based selected channel numbers from an input string specificaiton.|
|fromtextfile||Create a region dictionary from a region text file.|
|fromtext||Create a region dictionary from a region text string.|
|fromfiletorecord||Create a region record(s) from a file(s).|
|tofile||Create a region record file that can be read by from filetorecord.|
|fromrecordtotable||Save regions stored in a record into a Table|
|fromtabletorecord||Restore regions from a Table to a record|
|intersection||Create the intersection of some world regions|
|ispixelregion||Is this region a pixel region ?|
|isworldregion||Is this region a world region ?|
|namesintable||Find the names of the regions stored in a Table|
|setcoordinates||Set new default Coordinate System|
|makeunion||Create a union of world regions|
|wbox||Create a world box region|
|wpolygon||Create a world polygon region with quantities|
regionmanager.absreltype - Function
regionmanager.box - Function
regionmanager.frombcs - Function
regionmanager.complement - Function
regionmanager.concatenation - Function
regionmanager.deletefromtable - Function
regionmanager.difference - Function
regionmanager.done - Function
regionmanager.selectedchannels - Function
regionmanager.fromtextfile - Function
regionmanager.fromtext - Function
regionmanager.fromfiletorecord - Function
regionmanager.tofile - Function
regionmanager.fromrecordtotable - Function
regionmanager.fromtabletorecord - Function
regionmanager.intersection - Function
regionmanager.ispixelregion - Function
regionmanager.isworldregion - Function
regionmanager.namesintable - Function
regionmanager.setcoordinates - Function
regionmanager.makeunion - Function
regionmanager.wbox - Function
regionmanager.wpolygon - Function
More information about CASA may be found at the CASA web page
Copyright © 2016 Associated Universities Inc., Washington, D.C.