add python scripts

3 files changed, 428 insertions, 0 deletions

diff --git a/python/radialPSD2d.py b/python/radialPSD2d.py
new file mode 100755
index 0000000..cc5dd85
--- /dev/null
+++ b/python/radialPSD2d.py
@@ -0,0 +1,173 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+#
+# Aaron LI <aaronly.me@gmail.com>
+# 2015/04/22
+#
+
+"""
+Computes the radially averaged power spectral density (power spectrum).
+"""
+
+
+import numpy as np
+from scipy import fftpack
+
+
+def PSD2d( img, normalize=True ):
+    """
+    Computes the 2D power spectrum of the given image.
+
+    Reference:
+    [1] raPsd2d.m by Evan Ruzanski
+        Radially averaged power spectrum of 2D real-valued matrix
+        https://www.mathworks.com/matlabcentral/fileexchange/23636-radially-averaged-power-spectrum-of-2d-real-valued-matrix
+    """
+    img = np.array( img )
+    rows, cols = img.shape
+    ## Compute power spectrum
+    # Perform the Fourier transform and shift the zero-frequency
+    # component to the center of the spectrum.
+    imgf = fftpack.fftshift( fftpack.fft2( img ) )
+    if normalize:
+        norm = rows * cols
+    else:
+        norm = 1.0  # Do not normalize
+    psd2d = ( np.abs( imgf ) / norm ) ** 2
+    return psd2d
+
+
+def radialPSD( psd2d ):
+    """
+    Computes the radially averaged power spectral density (power spectrum)
+    from the provided 2D power spectrum.
+
+    Reference:
+    [1] raPsd2d.m by Evan Ruzanski
+        Radially averaged power spectrum of 2D real-valued matrix
+        https://www.mathworks.com/matlabcentral/fileexchange/23636-radially-averaged-power-spectrum-of-2d-real-valued-matrix
+    """
+    psd2d = np.array( psd2d )
+    rows, cols = psd2d.shape
+    ## Adjust the PSD array size
+    dim_diff = np.abs( rows - cols )
+    dim_max = max( rows, cols )
+    # Pad the PSD array to be sqaure
+    if rows > cols:
+        # pad columns
+        if np.mod( dim_diff, 2 ) == 0:
+            cols_left = np.zeros( (rows, dim_diff/2) )
+            cols_left[:] = np.nan
+            cols_right = np.zeros( (rows, dim_diff/2) )
+            cols_right[:] = np.nan
+            psd2d = np.hstack( (cols_left, psd2d, cols_right) )
+        else:
+            cols_left = np.zeros( (rows, np.floor(dim_diff/2)) )
+            cols_left[:] = np.nan
+            cols_right = np.zeros( (rows, np.floor(dim_diff/2)+1) )
+            cols_right[:] = np.nan
+            psd2d = np.hstack( (cols_left, psd2d, cols_right) )
+    elif rows < cols:
+        # pad rows
+        if np.mod( dim_diff, 2 ) == 0:
+            rows_top = np.zeros( (dim_diff/2, cols) )
+            rows_top[:] = np.nan
+            rows_bottom = np.zeros( (dim_diff/2, cols) )
+            rows_bottom[:] = np.nan
+            psd2d = np.vstack( (rows_top, psd2d, rows_bottom) )
+        else:
+            rows_top = np.zeros( (np.floor(dim_diff/2), cols) )
+            rows_top[:] = np.nan
+            rows_bottom = np.zeros( (np.floor(dim_diff/2)+1, cols) )
+            rows_bottom[:] = np.nan
+            psd2d = np.vstack( (rows_top, psd2d, rows_bottom) )
+    ## Compute radially average power spectrum
+    px = np.arange( -dim_max/2, dim_max/2 )
+    x, y = np.meshgrid( px, px )
+    rho, phi = cart2pol( x, y )
+    rho = np.around( rho ).astype(int)
+    dim_half = np.floor( dim_max/2 ) + 1
+    radial_psd = np.zeros( dim_half )
+    radial_psd_err = np.zeros( dim_half ) # standard error
+    for r in np.arange( dim_half, dtype=int ):
+        # Get the indices of the elements satisfying rho[i,j]==r
+        ii, jj = (rho == r).nonzero()
+        # Calculate the mean value at a given radii
+        data = psd2d[ii, jj]
+        radial_psd[r] = np.nanmean( data )
+        radial_psd_err[r] = np.nanstd( data )
+    # Calculate frequencies
+    f = fftpack.fftfreq( dim_max, d=1 ) # sample spacing: set to 1 pixel
+    freqs = np.abs( f[:dim_half] )
+    #
+    return (freqs, radial_psd, radial_psd_err)
+
+
+def plotRadialPSD( freqs, radial_psd, radial_psd_err=None ):
+    """
+    Make a plot of the radial 1D PSD with matplotlib.
+    """
+    try:
+        import matplotlib.pyplot as plt
+    except ImportError:
+        import sys
+        print( "Error: matplotlib.pyplot cannot be imported!",
+                file=sys.stderr )
+        sys.exit( 30 )
+    dim_half = radial_psd.size
+    # plot
+    plt.figure()
+    plt.loglog( freqs, radial_psd )
+    plt.title( "Radially averaged power spectrum" )
+    plt.xlabel( "k (/pixel)" )
+    plt.ylabel( "Power" )
+    plt.show()
+
+
+def cart2pol( x, y ):
+    """
+    Convert Cartesian coordinates to polar coordinates.
+    """
+    rho = np.sqrt( x**2 + y**2 )
+    phi = np.arctan2( y, x )
+    return (rho, phi)
+
+def pol2cart( rho, phi ):
+    """
+    Convert polar coordinates to Cartesian coordinates.
+    """
+    x = rho * np.cos( phi )
+    y = rho * np.sin( phi )
+    return (x, y)
+
+
+def loadData( filename, ftype="fits" ):
+    """
+    Load data from file into numpy array.
+    """
+    if ftype == "fits":
+        try:
+            from astropy.io import fits
+        except ImportError:
+            import sys
+            print( "Error: astropy.io.fits cannot be imported!",
+                    file=sys.stderr )
+            sys.exit( 20 )
+        ffile = fits.open( filename )
+        data = ffile[0].data.astype( float )
+        ffile.close()
+    else:
+        print( "Error: not implemented yet!",
+                file=sys.stderr )
+        sys.exit( 10 )
+    #
+    return data
+
+
+def main():
+    pass
+
+
+if __name__ == "__main__":
+    main()
+
diff --git a/python/splitBoxRegion.py b/python/splitBoxRegion.py
new file mode 100755
index 0000000..5254686
--- /dev/null
+++ b/python/splitBoxRegion.py
@@ -0,0 +1,148 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8-
+#
+# Split the strip-shaped CCD gaps regions into a series of small
+# square regions, which are used as the input regions of 'roi' to
+# determine the corresponding background regions, and finally providied
+# to 'dmfilth' in order to fill in the CCD gaps.
+#
+# Aaron LI
+# 2015/08/12
+#
+# Changelogs:
+# v0.1.0, 2015/08/12
+#   * initial version
+#
+
+
+__version__ = "0.1.0"
+__date__    = "2015/08/12"
+
+
+import os
+import sys
+import re
+import math
+import argparse
+from io import TextIOWrapper
+
+
+## BoxRegion {{{
+class BoxRegion(object):
+    """
+    CIAO/DS9 "rotbox"/"box" region class.
+
+    rotbox/box format:
+        rotbox(xc, yc, width, height, rotation)
+        box(xc, yc, width, height, rotation)
+    Notes:
+        rotation: [0, 360) (degree)
+    """
+    def __init__(self, xc=None, yc=None,
+                 width=None, height=None, rotation=None):
+        self.regtype  = "rotbox"
+        self.xc       = xc
+        self.yc       = yc
+        self.width    = width
+        self.height   = height
+        self.rotation = rotation
+
+    def __str__(self):
+        return "%s(%s,%s,%s,%s,%s)" % (self.regtype, self.xc, self.yc,
+                self.width, self.height, self.rotation)
+
+    @classmethod
+    def parse(cls, regstr):
+        """
+        Parse region string.
+        """
+        regex_box = re.compile(r'^\s*(box|rotbox)\(([0-9. ]+),([0-9. ]+),([0-9. ]+),([0-9. ]+),([0-9. ]+)\)\s*$', re.I)
+        m = regex_box.match(regstr)
+        if m:
+            regtype  = m.group(1)
+            xc       = float(m.group(2))
+            yc       = float(m.group(3))
+            width    = float(m.group(4))
+            height   = float(m.group(5))
+            rotation = float(m.group(6))
+            return cls(xc, yc, width, height, rotation)
+        else:
+            return None
+
+    def split(self, filename=None):
+        """
+        Split strip-shaped box region into a series small square regions.
+        """
+        angle = self.rotation * math.pi / 180.0
+        # to record the center coordinates of each split region
+        centers = []
+        if self.width > self.height:
+            # number of regions after split
+            nreg = math.ceil(self.width / self.height)
+            # width & height of the split region
+            width = self.width / nreg
+            height = self.height
+            # position of the left-most region
+            x_l = self.xc - 0.5*self.width * math.cos(angle)
+            y_l = self.yc - 0.5*self.width * math.sin(angle)
+            for i in range(nreg):
+                x = x_l + (0.5 + i) * width * math.cos(angle)
+                y = y_l + (0.5 + i) * width * math.sin(angle)
+                centers.append((x, y))
+        else:
+            # number of regions after split
+            nreg = math.ceil(self.height / self.width)
+            # width & height of the split region
+            width = self.width
+            height = self.height / nreg
+            # position of the left-most region
+            x_l = self.xc + 0.5*self.height * math.cos(angle + math.pi/2)
+            y_l = self.yc + 0.5*self.height * math.sin(angle + math.pi/2)
+            for i in range(nreg):
+                x = x_l - (0.5 + i) * height * math.cos(angle + math.pi/2)
+                y = y_l - (0.5 + i) * height * math.sin(angle + math.pi/2)
+                centers.append((x, y))
+        # create split regions
+        regions = []
+        for (x, y) in centers:
+            regions.append(self.__class__(x, y, width+2, height+2,
+                                          self.rotation))
+        # write split regions into file if specified
+        if isinstance(filename, str):
+            regout = open(filename, "w")
+            regout.write("\n".join(map(str, regions)))
+            regout.close()
+        else:
+            return regions
+
+## BoxRegion }}}
+
+
+def main():
+    # command line arguments
+    parser = argparse.ArgumentParser(
+            description="Split strip-shaped rotbox region into " + \
+                    "a series of small square regions.",
+            epilog="Version: %s (%s)" % (__version__, __date__))
+    parser.add_argument("-V", "--version", action="version",
+            version="%(prog)s " + "%s (%s)" % (__version__, __date__))
+    parser.add_argument("infile", help="input rotbox region file")
+    parser.add_argument("outfile", help="output file of the split regions")
+    args = parser.parse_args()
+
+    outfile = open(args.outfile, "w")
+    regex_box = re.compile(r'^\s*(box|rotbox)\([0-9., ]+\)\s*$', re.I)
+    for line in open(args.infile, "r"):
+        if regex_box.match(line):
+            reg = BoxRegion.parse(line)
+            split_regs = reg.split()
+            outfile.write("\n".join(map(str, split_regs)) + "\n")
+        else:
+            outfile.write(line)
+
+    outfile.close()
+
+
+if __name__ == "__main__":
+    main()
+
diff --git a/python/splitCCDgaps.py b/python/splitCCDgaps.py
new file mode 100644
index 0000000..bc26c29
--- /dev/null
+++ b/python/splitCCDgaps.py
@@ -0,0 +1,107 @@
+# -*- coding: utf-8-
+#
+# Aaron LI
+# 2015/08/12
+#
+
+"""
+Split the long-strip-shaped CCD gaps regions into a series of little
+square regions, which are used as the input regions of 'roi' to
+determine the corresponding background regions, and finally providied
+to 'dmfilth' in order to fill in the CCD gaps.
+"""
+
+
+import re
+import math
+from io import TextIOWrapper
+
+
+class BoxRegion(object):
+    """
+    CIAO/DS9 "rotbox"/"box" region class.
+
+    rotbox/box format:
+        rotbox(xc, yc, width, height, rotation)
+        box(xc, yc, width, height, rotation)
+    Notes:
+        rotation: [0, 360) (degree)
+    """
+    def __init__(self, xc=None, yc=None,
+                 width=None, height=None, rotation=None):
+        self.regtype  = "rotbox"
+        self.xc       = xc
+        self.yc       = yc
+        self.width    = width
+        self.height   = height
+        self.rotation = rotation
+
+    def __str__(self):
+        return "%s(%s,%s,%s,%s,%s)" % (self.regtype, self.xc, self.yc,
+                self.width, self.height, self.rotation)
+
+    @classmethod
+    def parse(cls, regstr):
+        """
+        Parse region string.
+        """
+        regex_box = re.compile(r'^(box|rotbox)\(([0-9.]+),([0-9.]+),([0-9.]+),([0-9.]+),([0-9.]+)\)$', re.I)
+        m = regex_box.match(regstr)
+        if m:
+            regtype  = m.group(1)
+            xc       = float(m.group(2))
+            yc       = float(m.group(3))
+            width    = float(m.group(4))
+            height   = float(m.group(5))
+            rotation = float(m.group(6))
+            return cls(xc, yc, width, height, rotation)
+        else:
+            return None
+
+    def split(self, filename=None):
+        """
+        Split long-strip-shaped box region into a series square box regions.
+        """
+        angle = self.rotation * math.pi / 180.0
+        # to record the center coordinates of each split region
+        centers = []
+        if self.width > self.height:
+            # number of regions after split
+            nreg = math.ceil(self.width / self.height)
+            # width & height of the split region
+            width = self.width / nreg
+            height = self.height
+            # position of the left-most region
+            x_l = self.xc - 0.5*self.width * math.cos(angle)
+            y_l = self.yc - 0.5*self.width * math.sin(angle)
+            for i in range(nreg):
+                x = x_l + (0.5 + i) * width * math.cos(angle)
+                y = y_l + (0.5 + i) * width * math.sin(angle)
+                centers.append((x, y))
+        else:
+            # number of regions after split
+            nreg = math.ceil(self.height / self.width)
+            # width & height of the split region
+            width = self.width
+            height = self.height / nreg
+            # position of the left-most region
+            x_l = self.xc + 0.5*self.height * math.cos(angle + math.pi/2)
+            y_l = self.yc + 0.5*self.height * math.sin(angle + math.pi/2)
+            for i in range(nreg):
+                x = x_l - (0.5 + i) * height * math.cos(angle + math.pi/2)
+                y = y_l - (0.5 + i) * height * math.sin(angle + math.pi/2)
+                centers.append((x, y))
+        # create split regions
+        regions = []
+        for (x, y) in centers:
+            regions.append(self.__class__(x, y, width+2, height+2,
+                                          self.rotation))
+        # write split regions into file if specified
+        if isinstance(filename, str):
+            regout = open(filename, "w")
+            regout.write("\n".join(map(str, regions)))
+            regout.close()
+        else:
+            return regions
+
+