【image】Bad Pixel Mask using DS9 Region Files

The basic idea to create the mask image is simple. Firstly, you need to define the bad pixels with ds9 region like box and polygon (one can only use these two type regions only usually). Usually,  there is two main steps to realize the idea.
(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 numpymatplotlib.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

[From: http://stackoverflow.com/questions/3654289/scipy-create-2d-polygon-mask]

(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    

=======================================================================

Two example to realize the idea:

(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

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

 USAGE: $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);