diff options
| -rwxr-xr-x | python/radialPSD2d.py | 173 | ||||
| -rwxr-xr-x | python/splitBoxRegion.py | 148 | ||||
| -rw-r--r-- | python/splitCCDgaps.py | 107 |
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 + + |