Merge pull request #65 from carlgogo/master

VIP v0.6.0

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.          ")

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 *
 
+

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

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,17 +212,17 @@ 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"    
 
 
 def append_extension(filename, array):
     """Appends an extension to fits file. 
     """
-    fits.append(filename, array)
+    ap_fits.append(filename, array)
     print "\nFits extension appended"
 
 

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)
@@ -456,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):

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:  

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)