diff --git a/.gitignore b/.gitignore index 8b0e045..6314d3e 100644 --- a/.gitignore +++ b/.gitignore @@ -1,14 +1,148 @@ -# compiled python modules -*.pyc +# Byte-compiled / optimized / DLL files +__pycache__/ +*.py[cod] +*$py.class -# emacs temporary save files -*~ +# C extensions +*.so + +# Distribution / packaging +.Python +build/ +develop-eggs/ +dist/ +downloads/ +eggs/ +.eggs/ +lib/ +lib64/ +parts/ +sdist/ +var/ +wheels/ +share/python-wheels/ +*.egg-info/ +.installed.cfg +*.egg +MANIFEST + +# PyInstaller +# Usually these files are written by a python script from a template +# before PyInstaller builds the exe, so as to inject date/other infos into it. +*.manifest +*.spec + +# Installer logs +pip-log.txt +pip-delete-this-directory.txt + +# Unit test / coverage reports +htmlcov/ +.tox/ +.nox/ +.coverage +.coverage.* +.cache +nosetests.xml +coverage.xml +*.cover +*.py,cover +.hypothesis/ +.pytest_cache/ +cover/ + +# Translations +*.mo +*.pot + +# Django stuff: +*.log +local_settings.py +db.sqlite3 +db.sqlite3-journal + +# Flask stuff: +instance/ +.webassets-cache + +# Scrapy stuff: +.scrapy + +# Sphinx documentation +docs/_build/ + +# PyBuilder +.pybuilder/ +target/ + +# Jupyter Notebook +.ipynb_checkpoints + +# IPython +profile_default/ +ipython_config.py -# setuptools distribution folder -/dist/ +# pyenv +# For a library or package, you might want to ignore these files since the code is +# intended to run in multiple environments; otherwise, check them in: +# .python-version -# python egg metadata, regenerated from source files by setuptools -/*.egg-info +# pipenv +# According to pypa/pipenv#598, it is recommended to include Pipfile.lock in version control. +# However, in case of collaboration, if having platform-specific dependencies or dependencies +# having no cross-platform support, pipenv may install dependencies that don't work, or not +# install all needed dependencies. +#Pipfile.lock -# build folder -/build/ +# poetry +# Similar to Pipfile.lock, it is generally recommended to include poetry.lock in version control. +# This is especially recommended for binary packages to ensure reproducibility, and is more +# commonly ignored for libraries. +# https://python-poetry.org/docs/basic-usage/#commit-your-poetrylock-file-to-version-control +#poetry.lock + +# PEP 582; used by e.g. github.com/David-OConnor/pyflow +__pypackages__/ + +# Celery stuff +celerybeat-schedule +celerybeat.pid + +# SageMath parsed files +*.sage.py + +# Environments +.env +.venv +env/ +venv/ +ENV/ +env.bak/ +venv.bak/ + +# Spyder project settings +.spyderproject +.spyproject + +# Rope project settings +.ropeproject + +# mkdocs documentation +/site + +# mypy +.mypy_cache/ +.dmypy.json +dmypy.json + +# Pyre type checker +.pyre/ + +# pytype static type analyzer +.pytype/ + +# Cython debug symbols +cython_debug/ + +# emacs temporary save files +*~ diff --git a/pyproject.toml b/pyproject.toml new file mode 100644 index 0000000..3b3aee9 --- /dev/null +++ b/pyproject.toml @@ -0,0 +1,2 @@ +[tool.pyright] +typeCheckingMode = "off" diff --git a/xara/core.py b/xara/core.py index dba0bc6..5c1df07 100644 --- a/xara/core.py +++ b/xara/core.py @@ -480,13 +480,13 @@ def find_fourier_origin(img, mykpo, m2pix, bmax=6.0): # Find Fourier plane coordinates which will be considered for the re-centering uv_dist = np.sqrt(mykpo.kpi.UVC[:, 0]**2+mykpo.kpi.UVC[:, 1]**2) - uv_cutoff = np.where(uv_dist < float(r_cutoff))[0] + uv_cutoff = np.where(uv_dist < float(bmax))[0] # Find best sub-pixel shift img_fft = np.fft.rfft2(img_cent) best_xy_shift = leastsq(func=fourier_phase_resid_2d, x0=np.array([0., 0.]), - args=(img_fft, m2pix, uv, uv_cutoff), + args=(img_fft, mykpo, m2pix, uv, uv_cutoff), ftol=1E-1)[0] # Return best shift @@ -514,7 +514,8 @@ def fourier_phase_resid_2d(xy, img_fft, mykpo, m2pix, uv, uv_cutoff): # ========================================================================= -def determine_origin(img, mask=None, algo="BCEN", verbose=True, wmin=10.0): +def determine_origin(img, mask=None, algo="BCEN", verbose=True, wmin=10.0, + mykpo=None, m2pix=None, bmax=None): ''' ------------------------------------------------------------ Determines the origin of the image, using among possible algorithms. @@ -540,6 +541,11 @@ def determine_origin(img, mask=None, algo="BCEN", verbose=True, wmin=10.0): if "cog" in algo.lower(): (x0, y0) = centroid(img1, verbose) + elif "fpnm" in algo.lower(): + if mykpo is None or m2pix is None or bmax is None: + raise ValueError("Need mykpo, m2pix, and bmax parameters") + (x0, y0) = find_fourier_origin(img1, mykpo, m2pix, bmax) + else: (x0, y0) = find_psf_center(img1, verbose, nbit=10, wmin=wmin) @@ -548,7 +554,8 @@ def determine_origin(img, mask=None, algo="BCEN", verbose=True, wmin=10.0): # ========================================================================= def recenter(im0, mask=None, algo="BCEN", subpix=True, between=False, - verbose=True): + verbose=True, return_center=False, dxdy=None, mykpo=None, + m2pix=None, bmax=None): ''' ------------------------------------------------------------ Re-centering algorithm of a 2D image im0 for kernel-analysis @@ -563,6 +570,18 @@ def recenter(im0, mask=None, algo="BCEN", subpix=True, between=False, - subpix: sub-pixel recentering (boolean, default=True) - between: center in between 4 pixels (boolean, default=False) - verbose: display some additional info (boolean, default=True) + - return_center: return center coordinate in original image (boolean, default=False) + - dxdy: center coordinates to use, origin is found with `algo` if None (Tuple[float], default=None) + - mykpo: KPO object, used to determine origin with FPNM (xara.kpo.KPO, default=None) + - m2pix: Meter to pixel scaling parameter for FPNM (float, default=None) + - bmax: Max baseline (in m) for FPNM (float, default=None) + + Returns: + ------- + - img: The recentered image + - If `return_center` is True: + + dx_temp: x coordinate of the center in original image + + dy_temp: y coordinate of the center in original image Remarks: ------- @@ -572,18 +591,25 @@ def recenter(im0, mask=None, algo="BCEN", subpix=True, between=False, ysz, xsz = im0.shape - (x0, y0) = determine_origin(im0, mask=mask, algo=algo, verbose=verbose) + if dxdy is None: + (x0, y0) = determine_origin(im0, mask=mask, algo=algo, verbose=verbose, + mykpo=mykpo, m2pix=m2pix, bmax=bmax) - dy, dx = (y0-ysz/2), (x0-xsz/2) - if between: - dy += 0.5 - dx += 0.5 + dy, dx = (y0-ysz/2), (x0-xsz/2) + if between: + dy += 0.5 + dx += 0.5 - if verbose: - print("centroid: dx=%+5.2f, dy=%+5.2f\n" % (dx, dy), - end="", flush=True) + if verbose: + print("centroid: dx=%+5.2f, dy=%+5.2f\n" % (dx, dy), + end="", flush=True) + + else: + dx, dy = dxdy # integer pixel recentering first + dx_temp = dx + dy_temp = dy im0 = np.roll(np.roll(im0, -int(round(dx)), axis=1), -int(round(dy)), axis=0) @@ -619,7 +645,10 @@ def recenter(im0, mask=None, algo="BCEN", subpix=True, between=False, offset = np.exp(1j*slope) dummy = np.real(shift(ifft(offset * fft(shift(im))))) im0 = dummy[oriy:oriy+ysz, orix:orix+xsz] - return im0 + if not return_center: + return im0 + else: + return im0, dx_temp, dy_temp # ========================================================================= @@ -780,6 +809,63 @@ def uv_phase_regrid_matrix(UVD, UVS, rad): return GG +def hexagon(dim, width, interp_edge=True): + """This function creates a hexagon. + + Adapted from: https://github.com/mikeireland/opticstools/blob/master/opticstools/utils.py#L164 + by Mike Ireland under MIT License + + Parameters + ---------- + dim: int + Size of the 2D array + width: int + flat-to-flat width of the hexagon + + Returns + ------- + pupil: float array (sz,sz) + 2D array hexagonal pupil mask + """ + x = np.arange(dim) - dim / 2.0 + xy = np.meshgrid(x, x) + xx = xy[1] + yy = xy[0] + hex = np.zeros((dim, dim)) + scale = 1.5 + offset = 0.5 + if interp_edge: + # !!! Not fully implemented yet. Need to compute the orthogonal distance + # from each line and accurately find fractional area of each pixel. + hex = ( + np.minimum(np.maximum(width / 2 - yy + offset, 0), 1) + * np.minimum(np.maximum(width / 2 + yy + offset, 0), 1) + * np.minimum( + np.maximum((width - np.sqrt(3) * xx - yy + offset) * scale, 0), 1 + ) + * np.minimum( + np.maximum((width - np.sqrt(3) * xx + yy + offset) * scale, 0), 1 + ) + * np.minimum( + np.maximum((width + np.sqrt(3) * xx - yy + offset) * scale, 0), 1 + ) + * np.minimum( + np.maximum((width + np.sqrt(3) * xx + yy + offset) * scale, 0), 1 + ) + ) + else: + w = np.where( + (yy < width / 2) + * (yy > (-width / 2)) + * (yy < (width - np.sqrt(3) * xx)) + * (yy > (-width + np.sqrt(3) * xx)) + * (yy < (width + np.sqrt(3) * xx)) + * (yy > (-width - np.sqrt(3) * xx)) + ) + hex[w] = 1.0 + return hex + + # ========================================================================= def create_discrete_model(apert, ppscale, step, binary=True, tmin=0.8): '''------------------------------------------------------------------ diff --git a/xara/kpi.py b/xara/kpi.py old mode 100644 new mode 100755 index 5bb3d92..bcab7d1 --- a/xara/kpi.py +++ b/xara/kpi.py @@ -41,7 +41,8 @@ from astropy.io import fits import pickle import gzip -from numpy.linalg import norm + +from .core import hexagon class KPI(object): @@ -57,7 +58,7 @@ class KPI(object): # ========================================================================= def __init__(self, fname=None, array=None, ndgt=5, - bmax=None, ID=""): + bmax=None, hexa=False, ID=""): ''' Default instantiation of a KerPhase_Relation object: ------------------------------------------------------------------- @@ -115,7 +116,7 @@ def __init__(self, fname=None, array=None, ndgt=5, self.load_aperture_model(fname=fname) except OSError: raise Exception("Not a valid coordinate file") - self.rebuild_model(ndgt=ndgt, bmax=bmax) + self.rebuild_model(ndgt=ndgt, bmax=bmax, hexa=hexa) print("KPI data successfully loaded") @@ -123,7 +124,7 @@ def __init__(self, fname=None, array=None, ndgt=5, print("Attempting to build KPI from array") try: self.load_aperture_model(data=array) - self.rebuild_model(ndgt=ndgt, bmax=bmax) + self.rebuild_model(ndgt=ndgt, bmax=bmax, hexa=hexa) self.name = ID except: print("Problem using array %s" % (array,)) @@ -176,7 +177,7 @@ def load_aperture_model(self, fname=None, data=None): # ========================================================================= # ========================================================================= - def rebuild_model(self, ndgt=5, bmax=None): + def rebuild_model(self, ndgt=5, bmax=None, hexa=False): ''' ------------------------------------------------------------------ Builds or rebuilds the Fourier-phase model, from the information provided by the virtual aperture coordinate (VAC) table. @@ -189,6 +190,8 @@ def rebuild_model(self, ndgt=5, bmax=None): - bmax: a baseline filtering parameter (in meters): if not None, baselines larger than bmax are eliminated from the discrete model + - hexa: if True, masks out noisy baselines using a hexagonal shape + ------------------------------------------------------------------ ''' prec = 10**(-ndgt) @@ -240,19 +243,36 @@ def rebuild_model(self, ndgt=5, bmax=None): # 1.5. Special case: baseline filtering # ------------------------------------- if bmax is not None: - uv_sampl = self.UVC.copy() # copy previously identified baselines - # uvm = np.abs(self.UVC).max() # max baseline length + if hexa: + flat_to_flat = 2.*6.5 # m, size is doubled in uv-plane + mask = hexagon(1024, 512*2.*float(bmax)/flat_to_flat, interp_edge=False) # center is (512, 512) + uv_sampl = self.UVC.copy() # copy previously identified baselines + # uvm = np.abs(self.UVC).max() # max baseline length + + xx = uv_sampl[:, 0]/flat_to_flat*512.+512. + yy = uv_sampl[:, 1]/flat_to_flat*512.+512. + + keep = np.array([mask[int(round(yy[ii])), int(round(xx[ii]))] for ii in range(uv_sampl.shape[0])]) > 0.5 + self.UVC = uv_sampl[keep] + self.nbuv = (self.UVC.shape)[0] + + print("%d baselines were preserved after filtering" % (self.nbuv,)) + self.BMAX = bmax + self.BLEN = np.hypot(self.UVC[:, 0], self.UVC[:, 1]) + else: + uv_sampl = self.UVC.copy() # copy previously identified baselines + # uvm = np.abs(self.UVC).max() # max baseline length - blength = np.sqrt(np.abs(uv_sampl[:, 0])**2 + - np.abs(uv_sampl[:, 1])**2) + blength = np.sqrt(np.abs(uv_sampl[:, 0])**2 + + np.abs(uv_sampl[:, 1])**2) - keep = (blength < bmax) - self.UVC = uv_sampl[keep] - self.nbuv = (self.UVC.shape)[0] + keep = (blength < bmax) + self.UVC = uv_sampl[keep] + self.nbuv = (self.UVC.shape)[0] - print("%d baselines were preserved after filtering" % (self.nbuv,)) - self.BMAX = bmax - self.BLEN = np.hypot(self.UVC[:, 0], self.UVC[:, 1]) + print("%d baselines were preserved after filtering" % (self.nbuv,)) + self.BMAX = bmax + self.BLEN = np.hypot(self.UVC[:, 0], self.UVC[:, 1]) # 2. compute baseline mapping model + redundancy # ---------------------------------------------- @@ -345,34 +365,14 @@ def compute_KPM(self,): self.KPM[i, :] = (Vh)[KPhiCol[i], :] self.KPM = self.KPM.dot(np.diag(self.RED)) + norm = np.linalg.norm(self.KPM, axis=1) + self.KPM = np.divide(self.KPM.T, norm).T print('first %d singular values for this array:' % ( min(10, self.nbap-1))) print(np.round(S[:10], 5)) print(self) - self._KPM_norm = norm(self.KPM, axis=1) - self.set_normalize_KPM(True) # new default setting! - - # ========================================================================= - # ========================================================================= - - def set_normalize_KPM(self, setting=True): - ''' ------------------------------------------------------------------ - Normalize the KP matrix. - - Set in the hope that it'll facilitate the performance comparison of - different models and/or with strictly non-redundant scenarios. - - Keyword argument: - setting -- (default True) - ------------------------------------------------------------------ ''' - if setting is True: - self._KPM0 = self.KPM.copy() - self.KPM = np.diag(1/self._KPM_norm).dot(self._KPM0) - else: - self.KPM = self._KPM0.copy() - # ========================================================================= # ========================================================================= @@ -429,7 +429,7 @@ def load_from_fits(self, fname): # ------------------------------------ tmp = hdulist['UV-PLANE'].data self.UVC = np.array([tmp['UUC'], tmp['VVC']]).T - self.RED = np.array(tmp['RED']).astype(np.float) + self.RED = np.array(tmp['RED']).astype(float) self.KPM = hdulist['KER-MAT'].data self.BLM = hdulist['BLM-MAT'].data @@ -505,50 +505,49 @@ def __str__(self): # ========================================================================= # ========================================================================= - def plot_pupil_and_uv(self, xymax=None, figsize=(12, 6), plot_redun=False, - cmap=cm.rainbow, ssize=12, lw=0, alpha=1.0, marker='o'): + def plot_pupil_and_uv(self, xymax=4.0, figsize=(12, 6), plot_redun=False, + cmap=cm.gray, ssize=12, lw=0, alpha=1.0, marker='s'): '''Nice plot of the pupil sampling and matching uv plane. - Parameters: + -------------------------------------------------------------------- + Options: ---------- - - xymax: radius of pupil plot in meters (default=None) + - xymax: radius of pupil plot in meters (default=4.0) - figsize: matplotlib figure size (default=(12,6)) - plot_redun: bool add the redundancy information (default=False) - - cmap: matplotlib colormap (default:cm.rainbow) + - cmap: matplotlib colormap (default:cm.gray) - ssize: symbol size (default=12) - lw: line width for symbol outline (default=0) - alpha: gamma (transparency) (default=1) - - maker: matplotlib marker for sub-aperture (default='o') + - maker: matplotlib marker for sub-aperture (default='s') - ------------------------------------------------------------------- ''' - if xymax is not None: - xym = xymax - else: - xym = np.ceil(np.max(self.VAC[:, :1]) * 2) / 2 f0 = plt.figure(figsize=figsize) plt.clf() ax0 = plt.subplot(121) s1, s2 = ssize**2, (ssize/2)**2 - ax0.scatter(self.VAC[:, 0], self.VAC[:, 1], s=s1, c=self.VAC[:, 2], - cmap=cmap, alpha=alpha, marker=marker, lw=lw) - ax0.axis([-xym, xym, -xym, xym]) + p0 = ax0.scatter(self.VAC[:, 0], self.VAC[:, 1], s=s1, c=self.VAC[:, 2], + cmap=cmap, alpha=alpha, marker=marker, lw=lw) + ax0.axis([-xymax, xymax, -xymax, xymax]) ax0.set_aspect('equal') ax0.set_xlabel("Aperture x-coordinate (meters)") ax0.set_ylabel("Aperture y-coordinate (meters)") + plt.colorbar(p0, ax=ax0) ax1 = plt.subplot(122) - ax1.scatter(-self.UVC[:, 0], -self.UVC[:, 1], s=s2, c=self.RED, + p1 = ax1.scatter(-self.UVC[:, 0], -self.UVC[:, 1], s=s2, c=self.RED, cmap=cmap, alpha=alpha, marker=marker, lw=lw) ax1.scatter(self.UVC[:, 0], self.UVC[:, 1], s=s2, c=self.RED, cmap=cmap, alpha=alpha, marker=marker, lw=lw) - ax1.axis([-2*xym, 2*xym, -2*xym, 2*xym]) + ax1.axis([-2*xymax, 2*xymax, -2*xymax, 2*xymax]) ax1.set_aspect('equal') ax1.set_xlabel("Fourier u-coordinate (meters)") ax1.set_ylabel("Fourier v-coordinate (meters)") + plt.colorbar(p1, ax=ax1) # complete previous plot with redundancy of the baseline # ------------------------------------------------------- @@ -593,7 +592,7 @@ def package_as_fits(self, fname=None): hdr['SOFTWARE'] = 'XARA' hdr['KPI-ID'] = self.name[:8] hdr['GRID'] = (False, "True for integer grid mode") - hdr['G-STEP'] = (0.0, "Used for integer grid mode") + hdr['G-STEP'] = (0.0, "Used for integer grid mode") hdr.add_comment("File created by the XARA python pipeline") try: _ = self.BMAX @@ -624,7 +623,7 @@ def package_as_fits(self, fname=None): # KER-MAT HDU # ----------- - kpm_hdu = fits.ImageHDU(self.KerPhi) + kpm_hdu = fits.ImageHDU(self.KPM) kpm_hdu.header.add_comment("Kernel-phase Matrix") kpm_hdu.header['EXTNAME'] = 'KER-MAT' diff --git a/xara/kpo.py b/xara/kpo.py index 548d62a..b353375 100644 --- a/xara/kpo.py +++ b/xara/kpo.py @@ -28,8 +28,8 @@ from astropy.time import Time import copy -from . import core -from . import kpi +from xara import core +from xara import kpi shift = np.fft.fftshift fft = np.fft.fft2 @@ -47,7 +47,7 @@ class KPO(): set. ------------------------------------------------------------------- ''' - def __init__(self, fname=None, array=None, ndgt=5, bmax=None, ID=""): + def __init__(self, fname=None, array=None, ndgt=5, bmax=None, hexa=False, ID=""): ''' Default instantiation of a KerPhase_Relation object: ------------------------------------------------------------------- @@ -57,7 +57,7 @@ def __init__(self, fname=None, array=None, ndgt=5, bmax=None, ID=""): # Default instantiation. self.kpi = kpi.KPI(fname=fname, array=array, - ndgt=ndgt, bmax=bmax, ID=ID) + ndgt=ndgt, bmax=bmax, hexa=hexa, ID=ID) self.TARGET = [] # source names self.CVIS = [] # complex visibilities @@ -449,8 +449,10 @@ def extract_KPD_single_cube(self, cube, pscale, cwavel, target=None, # ========================================================================= # ========================================================================= def extract_KPD_single_frame(self, frame, pscale, cwavel, target=None, - recenter=False, wrad=None, method="LDFT1"): - """Handles the kernel processing of a single square frame + recenter=False, wrad=None, method="LDFT1", + algo_cent="BCEN", bmax_cent=None): + """ ---------------------------------------------------------------- + Handles the kernel processing of a single square frame Parameters ---------- @@ -480,26 +482,29 @@ def extract_KPD_single_frame(self, frame, pscale, cwavel, target=None, self.sgmask = core.super_gauss(isz, isz, wrad) img = frame.copy() - if recenter is True: + if recenter: ysz, xsz = frame.shape (x0, y0) = core.determine_origin(img, mask=self.sgmask, - algo="BCEN", verbose=False) + algo=algo_cent, verbose=False, + mykpo=self, m2pix=m2pix, + bmax=bmax_cent) dy, dx = (y0-ysz/2), (x0-xsz/2) + img = np.roll(np.roll(img, -int(round(dx)), axis=1), + -int(round(dy)), axis=0) + dx_temp = dx - int(round(dx)) + dy_temp = dy - int(round(dy)) + if self.sgmask is not None: # use apodization mask before extraction - if recenter is True: - img *= np.roll(np.roll(self.sgmask, int(round(dx)), axis=1), - int(round(dy)), axis=0) - else: - img *= self.sgmask + img *= self.sgmask # ----- complex visibility extraction ----- temp = self.extract_cvis_from_img(img, m2pix, method) # ---- sub-pixel recentering correction ----- - if recenter is True: + if recenter: uvc = self.kpi.UVC * self.M2PIX - corr = np.exp(i2pi * uvc.dot(np.array([dx, dy])/float(ysz))) + corr = np.exp(i2pi * uvc.dot(np.array([dx_temp, dy_temp])/float(ysz))) temp *= corr cvis.append(temp) @@ -519,7 +524,10 @@ def extract_KPD_single_frame(self, frame, pscale, cwavel, target=None, self.KPDT.append(np.array(kpdata)) self.DETPA.append(np.array(detpa).flatten()) self.MJDATE.append(np.array(mjdate)) - return + if recenter: + return dx, dy + else: + return # ========================================================================= # ========================================================================= @@ -594,7 +602,7 @@ def save_as_fits(self, fname): # KP-DATA HDU # ----------- - kpd_hdu = fits.ImageHDU(self.KPDT[ii].astype(np.float64)) + kpd_hdu = fits.ImageHDU(self.KPDT[ii].astype(float)) kpd_hdu.header.add_comment("Kernel-phase data") kpd_hdu.header['EXTNAME'] = 'KP-DATA%d' % (ii+1,) kpd_hdu.header['TARGET'] = self.TARGET[ii] @@ -626,7 +634,7 @@ def save_as_fits(self, fname): try: _ = self.kp_cov - kcv_hdu = fits.ImageHDU(self.kp_cov.astype(np.float64)) + kcv_hdu = fits.ImageHDU(self.kp_cov.astype(float)) kcv_hdu.header.add_comment( "Kernel-phase covariance matrix") kcv_hdu.header['EXTNAME'] = 'KP-COV' @@ -958,7 +966,7 @@ def scatter_uv_map(self, data, sym=True, title="", cbar=True, fsize=5, data2 = np.append(data, data) if sym else np.append(data, -data) ssz = ssize**2 # symbol size - dtitle = "Fourier map" if title is "" else title + dtitle = "Fourier map" if title == "" else title if not cbar: fig, ax = plt.subplots()