Merge pull request #78 from VChristiaens/master

Correction of minor bugs.

vip/negfc/mcmc_sampling.py

@@ -506,8 +506,7 @@ def mcmc_negfc_sampling(cubes, angs, psfn, ncomp, plsc, initial_state,
                                             cube_ref, svd_mode, scaling, fmerit,
                                             collapse]),
                                     threads=nproc)
-
-    duration_start = datetime.datetime.now()
+
     start = datetime.datetime.now()
 
     # #########################################################################

vip/preproc/cosmetics_ifs.py

@@ -153,7 +153,7 @@ def approx_stellar_position(cube, fwhm, return_test=False, verbose=False):
     #1/ Write a 2-columns array with indices of all max pixel values in the cube
     star_tmp_idx = np.zeros([n_z,2])
     star_approx_idx = np.zeros([n_z,2])
-    test_result = np.zeros(n_z)
+    test_result = np.ones(n_z)
     for zz in range(n_z):
         star_tmp_idx[zz] = peak_coordinates(obj_tmp[zz], fwhm[zz])
 
@@ -174,22 +174,22 @@ def approx_stellar_position(cube, fwhm, return_test=False, verbose=False):
     for zz in range(n_z):
         if ((star_tmp_idx[zz,0]<lim_inf_y) or (star_tmp_idx[zz,0]>lim_sup_y) or
             (star_tmp_idx[zz,1]<lim_inf_x) or (star_tmp_idx[zz,1]>lim_sup_x)):
-            test_result[zz] = 1
+            test_result[zz] = 0
 
     #3/ Replace by the median of neighbouring good coordinates if need be
     for zz in range(n_z):             
-        if test_result[zz] == 1:
+        if test_result[zz] == 0:
             ii= 1
             inf_neigh = max(0,zz-ii)
             sup_neigh = min(n_z-1,zz+ii)
-            while test_result[inf_neigh] == 1 and test_result[sup_neigh] == 1:
+            while test_result[inf_neigh] == 0 and test_result[sup_neigh] == 0:
                 ii=ii+1
                 inf_neigh = max(0,zz-ii)
                 sup_neigh = min(n_z-1,zz+ii)
-            if test_result[inf_neigh] == 0 and test_result[sup_neigh] == 0:
+            if test_result[inf_neigh] == 1 and test_result[sup_neigh] == 1:
                 star_approx_idx[zz] = np.floor((star_tmp_idx[sup_neigh]+ \
                                                 star_tmp_idx[inf_neigh])/2.)
-            elif test_result[inf_neigh] == 0: 
+            elif test_result[inf_neigh] == 1: 
                 star_approx_idx[zz] = star_tmp_idx[inf_neigh]
             else: star_approx_idx[zz] = star_tmp_idx[sup_neigh]
         else: star_approx_idx[zz] = star_tmp_idx[zz]

vip/preproc/recentering.py

