Merge pull request #315 from tomasstolker/optimize

optimize parameter FalsePositiveModule

docs/tutorial.rst

@@ -23,7 +23,7 @@ The pipeline works with two different components:
 2. The actual pipeline :class:`pynpoint.core.pypeline` which capsules a list of pipeline modules.
 
 .. important::
-   Pixel coordinates are zero-indexed, meaning that (x, y) = (0, 0) corresponds to the center of the pixel in the bottom-left corner of the image. The coordinate of the bottom-left corner is therefore (x, y) = (-0.5, -0.5). This means that there is an offset of 1.0 in both directions with respect to the pixel coordinates of DS9, for which the bottom-left corner is (x, y) = (0.5, 0.5).
+   Pixel coordinates are zero-indexed, meaning that (x, y) = (0, 0) corresponds to the center of the pixel in the bottom-left corner of the image. The coordinate of the bottom-left corner is therefore (x, y) = (-0.5, -0.5). This means that there is an offset of -1.0 in both directions with respect to the pixel coordinates of DS9, for which the bottom-left corner is (x, y) = (0.5, 0.5).
 
 .. _data-types:
 

pynpoint/processing/fluxposition.py

@@ -12,13 +12,12 @@
 import emcee
 
 from six.moves import range
-from scipy.stats import t
 from scipy.optimize import minimize
 from photutils import aperture_photometry, CircularAperture
 
 from pynpoint.core.processing import ProcessingModule
 from pynpoint.util.analysis import fake_planet, merit_function, false_alarm
-from pynpoint.util.image import create_mask, polar_to_cartesian
+from pynpoint.util.image import create_mask, polar_to_cartesian, get_image, cartesian_to_polar
 from pynpoint.util.mcmc import lnprob
 from pynpoint.util.module import progress, memory_frames, image_size_port, number_images_port, \
                                  rotate_coordinates
@@ -417,7 +416,7 @@ def _objective(arg):
                  x0=[pos_init[0], pos_init[1], self.m_magnitude],
                  method="Nelder-Mead",
                  tol=None,
-                 options={'xatol': self.m_tolerance, 'fatol': float("inf")})
+                 options={'xatol':self.m_tolerance, 'fatol':float("inf")})
 
         sys.stdout.write(" [DONE]\n")
         sys.stdout.flush()
@@ -434,7 +433,8 @@ def _objective(arg):
 class FalsePositiveModule(ProcessingModule):
     """
     Module to calculate the signal-to-noise ratio (SNR) and false positive fraction (FPF) at a
-    specified location in an image by using the Student's t-test (Mawet et al. 2014).
+    specified location in an image by using the Student's t-test (Mawet et al. 2014). Optionally,
+    the SNR can be optimized with the aperture position as free parameter.
     """
 
     def __init__(self,
@@ -443,12 +443,14 @@ def __init__(self,
                  ignore=False,
                  name_in="snr",
                  image_in_tag="im_arr",
-                 snr_out_tag="snr_fpf"):
+                 snr_out_tag="snr_fpf",
+                 optimize=False,
+                 **kwargs):
         """
         Constructor of FalsePositiveModule.
 
         :param position: The x and y position (pix) where the SNR and FPF is calculated. Note that
-                         the bottom left of the image is defined as (0, 0) so there is a -0.5
+                         the bottom left of the image is defined as (-0.5, -0.5) so there is a -1.0
                          offset with respect to the DS9 coordinate system. Aperture photometry
                          corrects for the partial inclusion of pixels at the boundary.
         :type position: (float, float)
@@ -467,10 +469,24 @@ def __init__(self,
                             counterclockwise direction with respect to the upward direction (i.e.,
                             East of North).
         :type snr_out_tag: str
+        :param optimize: Optimize the SNR. The aperture position is written in the history. The
+                         size of the aperture is kept fixed.
+        :type optimize: bool
+        :param kwargs:
+            See below.
+
+        :Keyword arguments:
+            **tolerance** (*float*) -- The absolute tolerance (pix) on the position for the
+                                       optimization to end. Default is set to 0.01 pix.
 
         :return: None
         """
 
