aboutsummaryrefslogtreecommitdiffstats
path: root/python
diff options
context:
space:
mode:
Diffstat (limited to 'python')
-rwxr-xr-xpython/radialPSD2d.py173
-rwxr-xr-xpython/splitBoxRegion.py148
-rw-r--r--python/splitCCDgaps.py107
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
+
+