Merge pull request #865 from jni/profile-line

Move `profile_line()` out of viewer and refactor

skimage/measure/__init__.py

@@ -4,6 +4,7 @@
 from ._structural_similarity import structural_similarity
 from ._polygon import approximate_polygon, subdivide_polygon
 from ._moments import moments, moments_central, moments_normalized, moments_hu
+from .profile import profile_line
 from .fit import LineModel, CircleModel, EllipseModel, ransac
 from .block import block_reduce
 
@@ -24,4 +25,5 @@
            'moments_normalized',
            'moments_hu',
            'marching_cubes',
-           'mesh_surface_area']
+           'mesh_surface_area',
+           'profile_line']

skimage/measure/profile.py

@@ -0,0 +1,87 @@
+import numpy as np
+import scipy.ndimage as nd
+
+
+def profile_line(img, src, dst, linewidth=1,
+                 order=1, mode='constant', cval=0.0):
+    """Return the intensity profile of an image measured along a scan line.
+
+    Parameters
+    ----------
+    img : numeric array, shape (M, N[, C])
+        The image, either grayscale (2D array) or multichannel 
+        (3D array, where the final axis contains the channel
+        information).
+    src : 2-tuple of numeric scalar (float or int)
+        The start point of the scan line.
+    dst : 2-tuple of numeric scalar (float or int)
+        The end point of the scan line.
+    linewidth : int, optional
+        Width of the scan, perpendicular to the line
+    order : int in {0, 1, 2, 3, 4, 5}, optional
+        The order of the spline interpolation to compute image values at
+        non-integer coordinates. 0 means nearest-neighbor interpolation.
+    mode : string, one of {'constant', 'nearest', 'reflect', 'wrap'}, optional
+        How to compute any values falling outside of the image.
+    cval : float, optional
+        If `mode` is 'constant', what constant value to use outside the image.
+
+    Returns
+    -------
+    return_value : array
+        The intensity profile along the scan line. The length of the profile
+        is the ceil of the computed length of the scan line.
+
+    Examples
+    --------
+    >>> x = np.array([[1, 1, 1, 2, 2, 2]])
+    >>> img = np.vstack([np.zeros_like(x), x, x, x, np.zeros_like(x)])
+    >>> img
+    array([[0, 0, 0, 0, 0, 0],
+           [1, 1, 1, 2, 2, 2],
+           [1, 1, 1, 2, 2, 2],
+           [1, 1, 1, 2, 2, 2],
+           [0, 0, 0, 0, 0, 0]])
+    >>> profile_line(img, (2, 1), (2, 4))
+    array([ 1.,  1.,  2.,  2.])
+
+    Notes
+    -----
+    The destination point is included in the profile, in contrast to
+    standard numpy indexing.
+    """
+    src_row, src_col = src = np.asarray(src, dtype=float)
+    dst_row, dst_col = dst = np.asarray(dst, dtype=float)
+    d_row, d_col = dst - src
+    theta = np.arctan2(d_row, d_col)
+
+    length = np.ceil(np.hypot(d_row, d_col) + 1)
+    # we add one above because we include the last point in the profile
+    # (in contrast to standard numpy indexing)
+    line_col = np.linspace(src_col, dst_col, length)
+    line_row = np.linspace(src_row, dst_row, length)
+
+    # we subtract 1 from linewidth to change from pixel-counting
+    # (make this line 3 pixels wide) to point distances (the
+    # distance between pixel centers)
+    col_width = (linewidth - 1) * np.sin(-theta) / 2
+    row_width = (linewidth - 1) * np.cos(theta) / 2
+    perp_rows = np.array([np.linspace(row_i - row_width, row_i + row_width,
+                                      linewidth) for row_i in line_row])
+    perp_cols = np.array([np.linspace(col_i - col_width, col_i + col_width,
+                                      linewidth) for col_i in line_col])
+    perp_lines = np.array([perp_rows, perp_cols])
+
+    if img.ndim == 3:
+        pixels = [nd.map_coordinates(img[..., i], perp_lines,
+                                     order=order, mode=mode, cval=cval)
+                  for i in range(img.shape[2])]
+        pixels = np.transpose(np.asarray(pixels), (1, 2, 0))
+    else:
+        pixels = nd.map_coordinates(img, perp_lines,
+                                    order=order, mode=mode, cval=cval)
+
+    intensities = pixels.mean(axis=1)
+
+    return intensities
+

skimage/measure/tests/test_profile.py

