Add refine function and expand unit tests

trackpy/artificial.py

@@ -3,6 +3,7 @@
 import six
 import numpy as np
 import pandas as pd
+import warnings
 from trackpy.find import drop_close
 from trackpy.utils import validate_tuple
 from trackpy.preprocessing import bandpass
@@ -474,50 +475,35 @@ def f(self, noise=0):
                 result[col] = float(s)
         return result
 
-def feat_brightfield(r, ndim, radii, dark_value, bright_value):
-    """ Brightfield particle at r = 0. """
+def feat_brightfield(r, ndim, radii, dark_value, bright_value, dip):
+    """ Brightfield particle with intensity dip in center at r = 0. """
     image = np.zeros_like(r)
 
-    # TODO: anisotropic sizes
+    if not np.isclose(radii, radii[0]).all():
+        warnings.warn("Anisotropic particle sizes for brightfield features" +
+                      " are not supported! I'll use the first value as radius")
+
     radius = radii[0]-0.5
 
-    thickness = 3
-    mask = r <= radius
+    thickness = radius*0.25
+
+    mask = r < radius
     mask_ring = mask & (r > (radius-thickness))
-    image[mask_ring] += dark_value
-    image[mask] += bright_value*np.exp(-(r[mask]/((radius-thickness/3.0)**2))**2)
 
-    return image
+    if dip:
+        mask_dip = r < thickness
+        image[mask_dip] += 1.5*dark_value
 
-def draw_feature_brightfield(image, position, size, **kwargs):
-    """ Draws a brightfield feature and adds it to the image at given
-    position. The given function will be evaluated at each pixel coordinate,
-    no averaging or convolution is done.
+    image[mask_ring] += dark_value
 
-    Parameters
-    ----------
-    image : ndarray
-        image to draw features on
-    position : iterable
-        coordinates of feature position
-    size : number
-        the size of the feature (meaning depends on feature, for feat_gauss,
-        it is the radius of gyration)
-    kwargs : keyword arguments are passed to draw_feature
+    gauss_radius = radius-1.6*thickness
+    image[mask] += bright_value*np.exp(-(r[mask]/(gauss_radius**2))**2)
 
-    See also
-    --------
-    draw_feature
-    """
-    kwargs['feat_func'] = feat_brightfield
-    kwargs['radii'] = np.array(size)/2.0
-    kwargs['dark_value'] = -0.3
-    kwargs['bright_value'] = 0.8
-    draw_feature(image, position, size, **kwargs)
+    return image
 
 
-def draw_features_brightfield(shape, positions, size, noise_level=0, 
-                              bitdepth=8, background=0.5, **kwargs):
+def draw_features_brightfield(shape, positions, size, noise_level=0,
+                              bitdepth=8, background=0.5, dip=False, **kwargs):
     """ Generates an image with features at given positions. A feature with
     position x will be centered around pixel x. In other words, the origin of
     the output image is located at the center of pixel (0, 0).
@@ -571,14 +557,15 @@ def draw_features_brightfield(shape, positions, size, noise_level=0,
     if np.max(image) > max_brightness:
         image = image/np.max(image)*max_brightness
 
+    kwargs['feat_func'] = feat_brightfield
+    kwargs['radii'] = np.array(size)/2.0
+    kwargs['dark_value'] = -0.3
+    kwargs['bright_value'] = 0.8
+    kwargs['dip'] = dip
+
     for pos in positions:
-        draw_feature_brightfield(image, pos, size, max_value=max_brightness,
-                                 **kwargs)
+        draw_feature(image, pos, size, max_value=max_brightness, **kwargs)
     result = image.clip(0, 2**bitdepth - 1).astype(dtype)
 
-    #import matplotlib.pyplot as plt
-    #plt.imshow(result, cmap='gray')
-    #plt.show(block=True)
-
     return result
 

trackpy/locate/brightfield_ring.py

@@ -7,12 +7,14 @@
 from pandas import (DataFrame,)
 
 from ..find import (grey_dilation, where_close)
+from ..refine import (refine_brightfield_ring,)
 from ..utils import (validate_tuple, default_pos_columns)
 from ..preprocessing import (bandpass, convert_to_int)
 
 
 def locate_brightfield_ring(raw_image, diameter, separation=None, noise_size=1,
-                            smoothing_size=None, threshold=None):
+                            smoothing_size=None, threshold=None,
+                            previous_coords=None):
     """Locate particles imaged in brightfield mode of some approximate size in
     an image.
 