@@ -753,7 +753,7 @@ def cube_recenter_gauss2d_fit(array, xy, fwhm=4, subi_size=5, nproc=1,
         Whether to save the shifts to a file in disk.
     offset : tuple of floats, optional
         If None the region of the frames used for the 2d Gaussian fit is shifted
-        to the center of the images (2d arrays). If a tuple is given is serves
+        to the center of the images (2d arrays). If a tuple is given it serves
         as the offset of the fitted area wrt the center of the 2d arrays.
     negative : {False, True}, optional
         If True a negative 2d Gaussian fit is performed.
@@ -773,9 +773,13 @@ def cube_recenter_gauss2d_fit(array, xy, fwhm=4, subi_size=5, nproc=1,
     if not array.ndim == 3:
         raise TypeError('Input array is not a cube or 3d array')
 
-    if isinstance(fwhm,int):  pass
-    elif isinstance(fwhm,float):  fwhm = int(np.round(fwhm))
-    else:  raise TypeError('fwhm is not a float or int')
+    n_frames = array.shape[0]
+    if isinstance(fwhm,int) or isinstance(fwhm,float):  
+        fwhm_tmp = fwhm 
+        fwhm = np.zeros(n_frames)
+        fwhm[:] = fwhm_tmp
+    subfr_sz = subi_size*fwhm
+    subfr_sz = subfr_sz.astype(int)
 
     if debug and array.shape[0]>20:
         msg = 'Debug with a big array will produce a very long output. '
@@ -795,7 +799,6 @@ def cube_recenter_gauss2d_fit(array, xy, fwhm=4, subi_size=5, nproc=1,
 
     if verbose:  start_time = time_ini()
 
-    n_frames = array.shape[0]
     cy, cx = frame_center(array[0])
     array_recentered = np.empty_like(array)  
 
@@ -806,18 +809,21 @@ def cube_recenter_gauss2d_fit(array, xy, fwhm=4, subi_size=5, nproc=1,
         bar = pyprind.ProgBar(n_frames, stream=1, 
                               title='2d Gauss-fitting, looping through frames')
         for i in range(n_frames):
-            res.append(_centroid_2dg_frame(array, i, subi_size*fwhm, 
-                                           pos_y, pos_x, negative, debug))
+            res.append(_centroid_2dg_frame(array, i, subfr_sz[i], 
+                                           pos_y, pos_x, negative, debug, fwhm[i]))
             bar.update()
         res = np.array(res)
     elif nproc>1:
         pool = Pool(processes=int(nproc))  
         res = pool.map(EFT, itt.izip(itt.repeat(_centroid_2dg_frame),
                                      itt.repeat(array),
                                      range(n_frames),
-                                     subi_size*fwhm, pos_y, pos_x,
+                                     subfr_sz, 
+                                     itt.repeat(pos_y), 
+                                     itt.repeat(pos_x),
                                      itt.repeat(negative),
-                                     itt.repeat(debug))) 
+                                     itt.repeat(debug),
+                                     fwhm)) 
         res = np.array(res)
         pool.close()
     y = cy - res[:,0]
@@ -883,7 +889,7 @@ def cube_recenter_moffat2d_fit(array, pos_y, pos_x, fwhm=4, subi_size=5,
         Whether to print to stdout the shifts or not. 
     unmoving_star : {False, True}, bool optional
         Whether the star centroid is expected to not move a lot within the 
-        frames of the input cube. If False, then an additional test is done to 
+        frames of the input cube. If True, then an additional test is done to 
         be sure the centroid fit returns a reasonable index value (close to the 
         median of the centroid indices in the other frames) - hence not taking 
         noise or a clump of uncorrected bad pixels.
@@ -954,7 +960,7 @@ def cube_recenter_moffat2d_fit(array, pos_y, pos_x, fwhm=4, subi_size=5,
         for i in range(n_frames):
             res.append(_centroid_2dm_frame(array, i, size[i], pos_y[i], 
                                            pos_x[i], star_approx_coords[i], 
-                                           star_not_present[i], negative))
+                                           star_not_present[i], negative, fwhm[i]))
             bar.update()
         res = np.array(res)
     elif nproc>1:
@@ -963,7 +969,7 @@ def cube_recenter_moffat2d_fit(array, pos_y, pos_x, fwhm=4, subi_size=5,
                                     itt.repeat(array), range(n_frames),
                                     size.tolist(), pos_y.tolist(), 
                                     pos_x.tolist(), star_approx_coords,
-                                    star_not_present, negative))
+                                    star_not_present, negative, fwhm))
         res = np.array(res)
         pool.close()
     y = cy - res[:,0]