@@ -0,0 +1,110 @@
+from numpy.testing import assert_equal, assert_almost_equal
+import numpy as np
+
+from skimage.measure import profile_line
+
+image = np.arange(100).reshape((10, 10)).astype(np.float)
+
+def test_horizontal_rightward():
+    prof = profile_line(image, (0, 2), (0, 8), order=0)
+    expected_prof = np.arange(2, 9)
+    assert_equal(prof, expected_prof)
+
+
+def test_horizontal_leftward():
+    prof = profile_line(image, (0, 8), (0, 2), order=0)
+    expected_prof = np.arange(8, 1, -1)
+    assert_equal(prof, expected_prof)
+
+
+def test_vertical_downward():
+    prof = profile_line(image, (2, 5), (8, 5), order=0)
+    expected_prof = np.arange(25, 95, 10)
+    assert_equal(prof, expected_prof)
+
+
+def test_vertical_upward():
+    prof = profile_line(image, (8, 5), (2, 5), order=0)
+    expected_prof = np.arange(85, 15, -10)
+    assert_equal(prof, expected_prof)
+
+
+def test_45deg_right_downward():
+    prof = profile_line(image, (2, 2), (8, 8), order=0)
+    expected_prof = np.array([22, 33, 33, 44, 55, 55, 66, 77, 77, 88])
+    # repeats are due to aliasing using nearest neighbor interpolation.
+    # to see this, imagine a diagonal line with markers every unit of
+    # length traversing a checkerboard pattern of squares also of unit
+    # length. Because the line is diagonal, sometimes more than one
+    # marker will fall on the same checkerboard box.
+    assert_almost_equal(prof, expected_prof)
+
+
+def test_45deg_right_downward_interpolated():
+    prof = profile_line(image, (2, 2), (8, 8), order=1)
+    expected_prof = np.linspace(22, 88, 10)
+    assert_almost_equal(prof, expected_prof)
+
+
+def test_45deg_right_upward():
+    prof = profile_line(image, (8, 2), (2, 8), order=1)
+    expected_prof = np.arange(82, 27, -6)
+    assert_almost_equal(prof, expected_prof)
+
+
+def test_45deg_left_upward():
+    prof = profile_line(image, (8, 8), (2, 2), order=1)
+    expected_prof = np.arange(88, 21, -22. / 3)
+    assert_almost_equal(prof, expected_prof)
+
+
+def test_45deg_left_downward():
+    prof = profile_line(image, (2, 8), (8, 2), order=1)
+    expected_prof = np.arange(28, 83, 6)
+    assert_almost_equal(prof, expected_prof)
+
+
+def test_pythagorean_triangle_right_downward():
+    prof = profile_line(image, (1, 1), (7, 9), order=0)
+    expected_prof = np.array([11, 22, 23, 33, 34, 45, 56, 57, 67, 68, 79])
+    assert_equal(prof, expected_prof)
+
+
+def test_pythagorean_triangle_right_downward_interpolated():
+    prof = profile_line(image, (1, 1), (7, 9), order=1)
+    expected_prof = np.linspace(11, 79, 11)
+    assert_almost_equal(prof, expected_prof)
+
+pyth_image = np.zeros((6, 7), np.float)
+line = ((1, 2, 2, 3, 3, 4), (1, 2, 3, 3, 4, 5))
+below = ((2, 2, 3, 4, 4, 5), (0, 1, 2, 3, 4, 4))
+above = ((0, 1, 1, 2, 3, 3), (2, 2, 3, 4, 5, 6))
+pyth_image[line] = 1.8
+pyth_image[below] = 0.6
+pyth_image[above] = 0.6
+
+
+def test_pythagorean_triangle_right_downward_linewidth():
+    prof = profile_line(pyth_image, (1, 1), (4, 5), linewidth=3, order=0)
+    expected_prof = np.ones(6)
+    assert_almost_equal(prof, expected_prof)
+
+
+def test_pythagorean_triangle_right_upward_linewidth():
+    prof = profile_line(pyth_image[::-1, :], (4, 1), (1, 5),
+                        linewidth=3, order=0)
+    expected_prof = np.ones(6)
+    assert_almost_equal(prof, expected_prof)
+
+
+def test_pythagorean_triangle_transpose_left_down_linewidth():
+    prof = profile_line(pyth_image.T[:, ::-1], (1, 4), (5, 1),
+                        linewidth=3, order=0)
+    expected_prof = np.ones(6)
+    assert_almost_equal(prof, expected_prof)
+
+
+if __name__ == "__main__":
+    from numpy.testing import run_module_suite
+    run_module_suite()
+

skimage/viewer/plugins/lineprofile.py

@@ -1,8 +1,9 @@
 import warnings
 
 import numpy as np
-import scipy.ndimage as ndi
 from skimage.util.dtype import dtype_range
+from skimage import draw
+from skimage import measure
 
 from .plotplugin import PlotPlugin
 from ..canvastools import ThickLineTool
