(1) convert the ds9 region to a uniform format. Both the box and polygon regions can be treated as polygons.
 |
| an example of an image mask. The green outline is the polygon drawn into SAOimage DS9. [FROM the website of galfit] |
(2) solve the PIP (point in polygon) problems (http://en.wikipedia.org/wiki/Point_in_polygon) and mask all the points in the polygon. I attached one useful python script to get rid of the PIP problem.
If you want to make two loops in python to find all the inside points, you may find the speed is too low! Please avoid to use 'For' in python. The fast way is to use the exist program in numpy
matplotlib.nxutils.points_inside_pol as following:import numpy as np
from matplotlib.nxutils import points_inside_poly
nx, ny = 10, 10
poly_verts = [(1,1), (5,1), (5,9),(3,2),(1,1)]
# Create vertex coordinates for each grid cell...
# (<0,0> is at the top left of the grid in this system)
x, y = np.meshgrid(np.arange(nx), np.arange(ny))
x, y = x.flatten(), y.flatten()
points = np.vstack((x,y)).T
grid = points_inside_poly(points, poly_verts)
grid = grid.reshape((ny,nx))
print grid(3) save your masked array !
=======================================================================
Python script to decide whether the point is inside the polygon.
def point_in_poly(x,y,poly):
'''
Decide whether the point is inside the polygon.
'''
n = len(poly)
inside = False
p1x,p1y = poly[0]
for i in range(n+1):
p2x,p2y = poly[i % n]
if y > min(p1y,p2y):
if y <= max(p1y,p2y):
if x <= max(p1x,p2x):
if p1y != p2y:
xints = (y-p1y)*(p2x-p1x)/(p2y-p1y)+p1x
if p1x == p2x or x <= xints:
inside = not inside
p1x,p1y = p2x,p2y
return inside
=======================================================================
(1) A copy from here ( http://www.cfht.hawaii.edu/~chyan/wircam/badpix/badpix.html)
The program "wircampolyregion" and other IDL programs are attached in the end of the website.
(2) Another way and program is mentioned in GALFIT (http://users.obs.carnegiescience.edu/peng/work/galfit/MASKING.html).
=======================================================================
5. IDL Source Code
; NAME:
; MASKREGION
; PURPOSE:
; Produce a new bad pixel mask based on an exist one and ds9 region file
;
; CALLING SEQUENCE:
; MASKREGION, BADPIX = badpix, REGION = ds9.reg, NEWBADPIX= newbadpix
;
; INPUTS:
; BADPIX - string scalar, the name of the original bad pixel mask.
; REGION - DS9 region file, only polygon and box allowed.
; NEWBADPIX - string scalar, the new bad pixel mask.
;
;
; REVISION HISTORY:
; initial version Chi-hung Yan March 2007
;-
PRO RDWIRCAMCUBE, file, im
im=fltarr(2048,2048,4)
for i=0,3 do begin
im[*,*,i]=mrdfits(file,i+1)
endfor
END
;
PRO MASKPOLYREGION, im, ext, x, y
ind=polyfillv(x,y,2048,2048)
tim=im[*,*,ext-1]
tim[ind]=!values.f_nan
im[*,*,ext-1]=tim[*,*]
END
PRO WRWIRCAMCUBE, im, fits, new
main_hd=headfits(fits,exten=0)
mwrfits, 3,new , main_hd, /create
for i=1,4 do begin
hdr = headfits(fits, exten=i)
mwrfits,im[*,*,i-1],new,hdr
endfor
END
PRO MASKREGION, badpix=badpix, region=region, newbadpix=newbadpix
file=region
fits=badpix
; Read image
rdwircamcube,fits,im
spawn,"wircampolyregion "+file+" /tmp/region.txt"
readcol,"/tmp/region.txt",id,ext,x,y
n=0
xx=fltarr(1)
yy=fltarr(1)
for i=0,n_elements(x)-1 do begin
if i eq 0 then begin
temp=id[i]
endif else begin
temp=id[i-1]
endelse
if i eq 0 then begin
xx=x[i]
yy=y[i]
endif else if id[i] eq temp then begin
xx=[xx,x[i]]
yy=[yy,y[i]]
endif else begin
maskpolyregion,im,ext[i-1],xx,yy
xx=x[i]
yy=y[i]
n=n+1
endelse
end
spawn,'rm -rf /tmp/region.txt'
wrwircamcube, im, fits, newbadpix
END
; MASKREGION
; PURPOSE:
; Produce a new bad pixel mask based on an exist one and ds9 region file
;
; CALLING SEQUENCE:
; MASKREGION, BADPIX = badpix, REGION = ds9.reg, NEWBADPIX= newbadpix
;
; INPUTS:
; BADPIX - string scalar, the name of the original bad pixel mask.
; REGION - DS9 region file, only polygon and box allowed.
; NEWBADPIX - string scalar, the new bad pixel mask.
;
;
; REVISION HISTORY:
; initial version Chi-hung Yan March 2007
;-
PRO RDWIRCAMCUBE, file, im
im=fltarr(2048,2048,4)
for i=0,3 do begin
im[*,*,i]=mrdfits(file,i+1)
endfor
END
;
PRO MASKPOLYREGION, im, ext, x, y
ind=polyfillv(x,y,2048,2048)
tim=im[*,*,ext-1]
tim[ind]=!values.f_nan
im[*,*,ext-1]=tim[*,*]
END
PRO WRWIRCAMCUBE, im, fits, new
main_hd=headfits(fits,exten=0)
mwrfits, 3,new , main_hd, /create
for i=1,4 do begin
hdr = headfits(fits, exten=i)
mwrfits,im[*,*,i-1],new,hdr
endfor
END
PRO MASKREGION, badpix=badpix, region=region, newbadpix=newbadpix
file=region
fits=badpix
; Read image
rdwircamcube,fits,im
spawn,"wircampolyregion "+file+" /tmp/region.txt"
readcol,"/tmp/region.txt",id,ext,x,y
n=0
xx=fltarr(1)
yy=fltarr(1)
for i=0,n_elements(x)-1 do begin
if i eq 0 then begin
temp=id[i]
endif else begin
temp=id[i-1]
endelse
if i eq 0 then begin
xx=x[i]
yy=y[i]
endif else if id[i] eq temp then begin
xx=[xx,x[i]]
yy=[yy,y[i]]
endif else begin
maskpolyregion,im,ext[i-1],xx,yy
xx=x[i]
yy=y[i]
n=n+1
endelse
end
spawn,'rm -rf /tmp/region.txt'
wrwircamcube, im, fits, newbadpix
END
6. Perl Script
#! /usr/bin/perl
#----------------------------------------------------------------------------
# wircampolyregion
#----------------------------------------------------------------------------
# Contents: This script parse the ds9 region file into better format. The out
# out file can be easily used for IDL program.
#----------------------------------------------------------------------------
# Part of: WIRcam C pipeline
# From: ASIAA
# Author: Chi-hung Yan(chyan@asiaa.sinica.edu.tw)
#----------------------------------------------------------------------------
use warnings;
use strict;
# Define global variable
our $SCRIPT_NAME = 'wircampolyregion';
sub trim {
my $string=$_;
for ($string) {
s/^\s+//;
s/\s+$//;
}
return $string;
}
# This finction is used to return program usage.
sub usage(){
print<<EOF
EOF
}
if(@ARGV == 0 || @ARGV != 2){
usage();
exit 0;
}
my $file = $ARGV[0];
my $i=0;
my $j=0;
my $ext;
my $reg=0;
open(OUT,">$ARGV[1]");
open(FILE,$file);
while(<FILE>){
if (/^# tile/){
my @chip=split(" ",$_);
$ext = $chip[2];
}
if (/^polygon/){
$reg=$reg+1;
my @word=split(/[\(\)]/,$_);
my @coord=split(",",$word[1]);
my $n=@coord/2;
for ($i=0;$i<$n;$i++){
printf(OUT "%d %d %7.2f %7.2f\n",$reg,$ext,$coord[2*$i],$coord[2*$i+1]);
}
}
if (/^box/){
$reg=$reg+1;
my @word=split(/[\(\)]/,$_);
my @para=split(",",$word[1]);
printf("$word[1]\n");
printf(OUT "%d %d %7.2f %7.2f\n",$reg,$ext,$para[0]-$para[2]/2,$para[1]-$para[3]/2);
printf(OUT "%d %d %7.2f %7.2f\n",$reg,$ext,$para[0]+$para[2]/2,$para[1]-$para[3]/2);
printf(OUT "%d %d %7.2f %7.2f\n",$reg,$ext,$para[0]+$para[2]/2,$para[1]+$para[3]/2);
printf(OUT "%d %d %7.2f %7.2f\n",$reg,$ext,$para[0]-$para[2]/2,$para[1]+$para[3]/2);
}
}
close(OUT);
#----------------------------------------------------------------------------
# wircampolyregion
#----------------------------------------------------------------------------
# Contents: This script parse the ds9 region file into better format. The out
# out file can be easily used for IDL program.
#----------------------------------------------------------------------------
# Part of: WIRcam C pipeline
# From: ASIAA
# Author: Chi-hung Yan(chyan@asiaa.sinica.edu.tw)
#----------------------------------------------------------------------------
use warnings;
use strict;
# Define global variable
our $SCRIPT_NAME = 'wircampolyregion';
sub trim {
my $string=$_;
for ($string) {
s/^\s+//;
s/\s+$//;
}
return $string;
}
# This finction is used to return program usage.
sub usage(){
print<<EOF
USAGE: $SCRIPT_NAME <ds9.reg> <output.txt>
EXAMPLES:
$SCRIPT_NAME ds9.reg output.txt
EXAMPLES:
$SCRIPT_NAME ds9.reg output.txt
EOF
}
if(@ARGV == 0 || @ARGV != 2){
usage();
exit 0;
}
my $file = $ARGV[0];
my $i=0;
my $j=0;
my $ext;
my $reg=0;
open(OUT,">$ARGV[1]");
open(FILE,$file);
while(<FILE>){
if (/^# tile/){
my @chip=split(" ",$_);
$ext = $chip[2];
}
if (/^polygon/){
$reg=$reg+1;
my @word=split(/[\(\)]/,$_);
my @coord=split(",",$word[1]);
my $n=@coord/2;
for ($i=0;$i<$n;$i++){
printf(OUT "%d %d %7.2f %7.2f\n",$reg,$ext,$coord[2*$i],$coord[2*$i+1]);
}
}
if (/^box/){
$reg=$reg+1;
my @word=split(/[\(\)]/,$_);
my @para=split(",",$word[1]);
printf("$word[1]\n");
printf(OUT "%d %d %7.2f %7.2f\n",$reg,$ext,$para[0]-$para[2]/2,$para[1]-$para[3]/2);
printf(OUT "%d %d %7.2f %7.2f\n",$reg,$ext,$para[0]+$para[2]/2,$para[1]-$para[3]/2);
printf(OUT "%d %d %7.2f %7.2f\n",$reg,$ext,$para[0]+$para[2]/2,$para[1]+$para[3]/2);
printf(OUT "%d %d %7.2f %7.2f\n",$reg,$ext,$para[0]-$para[2]/2,$para[1]+$para[3]/2);
}
}
close(OUT);