+        if "tolerance" in kwargs:
+            self.m_tolerance = kwargs["tolerance"]
+        else:
+            self.m_tolerance = 1e-2
+
         super(FalsePositiveModule, self).__init__(name_in)
 
         self.m_image_in_port = self.add_input_port(image_in_tag)
@@ -479,62 +495,64 @@ def __init__(self,
         self.m_position = position
         self.m_aperture = aperture
         self.m_ignore = ignore
+        self.m_optimize = optimize
 
     def run(self):
         """
         Run method of the module. Calculates the SNR and FPF for a specified position in a post-
-        processed image with the Student's t-test (Mawet et al. 2014). This approach accounts
-        for small sample statistics.
+        processed image with the Student's t-test (Mawet et al. 2014). This approach assumes
+        Gaussian noise but accounts for small sample statistics.
 
         :return: None
         """
 
+        def _fpf_minimize(arg):
+            pos_x, pos_y = arg
+
+            _, _, fpf = false_alarm(image=image,
+                                    x_pos=pos_x,
+                                    y_pos=pos_y,
+                                    size=self.m_aperture,
+                                    ignore=self.m_ignore)
+
+            return fpf
+
         self.m_snr_out_port.del_all_data()
         self.m_snr_out_port.del_all_attributes()
 
         pixscale = self.m_image_in_port.get_attribute("PIXSCALE")
         self.m_aperture /= pixscale
 
         nimages = number_images_port(self.m_image_in_port)
-        center = self.m_image_in_port.get_shape()[-1]/2.
-
-        sep = math.sqrt((center-self.m_position[0])**2.+(center-self.m_position[1])**2.)
-        ang = (math.atan2(self.m_position[1]-center,
-                          self.m_position[0]-center)*180./math.pi - 90.)%360.
-
-        num_ap = int(math.pi*sep/self.m_aperture)
-        ap_theta = np.linspace(0, 2.*math.pi, num_ap, endpoint=False)
-
-        if self.m_ignore:
-            num_ap -= 2
-            ap_theta = np.delete(ap_theta, [1, np.size(ap_theta)-1])
 
         for j in range(nimages):
             progress(j, nimages, "Running FalsePositiveModule...")
 
-            if nimages == 1:
-                image = self.m_image_in_port.get_all()
-                if image.ndim == 3:
-                    image = np.squeeze(image, axis=0)
+            image = get_image(self.m_image_in_port, j, nimages)
 
-            else:
-                image = self.m_image_in_port[j, ]
+            if self.m_optimize:
+                result = minimize(fun=_fpf_minimize,
+                                  x0=[self.m_position[0], self.m_position[1]],
+                                  method="Nelder-Mead",
+                                  tol=None,
+                                  options={'xatol':self.m_tolerance, 'fatol':float("inf")})
 
-            ap_phot = np.zeros(num_ap)
-            for i, theta in enumerate(ap_theta):
-                x_tmp = center + (self.m_position[0]-center)*math.cos(theta) - \
-                                 (self.m_position[1]-center)*math.sin(theta)
-                y_tmp = center + (self.m_position[0]-center)*math.sin(theta) + \
-                                 (self.m_position[1]-center)*math.cos(theta)
+                _, snr, fpf = false_alarm(image=image,
+                                          x_pos=result.x[0],
+                                          y_pos=result.x[1],
+                                          size=self.m_aperture,
+                                          ignore=self.m_ignore)
 
-                aperture = CircularAperture((x_tmp, y_tmp), self.m_aperture)
-                phot_table = aperture_photometry(image, aperture, method='exact')
-                ap_phot[i] = phot_table['aperture_sum']
+                sep, ang = cartesian_to_polar(image, result.x[0], result.x[1])
 
-            snr = (ap_phot[0] - np.mean(ap_phot[1:])) / \
-                  (np.std(ap_phot[1:]) * math.sqrt(1.+1./float(num_ap-1)))
+            else:
+                _, snr, fpf = false_alarm(image=image,
+                                          x_pos=self.m_position[0],
+                                          y_pos=self.m_position[1],
+                                          size=self.m_aperture,
+                                          ignore=self.m_ignore)
 
-            fpf = 1. - t.cdf(snr, num_ap-2)
+                sep, ang = cartesian_to_polar(image, self.m_position[0], self.m_position[1])
 
             result = np.column_stack((self.m_position[0],
                                       self.m_position[1],

pynpoint/util/analysis.py

@@ -14,7 +14,7 @@
 from photutils import aperture_photometry, CircularAperture, EllipticalAperture
 from six.moves import range
 
-from pynpoint.util.image import shift_image
+from pynpoint.util.image import shift_image, image_center_subpixel
 
 
 def false_alarm(image,
@@ -23,14 +23,16 @@ def false_alarm(image,
                 size,
                 ignore):
     """
-    Function for the formal t-test for high-contrast imaging at small working angles, as well as
-    the related false positive fraction (Mawet et al. 2014).
+    Function for the formal t-test for high-contrast imaging at small working angles and the
+    related false positive fraction (Mawet et al. 2014).
 
     :param image: Input image.
     :type image: numpy.ndarray
-    :param x_pos: Position (pix) along the x-axis.
+    :param x_pos: Position (pix) along the horizontal axis. The pixel coordinates of the
+                  bottom-left corner of the image are (-0.5, -0.5).
     :type x_pos: float
-    :param y_pos: Position (pix) along the y-axis.
+    :param y_pos: Position (pix) along the vertical axis. The pixel coordinates of the
+                  bottom-left corner of the image are (-0.5, -0.5).
     :type y_pos: float
     :param size: Aperture radius (pix).
     :type size: float
@@ -41,7 +43,7 @@ def false_alarm(image,
     :rtype: float, float, float
     """
 
-    center = (np.size(image, 0)/2., np.size(image, 1)/2.)
+    center = image_center_subpixel(image)
     radius = math.sqrt((center[0]-y_pos)**2.+(center[1]-x_pos)**2.)
 
     num_ap = int(math.pi*radius/size)
@@ -53,13 +55,13 @@ def false_alarm(image,
 
     if num_ap < 3:
         raise ValueError("Number of apertures (num_ap=%s) is too small to calculate the "
-                         "false positive fraction. Increase the lower limit of the "
-                         "separation argument." % num_ap)
+                         "false positive fraction." % num_ap)
 
     ap_phot = np.zeros(num_ap)
     for i, theta in enumerate(ap_theta):
         x_tmp = center[1] + (x_pos-center[1])*math.cos(theta) - \
                             (y_pos-center[0])*math.sin(theta)
+
         y_tmp = center[0] + (x_pos-center[1])*math.sin(theta) + \
                             (y_pos-center[0])*math.cos(theta)
 

pynpoint/util/image.py

@@ -18,7 +18,7 @@ def image_center_pixel(image):
     """
     Function to get the pixel position of the image center. Note that this position
     can not be unambiguously defined for an even-sized image. Python indexing starts
-    at 0 so the coordinates of the pixel in the bottom left corner are (0, 0).
+    at 0 so the coordinates of the pixel in the bottom-left corner are (0, 0).
 
     :param image: Input image (2D or 3D).
     :type image: numpy.ndarray
@@ -88,11 +88,11 @@ def crop_image(image,
     if size%2 == 0:
         size += 1
 
-    x_start = center[1] - (size-1) // 2
-    x_end = center[1] + (size-1) // 2 + 1
+    x_start = center[1] - (size-1)//2
+    x_end = center[1] + (size-1)//2 + 1
 
-    y_start = center[0] - (size-1) // 2
-    y_end = center[0] + (size-1) // 2 + 1
+    y_start = center[0] - (size-1)//2
+    y_end = center[0] + (size-1)//2 + 1
 
     if x_start < 0 or y_start < 0 or x_end > image.shape[1] or y_end > image.shape[0]:
         raise ValueError("Target image resolution does not fit inside the input image resolution.")
@@ -228,6 +228,32 @@ def scale_image(image,
 
     return im_scale * (sum_before / sum_after)
 
+def cartesian_to_polar(image, x_pos, y_pos):
+    """
+    Function to convert pixel coordinates to polar coordinates.
+
+    :param image: Input image (2D or 3D).
+    :type image: numpy.ndarray
+    :param x_pos: Pixel coordinate along the horizontal axis. The bottom left corner of the image
+                  is (-0.5, -0.5).
+    :type x_pos: float
+    :param y_pos: Pixel coordinate along the vertical axis. The bottom left corner of the image
+                  is (-0.5, -0.5).
+    :type y_pos: float
+
+    :return: Polar coordinates as separation (pix) and position angle (deg). The angle is measured
+             counterclockwise with respect to the positive y-axis.
+    :rtype: float, float
+    """
+
+    center = image_center_subpixel(image) # (y, x)
+
+    sep = math.sqrt((center[1]-x_pos)**2.+(center[0]-y_pos)**2.)
+    ang = math.atan2(y_pos-center[1], x_pos-center[0])
+    ang = (math.degrees(ang)-90.)%360.
+
+    return sep, ang
+
 def polar_to_cartesian(image, sep, ang):
     """
     Function to convert polar coordinates to pixel coordinates.
@@ -250,3 +276,29 @@ def polar_to_cartesian(image, sep, ang):
     y_pos = center[0] + sep*math.sin(math.radians(ang+90.))
 
     return x_pos, y_pos
+
+def get_image(input_port, index, nimages):
+    """
+    Function to get an image from an input port based on its index.
+
+    :param input_port: Input port with the stack of images.
+    :type input_port: pynpoint.core.dataio.InputPort
+    :param index: Index in the stack of images.
+    :type index: int
+    :param nimages: Total number of images in the stack.
+    :type nimages: int
+
+    :return: Selected image.
+    :rtype: numpy.ndarray
+    """
+
+    if nimages == 1:
+        image = input_port.get_all()
+
+        if image.ndim == 3:
+            image = np.squeeze(image, axis=0)
+
+    else:
+        image = input_port[index, ]
+
+    return image

tests/test_processing/test_fluxposition.py

@@ -176,10 +176,12 @@ def test_false_positive(self):
         self.pipeline.run_module("false")
 
         data = self.pipeline.get_data("snr_fpf")
-        assert np.allclose(data[0, 2], 0.5280553948214145, rtol=limit, atol=0.)
-        assert np.allclose(data[0, 3], 94.39870535499551, rtol=limit, atol=0.)
-        assert np.allclose(data[0, 4], 8.542166952478182, rtol=limit, atol=0.)
-        assert np.allclose(data[0, 5], 9.54772666372783e-07, rtol=limit, atol=0.)
+        assert np.allclose(data[0, 0], 31.0, rtol=limit, atol=0.)
+        assert np.allclose(data[0, 1], 49.0, rtol=limit, atol=0.)
+        assert np.allclose(data[0, 2], 0.513710034941892, rtol=limit, atol=0.)
+        assert np.allclose(data[0, 3], 93.01278750418334, rtol=limit, atol=0.)
+        assert np.allclose(data[0, 4], 7.633199090133858, rtol=limit, atol=0.)
+        assert np.allclose(data[0, 5], 3.029521252528866e-06, rtol=limit, atol=0.)
 
     def test_simplex_minimization(self):
 

tests/test_processing/test_limits.py

@@ -89,7 +89,7 @@ def test_contrast_curve(self):
 
             data = self.pipeline.get_data("limits_"+item)
             assert np.allclose(data[0, 0], 5.00000000e-01, rtol=limit, atol=0.)
-            assert np.allclose(data[0, 1], 6.647211555936695, rtol=limit, atol=0.)
-            assert np.allclose(data[0, 2], 0.07837070875760642, rtol=limit, atol=0.)
+            assert np.allclose(data[0, 1], 6.791207747570333, rtol=limit, atol=0.)
+            assert np.allclose(data[0, 2], 0.09527007675032029, rtol=limit, atol=0.)
             assert np.allclose(data[0, 3], 0.0002012649090622487, rtol=limit, atol=0.)
             assert data.shape == (1, 4)