Some small bugs and suggestions

[VChristiaens]

1/ in frame_center_radon() (calib/recentering.py):

• cost_bound variable is only defined when plot=True, but is then later called even if plot=False.
• (added in the edit) When I want to use the wavelet filter, I get systematically an error of this type: "ValueError: Level value of 100 is too high. Maximum allowed is 6." I see that the value of 100 is actually hardcoded at l.432. When I change the value to say 5, I get an other error stemming from wavelet_denoise(): "AttributeError: 'module' object has no attribute 'thresholding'"

2/ in contrast_curve() (phot/contr_curve.py):

• l.173-175: when noise_samp size is < 20 (i.e. when your cropped frames can carry less than 20 concentric annuli), an error is raised by the savgol_filter function because the polynomial order is set to 2 and is less or equal than the filter window length (set to 0.1*noise_samp size). Maybe we can allow a bit more flexibility? A contrast curve could still make sense with less than 20 concentric annuli, there is maybe just no need to use a savgol_filter in that kind of case.
• when the contrast curve is saved, the pdf appears truncated => simple solution: figsize set to (10,6) instead of (8,4)

3/ in _centroid_2dm_frame() (calib/recentering.py): I'd remove the condition "if np.allclose(miny, cy, atol=2) and np.allclose(minx, cx, atol=2):" which I find arbitrary to detect automatically if the data were taken with the AGPM or not! With some recent NACO+AGPM data, half of the frames pass the condition, the other half not. It is probably better to implement it in the same way as _centroid_2dg_frame()

4/ in var/shapes.py: I implemented the equivalent of get_circle and get_annulus, but with ellipses. I found it useful in the case of inclined disks. I let you decide whether to include them in vip or not:

def get_ellipse(array, a, b, PA, output_values=False, cy=None, cx=None, output_indices=False):

"""
Returns a centered elliptical region from a 2d ndarray. All the rest 
pixels are set to zeros.

Parameters
----------
array : array_like
    Input 2d array or image. 
a : flt
    Semi-major axis.
b : flt
    Semi-minor axis.
PA : deg
    The PA of the semi-major axis.
output_values : {False, True}, optional
    If True returns the values of the pixels in the annulus.
cy, cx : int
    Coordinates of the circle center.
output_indices : {False, True}, optional
    If True returns the indices inside the annulus.

Returns
-------
Depending on output_values, output_indices:
values : array_like
    1d array with the values of the pixels in the circular region.
array_masked : array_like
    Input array with the circular mask applied.
y, x : array_like
    Coordinates of pixels in circle.
"""

if not array.ndim == 2:
    raise TypeError('Input array is not a frame or 2d array.')
sy, sx = array.shape
if not cy and not cx:
    cy, cx = frame_center(array, verbose=False)

# Definition of other parameters of the ellipse
f = np.sqrt(a**2-b**2) #distance between center and foci of the ellipse
PA_rad = np.deg2rad(PA)    
pos_f1 = (cy+f*np.cos(PA_rad),cx+f*np.sin(PA_rad)) #coords of first focus
pos_f2 = (cy-f*np.cos(PA_rad),cx-f*np.sin(PA_rad)) # coords of second focus

yy, xx = np.ogrid[:sy, :sx] 
# ogrid is a multidim mesh creator (faster than mgrid)
ellipse = dist(yy,xx,pos_f1[0],pos_f1[1]) + dist(yy,xx,pos_f2[0],pos_f2[1])       
ellipse_mask = ellipse < 2*a   # mask of 1's and 0's

if output_values and not output_indices:
    values = array[ellipse_mask]
    return values
elif output_indices and not output_values:      
    indices = np.array(np.where(ellipse_mask))
    y = indices[0]
    x = indices[1]
    return y, x
elif output_indices and output_values:
    msg =  'output_values and output_indices cannot be both True.'
    raise ValueError(msg)
else:
    array_masked = array*ellipse_mask
    return array_masked

def get_ell_annulus(array, a, b, PA, width, output_values=False, output_indices=False, cy=None, cx=None):

"""Returns a centered elliptical annulus from a 2d ndarray. All the rest pixels are 
set to zeros. 

Parameters
----------
array : array_like
    Input 2d array or image. 
a : flt
    Semi-major axis.
b : flt
    Semi-minor axis.
PA : deg
    The PA of the semi-major axis.
width : flt
    The size of the annulus along the semi-major axis; it is proportionnally 
    thinner along the semi-minor axis).
output_values : {False, True}, optional
    If True returns the values of the pixels in the annulus.
output_indices : {False, True}, optional
    If True returns the indices inside the annulus.
cy,cx: float, optional
    Location of the center of the annulus to be defined. If not provided, 
it assumes the annuli are centered on the frame.

Returns
-------
Depending on output_values, output_indices:
values : array_like
    1d array with the values of the pixels in the annulus.
array_masked : array_like
    Input array with the annular mask applied.
y, x : array_like
    Coordinates of pixels in annulus.

"""
if not array.ndim == 2:
    raise TypeError('Input array is not a frame or 2d array.')
if cy is None or cx is None:
    cy, cx = frame_center(array)
sy, sx = array.shape

width_a = width
width_b = width*b/a

# Definition of big ellipse
f_big = np.sqrt((a+width_a/2.)**2-(b+width_b/2.)**2) #distance between center and foci of the ellipse
PA_rad = np.deg2rad(PA)    
pos_f1_big = (cy+f_big*np.cos(PA_rad),cx+f_big*np.sin(PA_rad)) #coords of first focus
pos_f2_big = (cy-f_big*np.cos(PA_rad),cx-f_big*np.sin(PA_rad)) # coords of second focus

# Definition of small ellipse
f_sma = np.sqrt((a-width_a/2.)**2-(b-width_b/2.)**2) #distance between center and foci of the ellipse   
pos_f1_sma = (cy+f_sma*np.cos(PA_rad),cx+f_sma*np.sin(PA_rad)) #coords of first focus
pos_f2_sma = (cy-f_sma*np.cos(PA_rad),cx-f_sma*np.sin(PA_rad)) # coords of second focus

yy, xx = np.ogrid[:sy, :sx] 
big_ellipse = dist(yy,xx,pos_f1_big[0],pos_f1_big[1]) + dist(yy,xx,pos_f2_big[0],pos_f2_big[1])
small_ellipse = dist(yy,xx,pos_f1_sma[0],pos_f1_sma[1]) + dist(yy,xx,pos_f2_sma[0],pos_f2_sma[1])         
ell_ann_mask = (big_ellipse < 2*(a+width/2.)) & (small_ellipse >= 2*(a-width/2.))  # mask of 1's and 0's

if output_values and not output_indices:
    values = array[ell_ann_mask]
    return values
elif output_indices and not output_values:      
    indices = np.array(np.where(ell_ann_mask))
    y = indices[0]
    x = indices[1]
    return y, x
elif output_indices and output_values:
    msg =  'output_values and output_indices cannot be both True.'
    raise ValueError(msg)
else:
    array_masked = array*ell_ann_mask
    return array_masked

[carlos-gg]

@VChristiaens. Took me a while to address your comments, but it makes sense to polish VIP as much as possible now that the paper is accepted.

1) cost_bound variable problem fixed. Wavelet filtering removed. Probably the API of pywavelets changed and I prefer to remove the functionality and package dependence.

2) The smoothie is now optional (boolean argument)

3) Done

4) Included the functions. They seem to be working and are consistent with the functionality of VIP and the shapes module.

Cheers, C.