@@ -44,6 +46,9 @@ def locate_brightfield_ring(raw_image, diameter, separation=None, noise_size=1,
         Clip bandpass result below this value. Thresholding is done on the
         already background-subtracted image.
         By default, 1 for integer images and 1/255 for float images.
+    previous_coords : DataFrame([x, y, size])
+        Optional previous particle positions from the preceding frame to use as
+        starting point for the refinement instead of the intensity peaks.
 
     Returns
     -------
@@ -120,33 +125,37 @@ def locate_brightfield_ring(raw_image, diameter, separation=None, noise_size=1,
 
     pos_columns = default_pos_columns(image.ndim)
 
-    # Find local maxima.
-    # Define zone of exclusion at edges of image, avoiding
-    #   - Features with incomplete image data ("radius")
-    #   - Extended particles that cannot be explored during subpixel
-    #       refinement ("separation")
-    #   - Invalid output of the bandpass step ("smoothing_size")
-    margin = tuple([max(rad, sep // 2 - 1, sm // 2) for (rad, sep, sm) in
-                    zip(radius, separation, smoothing_size)])
-    # Find features with minimum separation distance of `separation`. This
-    # excludes detection of small features close to large, bright features
-    # using the `maxsize` argument.
-    coords = grey_dilation(image, separation, margin=margin, precise=False)
-
-    coords_df = DataFrame(columns=pos_columns, data=coords)
+    if previous_coords is None or len(previous_coords) == 0:
+        # Find local maxima.
+        # Define zone of exclusion at edges of image, avoiding
+        #   - Features with incomplete image data ("radius")
+        #   - Extended particles that cannot be explored during subpixel
+        #       refinement ("separation")
+        #   - Invalid output of the bandpass step ("smoothing_size")
+        margin = tuple([max(rad, sep // 2 - 1, sm // 2) for (rad, sep, sm) in
+                        zip(radius, separation, smoothing_size)])
+        # Find features with minimum separation distance of `separation`. This
+        # excludes detection of small features close to large, bright features
+        # using the `maxsize` argument.
+        coords = grey_dilation(image, separation, margin=margin, precise=False)
+
+        coords_df = DataFrame(columns=pos_columns, data=coords)
+    else:
+        coords_df = previous_coords
 
     if len(coords_df) == 0:
         warnings.warn("No particles found in the image before refinement.")
         return coords_df
 
+    coords_df = refine_brightfield_ring(image, radius, coords_df,
+                                        pos_columns=pos_columns)
+
     # Flat peaks return multiple nearby maxima. Eliminate duplicates.
     if np.all(np.greater(separation, 0)):
         to_drop = where_close(coords_df[pos_columns], separation)
         coords_df.drop(to_drop, axis=0, inplace=True)
         coords_df.reset_index(drop=True, inplace=True)
 
-    # TODO: refine particle positions
-
     if len(coords_df) == 0:
         warnings.warn("No particles found in the image after refinement.")
         return coords_df

trackpy/refine/__init__.py

@@ -1,2 +1,3 @@
 from .center_of_mass import refine_com, refine_com_arr
 from .least_squares import refine_leastsq
+from .brightfield_ring import refine_brightfield_ring

trackpy/refine/brightfield_ring.py

@@ -0,0 +1,36 @@
+import numpy as np
+import pandas as pd
+
+from ..utils import (validate_tuple, guess_pos_columns, default_pos_columns)
+
+
+def refine_brightfield_ring(image, radius, coords_df, pos_columns=None):
+    """Find the center of mass of a brightfield feature starting from an
+    estimate.
+
+    Parameters
+    ----------
+    image : array (any dimension)
+        processed image, used for locating center of mass
+    coords_df : DataFrame
+        estimated positions
+    pos_columns: list of strings, optional
+        Column names that contain the position coordinates.
+        Defaults to ``['y', 'x']`` or ``['z', 'y', 'x']``, if ``'z'`` exists.
+    """
+    if pos_columns is None:
+        pos_columns = guess_pos_columns(coords_df)
+
+    radius = validate_tuple(radius, image.ndim)
+
+    if pos_columns is None:
+        pos_columns = default_pos_columns(image.ndim)
+
+    columns = pos_columns + ['size']
+
+    if len(coords_df) == 0:
+        return pd.DataFrame(columns=columns)
+
+    refined = coords_df
+
+    return refined

trackpy/tests/locate/test_brightfield_ring.py

@@ -4,10 +4,12 @@
 from six.moves import range
 import os
 import warnings
+import logging
 
 import nose
 import numpy as np
 from pandas import DataFrame
+import matplotlib.pyplot as plt
 from numpy.testing import assert_allclose, assert_array_less
 from pandas.util.testing import assert_produces_warning
 
@@ -18,11 +20,9 @@
 from trackpy.tests.common import sort_positions, StrictTestCase
 from trackpy.locate.brightfield_ring import locate_brightfield_ring
 
-
 path, _ = os.path.split(os.path.abspath(__file__))
 
-
-def artificial_image(shape, count, radius, noise_level, **kwargs):
+def artificial_image(shape, count, radius, noise_level, dip=False, **kwargs):
     radius = tp.utils.validate_tuple(radius, len(shape))
     # tp.locate ignores a margin of size radius, take 1 px more to be safe
     margin = tuple([r + 1 for r in radius])
@@ -31,28 +31,59 @@ def artificial_image(shape, count, radius, noise_level, **kwargs):
     separation = tuple([d * 3 for d in diameter])
     cols = ['x', 'y', 'z'][:len(shape)][::-1]
     pos = gen_nonoverlapping_locations(shape, count, separation, margin)
-    image = draw_features_brightfield(shape, pos, size, noise_level)
+    image = draw_features_brightfield(shape, pos, size, noise_level, dip=dip)
     f = locate_brightfield_ring(image, diameter, **kwargs)
     actual = pandas_sort(f[cols], cols)
     expected = pandas_sort(DataFrame(pos, columns=cols), cols)
+
+    if True:
+        plot_result(image, actual, expected)
+
     return actual, expected
 
+def plot_result(image, actual, expected):
+    plt.imshow(image, cmap='gray')
+    plt.scatter(actual['x'], actual['y'], c='r', label='actual', s=3)
+    plt.scatter(expected['x'], expected['y'], c='g', label='expected', s=3)
+    plt.legend()
+    plt.show(block=True)
+
+
 
 class TestLocateBrightfieldRing(StrictTestCase):
+    def setUp(self):
+        mpl_logger = logging.getLogger('matplotlib')
+        mpl_logger.setLevel(logging.WARNING)
+
     def test_multiple_simple_sparse(self):
         actual, expected = artificial_image((200, 300), 4, 5, noise_level=0)
         assert_allclose(actual, expected, atol=0.5)
 
+    def test_multiple_simple_sparse_dip(self):
+        actual, expected = artificial_image((200, 300), 4, 5, noise_level=0,
+                                            dip=True)
+        assert_allclose(actual, expected, atol=0.5)
+
     def test_multiple_noisy_sparse(self):
         #  4% noise
         actual, expected = artificial_image((200, 300), 4, 5, noise_level=10)
         assert_allclose(actual, expected, atol=0.5)
 
+    def test_multiple_noisy_sparse_dip(self):
+        actual, expected = artificial_image((200, 300), 4, 5, noise_level=10,
+                                            dip=True)
+        assert_allclose(actual, expected, atol=0.5)
+
     def test_multiple_more_noisy_sparse(self):
         # 20% noise
         actual, expected = artificial_image((200, 300), 4, 5, noise_level=51)
         assert_allclose(actual, expected, atol=0.5)
 
+    def test_multiple_more_noisy_sparse_dip(self):
+        actual, expected = artificial_image((200, 300), 4, 5, noise_level=51,
+                                            dip=True)
+        assert_allclose(actual, expected, atol=0.5)
+
 
 if __name__ == '__main__':
     import nose