@@ -70,7 +71,11 @@ def attach(self, image_viewer):
                                        on_change=self.line_changed)
         self.line_tool.end_points = np.transpose([x, y])
 
-        scan_data = profile_line(image, self.line_tool.end_points)
+        scan_data = measure.profile_line(image, 
+                                         *self.line_tool.end_points[:, ::-1])
+        self.scan_data = scan_data
+        if scan_data.ndim == 1:
+            scan_data = scan_data[:, np.newaxis]
 
         self.reset_axes(scan_data)
 
@@ -104,8 +109,12 @@ def _autoscale_view(self):
     def line_changed(self, end_points):
         x, y = np.transpose(end_points)
         self.line_tool.end_points = end_points
-        scan = profile_line(self.image_viewer.original_image, end_points,
-                            linewidth=self.line_tool.linewidth)
+        scan = measure.profile_line(self.image_viewer.original_image,
+                                    *end_points[:, ::-1],
+                                    linewidth=self.line_tool.linewidth)
+        self.scan_data = scan
+        if scan.ndim == 1:
+            scan = scan[:, np.newaxis]
 
         if scan.shape[1] != len(self.profile):
             self.reset_axes(scan)
@@ -131,79 +140,20 @@ def reset_axes(self, scan_data):
                                         scan_data[:, 1], 'g-',
                                         scan_data[:, 2], 'b-')
 
+    def output(self):
+        """Return the drawn line and the resulting scan.
 
-def _calc_vert(img, x1, x2, y1, y2, linewidth):
-    # Quick calculation if perfectly horizontal
-    pixels = img[min(y1, y2): max(y1, y2) + 1,
-                 x1 - linewidth / 2: x1 + linewidth / 2 + 1]
-
-    # Reverse index if necessary
-    if y2 > y1:
-        pixels = pixels[::-1, :]
-
-    return pixels.mean(axis=1)[:, np.newaxis]
-
-
-def profile_line(img, end_points, linewidth=1):
-    """Return the intensity profile of an image measured along a scan line.
+        Returns
+        -------
+        line_image : (M, N) uint8 array, same shape as image
+            An array of 0s with the scanned line set to 255.
+        scan : (P,) or (P, 3) array of int or float
+            The line scan values across the image.
+        """
+        (x1, y1), (x2, y2) = self.line_tool.end_points
+        line_image = np.zeros(self.image_viewer.original_image.shape[:2],
+                              np.uint8)
+        rr, cc = draw.line(y1, x1, y2, x2)
+        line_image[rr, cc] = 255
+        return line_image, self.scan_data
 
-    Parameters
-    ----------
-    img : 2d or 3d array
-        The image, in grayscale (2d) or RGB (3d) format.
-    end_points: (2, 2) list
-        End points ((x1, y1), (x2, y2)) of scan line.
-    linewidth: int
-        Width of the scan, perpendicular to the line
-
-    Returns
-    -------
-    return_value : array
-        The intensity profile along the scan line. The length of the profile
-        is the ceil of the computed length of the scan line.
-    """
-    point1, point2 = end_points
-    x1, y1 = point1 = np.asarray(point1, dtype=float)
-    x2, y2 = point2 = np.asarray(point2, dtype=float)
-    dx, dy = point2 - point1
-    channels = 1
-    if img.ndim == 3:
-        channels = 3
-
-    # Quick calculation if perfectly vertical; shortcuts div0 error
-    if x1 == x2:
-        if channels == 1:
-            img = img[:, :, np.newaxis]
-
-        img = np.rollaxis(img, -1)
-        intensities = np.hstack([_calc_vert(im, x1, x2, y1, y2, linewidth)
-                                 for im in img])
-        return intensities
-
-    theta = np.arctan2(dy, dx)
-    a = dy / dx
-    b = y1 - a * x1
-    length = np.hypot(dx, dy)
-
-    line_x = np.linspace(x2, x1, np.ceil(length))
-    line_y = line_x * a + b
-    y_width = abs(linewidth * np.cos(theta) / 2)
-    perp_ys = np.array([np.linspace(yi - y_width,
-                                    yi + y_width, linewidth) for yi in line_y])
-    perp_xs = - a * perp_ys + (line_x + a * line_y)[:, np.newaxis]
-
-    perp_lines = np.array([perp_ys, perp_xs])
-    if img.ndim == 3:
-        pixels = [ndi.map_coordinates(img[..., i], perp_lines)
-                  for i in range(3)]
-        pixels = np.transpose(np.asarray(pixels), (1, 2, 0))
-    else:
-        pixels = ndi.map_coordinates(img, perp_lines)
-        pixels = pixels[..., np.newaxis]
-
-    intensities = pixels.mean(axis=1)
-
-    if intensities.ndim == 1:
-        return intensities[..., np.newaxis]
-    else:
-        return intensities