@@ -985,28 +991,29 @@ def cube_recenter_moffat2d_fit(array, pos_y, pos_x, fwhm=4, subi_size=5,
         return array_recentered
 
 
-def _centroid_2dg_frame(cube, frnum, size, pos_y, pos_x, negative, debug=False):
+def _centroid_2dg_frame(cube, frnum, size, pos_y, pos_x, negative, debug=False, 
+                        fwhm=4):
     """ Finds the centroid by using a 2d gaussian fitting in one frame from a 
-    cube. To be called from whitin cube_recenter_gauss2d_fit().
+    cube. To be called from within cube_recenter_gauss2d_fit().
     """
     sub_image, y1, x1 = get_square_robust(cube[frnum], size=size, y=pos_y, 
                                           x=pos_x, position=True)
 
     # negative gaussian fit
     if negative:  sub_image = -sub_image + np.abs(np.min(-sub_image))
 
-    y_i, x_i = fit_2dgaussian(sub_image, crop=False, threshold=False, 
-                              sigfactor=1, debug=debug)
+    y_i, x_i = fit_2dgaussian(sub_image, crop=False, fwhmx=fwhm, fwhmy=fwhm, 
+                              threshold=False, sigfactor=1, debug=debug)
     y_i = y1 + y_i
     x_i = x1 + x_i
     return y_i, x_i
 
 
 def _centroid_2dm_frame(cube, frnum, size, pos_y, pos_x, 
                         star_approx_coords=None, star_not_present=None,
-                        negative=False):
+                        negative=False, fwhm=4):
     """ Finds the centroid by using a 2d moffat fitting in one frame from a 
-    cube. To be called from whitin cube_recenter_moffat2d_fit().
+    cube. To be called from within cube_recenter_moffat2d_fit().
     """
     sub_image, y1, x1 = get_square_robust(cube[frnum], size=size+1, y=pos_y, 
                                           x=pos_x,position=True)
@@ -1018,9 +1025,9 @@ def _centroid_2dm_frame(cube, frnum, size, pos_y, pos_x,
         if star_not_present:
             y_i,x_i = star_approx_coords
         else:
-            y_i, x_i = fit_2dmoffat(sub_image, y1, x1, full_output=False)
+            y_i, x_i = fit_2dmoffat(sub_image, y1, x1, full_output=False, fwhm=fwhm)
     else:
-        y_i, x_i = fit_2dmoffat(sub_image, y1, x1, full_output=False)
+        y_i, x_i = fit_2dmoffat(sub_image, y1, x1, full_output=False, fwhm=fwhm)
     return y_i, x_i
 
 

vip/preproc/rescaling.py

@@ -334,7 +334,7 @@ def check_scal_vector(scal_vec):
     correct = False
 
     if isinstance(scal_vec, list):
-        scal_list = scal_vec.copy()
+        scal_list = scal_vec[:]
         nz = len(scal_list)
         scal_vec = np.zeros(nz)
         for ii in range(nz):

vip/var/fit_2d.py

@@ -145,7 +145,7 @@ def fit_2dgaussian(array, crop=False, cent=None, cropsize=15, fwhmx=4, fwhmy=4,
 
 
 
-def fit_2dmoffat(array, yy, xx, full_output=False):
+def fit_2dmoffat(array, yy, xx, full_output=False,fwhm=4):
     """Fits a star/planet with a 2D circular Moffat PSF.
     
     Parameters
@@ -158,6 +158,11 @@ def fit_2dmoffat(array, yy, xx, full_output=False):
     xx : int
         X integer position of the first pixel (0,0) of the subimage in the 
         whole image.
+    full_output: bool, opt
+        Whether to return floor, height, mean_y, mean_x, fwhm, beta, or just 
+        mean_y, mean_x
+    fwhm: float, opt
+        First estimate of the fwhm
     
     Returns
     -------

vip/var/shapes.py

@@ -137,7 +137,7 @@ def get_square(array, size, y, x, position=False):
         raise TypeError('Input array is not a frame or 2d array.')
 
     if size%2!=0:  size -= 1 # making it odd to get the wing
-    wing = size/2
+    wing = int(size/2)
     # wing is added to the sides of the subframe center. Note the +1 when 
     # closing the interval (python doesn't include the endpoint)
     array_view = array[int(y-wing):int(y+wing+1),