From 1b82cd0f82decdcd943c85b2f323d489f076d608 Mon Sep 17 00:00:00 2001 From: Carlos Gomez Date: Fri, 9 Dec 2016 13:21:24 +0100 Subject: [PATCH 1/3] no message --- vip/pca/pca_fullfr.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/vip/pca/pca_fullfr.py b/vip/pca/pca_fullfr.py index a1368c94..5aa591ba 100644 --- a/vip/pca/pca_fullfr.py +++ b/vip/pca/pca_fullfr.py @@ -210,7 +210,8 @@ def subtract_projection(cube, cube_ref, ncomp, scaling, mask_center_px, if not cube.ndim>2: raise TypeError('Input array is not a 3d or 4d array') if angle_list is not None: - if (cube.ndim==3 and not (cube.shape[0] == angle_list.shape[0])) or (cube.ndim==4 and not (cube.shape[1] == angle_list.shape[0])): + if (cube.ndim==3 and not (cube.shape[0] == angle_list.shape[0])) or \ + (cube.ndim==4 and not (cube.shape[1] == angle_list.shape[0])): msg = "Angle list vector has wrong length. It must equal the number" msg += " frames in the cube." raise TypeError(msg) From 4d140f5a1a30c119ab349dada56dabe80985d27b Mon Sep 17 00:00:00 2001 From: Carlos Gomez Date: Fri, 9 Dec 2016 15:44:24 +0100 Subject: [PATCH 2/3] VIP v0.6.0 changes: - Included PCA based sky subtraction - Fixed bugs in verify_fits, pca_optimize_snr and contrast_curve functions --- vip/__init__.py | 2 +- vip/calib/__init__.py | 3 + vip/calib/skysubtraction.py | 130 ++++++++++++++++++++++++++++++++++++ vip/fits/fits.py | 41 +++++++----- vip/pca/pca_fullfr.py | 2 +- vip/phot/contrcurve.py | 7 +- 6 files changed, 162 insertions(+), 23 deletions(-) create mode 100644 vip/calib/skysubtraction.py diff --git a/vip/__init__.py b/vip/__init__.py index 376e592f..1f58a80e 100644 --- a/vip/__init__.py +++ b/vip/__init__.py @@ -12,7 +12,7 @@ from . import stats from . import var -__version__ = "0.5.5" +__version__ = "0.6.0" print("---------------------------------------------------") print(" oooooo oooo ooooo ooooooooo. ") diff --git a/vip/calib/__init__.py b/vip/calib/__init__.py index c094e608..b73da09a 100644 --- a/vip/calib/__init__.py +++ b/vip/calib/__init__.py @@ -17,6 +17,7 @@ - DFT upsampling or fourier cross-correlation (Guizar et al. 2008), - radon transform for broadband frames (Pueyo et al. 2014), - using satellite/waffle spots (fitting plus intersection). + - sky subtraction (PCA method). Astronomical calibration functionality like flat fielding and dark-sky subtraction, in spite of its simplicity was not included in VIP because of the @@ -33,5 +34,7 @@ from parangles import * from recentering import * from rescaling import * +from skysubtraction import * from subsampling import * + diff --git a/vip/calib/skysubtraction.py b/vip/calib/skysubtraction.py new file mode 100644 index 00000000..1745dd73 --- /dev/null +++ b/vip/calib/skysubtraction.py @@ -0,0 +1,130 @@ +#! /usr/bin/env python + +""" +Module with sky subtraction functionalities. +""" + +from __future__ import division + +__author__ = 'C. Gomez @ ULg' +__all__ = ['cube_subtract_sky_pca'] + +import numpy as np +from ..pca import prepare_matrix, svd_wrapper + + +def cube_subtract_sky_pca(sci_cube, sky_cube, mask, ref_cube=None, ncomp=2): + """ PCA based sky subtraction. + + Notes + ----- + MSF : masked science frame + MPC_1,...,MPC_k : PCs masked (the same way) + MPC : matrix whose columns are MPC_i + + The coefficients c_1,...,c_k are obtained by solving in the least square + sense. + MSF = sum_i(c_i * MPC_i) + + the solution vector c = (c_1,...,c_k)' is given by + c = inv(MPC' * MPC) * MPC' * MSF, + where MSF is in vector column form, ' denotes the matrix transpose and * the + matrix product. + + Note that MPC' * MSF is equal to PC' * MSF, but the masked PCs are not + orthonormal, hence MPC' * MPC is not the identity, therefore + inv(MPC' * MPC) * MPC' * MSF does not reduce to PC' * MSF. + + Parameters + ---------- + sci_cube : array_like + 3d array of science frames. + sky_cube : array_like + 3d array of sky frames. + mask : array_like + Mask indicating the region for the analysis. Can be created with the + function vip.var.create_ringed_spider_mask. + ref_cube : array_like or None + Reference cube. + ncomp : int + Sets the number of PCs you want to use in the sky subtraction. + + Returns + ------- + Sky subtracted cube. + + """ + if sci_cube.shape[1] != sky_cube.shape[1] or sci_cube.shape[2] != \ + sky_cube.shape[2]: + raise TypeError('Science and Sky frames sizes do not match') + if ref_cube is not None: + if sci_cube.shape[1] != ref_cube.shape[1] or sci_cube.shape[2] != \ + ref_cube.shape[2]: + raise TypeError('Science and Reference frames sizes do not match') + + # Getting the EVs from the sky cube + Msky = prepare_matrix(sky_cube, scaling=None, verbose=False) + sky_pcs = svd_wrapper(Msky, 'lapack', sky_cube.shape[0], False, + False) + sky_pcs_cube = sky_pcs.reshape(sky_cube.shape[0], sky_cube.shape[1], + sky_cube.shape[1]) + + # Masking the science cube + sci_cube_masked = np.zeros_like(sci_cube) + ind_masked = np.where(mask == 0) + for i in xrange(sci_cube.shape[0]): + masked_image = np.copy(sci_cube[i]) + masked_image[ind_masked] = 0 + sci_cube_masked[i] = masked_image + Msci_masked = prepare_matrix(sci_cube_masked, scaling=None, + verbose=False) + + # Masking the PCs learned from the skies + sky_pcs_cube_masked = np.zeros_like(sky_pcs_cube) + for i in xrange(sky_pcs_cube.shape[0]): + masked_image = np.copy(sky_pcs_cube[i]) + masked_image[ind_masked] = 0 + sky_pcs_cube_masked[i] = masked_image + + # Project the masked frames onto the sky PCs to get the coefficients + transf_sci = np.zeros((sky_cube.shape[0], Msci_masked.shape[0])) + for i in xrange(Msci_masked.shape[0]): + transf_sci[:, i] = np.inner(sky_pcs, Msci_masked[i].T) + + Msky_pcs_masked = prepare_matrix(sky_pcs_cube_masked, scaling=None, + verbose=False) + mat_inv = np.linalg.inv(np.dot(Msky_pcs_masked, Msky_pcs_masked.T)) + transf_sci_scaled = np.dot(mat_inv, transf_sci) + + # Obtaining the optimized sky and subtraction + sci_cube_skysub = np.zeros_like(sci_cube) + for i in xrange(Msci_masked.shape[0]): + sky_opt = np.array([np.sum( + transf_sci_scaled[j, i] * sky_pcs_cube[j] for j in range(ncomp))]) + sci_cube_skysub[i] = sci_cube[i] - sky_opt + + # Processing the reference cube (if any) + if ref_cube is not None: + ref_cube_masked = np.zeros_like(ref_cube) + for i in xrange(ref_cube.shape[0]): + masked_image = np.copy(ref_cube[i]) + masked_image[ind_masked] = 0 + ref_cube_masked[i] = masked_image + Mref_masked = prepare_matrix(ref_cube_masked, scaling=None, + verbose=False) + transf_ref = np.zeros((sky_cube.shape[0], Mref_masked.shape[0])) + for i in xrange(Mref_masked.shape[0]): + transf_ref[:, i] = np.inner(sky_pcs, Mref_masked[i].T) + + transf_ref_scaled = np.dot(mat_inv, transf_ref) + + ref_cube_skysub = np.zeros_like(ref_cube) + for i in xrange(Mref_masked.shape[0]): + sky_opt = np.array([np.sum( + transf_ref_scaled[j, i] * sky_pcs_cube[j] for j in + range(ncomp))]) + ref_cube_skysub[i] = ref_cube[i] - sky_opt + + return sci_cube_skysub, ref_cube_skysub + else: + return sci_cube_skysub diff --git a/vip/fits/fits.py b/vip/fits/fits.py index 8032e403..3970968f 100644 --- a/vip/fits/fits.py +++ b/vip/fits/fits.py @@ -15,7 +15,7 @@ import os -from astropy.io import fits +from astropy.io import fits as ap_fits def open_fits(fitsfilename, n=0, header=False, ignore_missing_end=False, @@ -45,7 +45,7 @@ def open_fits(fitsfilename, n=0, header=False, ignore_missing_end=False, """ if not fitsfilename.endswith('.fits'): fitsfilename = str(fitsfilename+'.fits') - hdulist = fits.open(fitsfilename, memmap=True, + hdulist = ap_fits.open(fitsfilename, memmap=True, ignore_missing_end=ignore_missing_end) data = hdulist[n].data #if data.dtype.name == 'float64': data = np.array(data, dtype='float32') @@ -116,7 +116,7 @@ def open_adicube(fitsfilename, verbose=True): """ if not fitsfilename.endswith('.fits'): fitsfilename = str(fitsfilename+'.fits') - hdulist = fits.open(fitsfilename, memmap=True) + hdulist = ap_fits.open(fitsfilename, memmap=True) data = hdulist[0].data if not data.ndim ==3: raise TypeError('Input fits file does not contain a cube or 3d array.') @@ -143,11 +143,6 @@ def byteswap_array(array): byte order array. Some operations require the data to be byteswaped before and will complain about it. This function will help in this cases. - - Notes - ----- - http://docs.scipy.org/doc/numpy-1.10.1/user/basics.byteswapping.html - Parameters ---------- array : array_like @@ -157,6 +152,12 @@ def byteswap_array(array): ------- array_out : array_like 2d resulting array after the byteswap operation. + + + Notes + ----- + http://docs.scipy.org/doc/numpy-1.10.1/user/basics.byteswapping.html + """ array_out = array.byteswap().newbyteorder() return array_out @@ -165,19 +166,25 @@ def byteswap_array(array): def info_fits(fitsfilename): """Prints the information about a fits file. """ - hdulist = fits.open(fitsfilename, memmap=True) + hdulist = ap_fits.open(fitsfilename, memmap=True) hdulist.info() -def verify_fits(fits): +def verify_fits(fitspath): """Verifies "the FITS standard" of a fits file or list of fits. + + Parameters + ---------- + fitspath : string + Path to the fits file or list with fits filename paths. + """ - if isinstance(fits, list): - for ffile in fits: - f = fits.open(ffile) + if isinstance(fitspath, list): + for ffile in fitspath: + f = ap_fits.open(ffile) f.verify() else: - f = fits.open(fits) + f = ap_fits.open(fitspath) f.verify() @@ -205,9 +212,9 @@ def write_fits(filename, array, header=None, dtype32=True, verbose=True): os.remove(filename) if verbose: print "\nFits file successfully overwritten" - fits.writeto(filename, array, header) + ap_fits.writeto(filename, array, header) else: - fits.writeto(filename, array, header) + ap_fits.writeto(filename, array, header) if verbose: print "\nFits file successfully saved" @@ -215,7 +222,7 @@ def write_fits(filename, array, header=None, dtype32=True, verbose=True): def append_extension(filename, array): """Appends an extension to fits file. """ - fits.append(filename, array) + ap_fits.append(filename, array) print "\nFits extension appended" diff --git a/vip/pca/pca_fullfr.py b/vip/pca/pca_fullfr.py index 5aa591ba..51dc2bf0 100644 --- a/vip/pca/pca_fullfr.py +++ b/vip/pca/pca_fullfr.py @@ -457,7 +457,7 @@ def subtract_projection(cube, cube_ref, ncomp, scaling, mask_center_px, def pca_optimize_snr(cube, angle_list, (source_xy), fwhm, cube_ref=None, - mode='full', annulus_width=20, range_pcs=None, + mode='fullfr', annulus_width=20, range_pcs=None, svd_mode='lapack', scaling=None, mask_center_px=None, fmerit='px', min_snr=0, collapse='median', verbose=True, full_output=False, debug=False, plot=True): diff --git a/vip/phot/contrcurve.py b/vip/phot/contrcurve.py index 754fc587..b4bc6b73 100644 --- a/vip/phot/contrcurve.py +++ b/vip/phot/contrcurve.py @@ -161,7 +161,8 @@ def contrast_curve(cube, angle_list, psf_template, fwhm, pxscale, starphot, radvec_trans = transmission[1] f2 = InterpolatedUnivariateSpline(radvec_trans, trans, k=1) trans_interp = f2(rad_samp) - + thruput_interp *= trans_interp + # smoothing the noise vector using a Savitzky-Golay filter win = int(noise_samp.shape[0]*0.1) if win%2==0.: win += 1 @@ -215,9 +216,7 @@ def contrast_curve(cube, angle_list, psf_template, fwhm, pxscale, starphot, cont_curve_samp_corr = ((sigma_corr * noise_samp_sm)/thruput_interp) cont_curve_samp_corr[np.where(cont_curve_samp_corr<0)] = 1 cont_curve_samp_corr[np.where(cont_curve_samp_corr>1)] = 1 - - if transmission is not None: cont_curve_samp = cont_curve_samp*trans_interp - + # plotting if plot or debug: if student: From 14fac52c0069a5c989ffc1d1515c9ee85d185119 Mon Sep 17 00:00:00 2001 From: Carlos Gomez Date: Fri, 9 Dec 2016 15:44:35 +0100 Subject: [PATCH 3/3] no message --- vip/var/shapes.py | 67 ++++++++++++++++++++++++++++++++++++++++++++--- 1 file changed, 64 insertions(+), 3 deletions(-) diff --git a/vip/var/shapes.py b/vip/var/shapes.py index 5d160851..863da91c 100644 --- a/vip/var/shapes.py +++ b/vip/var/shapes.py @@ -15,12 +15,73 @@ 'get_annulus', 'get_annulus_quad', 'get_annulus_cube', - 'mask_circle'] + 'mask_circle', + 'create_ringed_spider_mask'] + +import numpy as np +from skimage.draw import polygon, circle + + +def create_ringed_spider_mask(im_shape, ann_out, ann_in=0, sp_width=10, sp_angle=0): + """ + Mask out information is outside the annulus and inside the spiders (zeros). + + Parameters + ---------- + im_shape : tuple of int + Tuple of length two with 2d array shape (Y,X). + ann_out : int + Outer radius of the annulus. + ann_in : int + Inner radius of the annulus. + sp_width : int + Width of the spider arms (3 branches). + sp_angle : int + angle of the first spider arm (on the positive horizontal axis) in + counter-clockwise sense. + + Returns + ------- + mask : array_like + 2d array of zeros and ones. + + """ + mask = np.zeros(im_shape) + + s = im_shape[0] + r = s/2. + theta = np.arctan2(sp_width/2., r) + + t0 = np.array([theta,np.pi-theta,np.pi+theta,np.pi*2.-theta]) + t1 = t0 + sp_angle/180. * np.pi + t2 = t1 + np.pi/3. + t3 = t2 + np.pi/3. + + x1 = r * np.cos(t1) + s/2. + y1 = r * np.sin(t1) + s/2. + x2 = r * np.cos(t2) + s/2. + y2 = r * np.sin(t2) + s/2. + x3 = r * np.cos(t3) + s/2. + y3 = r * np.sin(t3) + s/2. + + rr1, cc1 = polygon(y1, x1) + rr2, cc2 = polygon(y2, x2) + rr3, cc3 = polygon(y3, x3) + + cy, cx = frame_center(mask) + rr0, cc0 = circle(cy, cx, min(ann_out, cy)) + rr4, cc4 = circle(cy, cx, ann_in) + + mask[rr0,cc0] = 1 + mask[rr1,cc1] = 0 + mask[rr2,cc2] = 0 + mask[rr3,cc3] = 0 + mask[rr4,cc4] = 0 + return mask -import numpy as np def dist(yc,xc,y1,x1): - """ Returns the distance between two points. + """ Returns the Euclidean distance between two points. """ return np.sqrt((yc-y1)**2+(xc-x1)**2)