diff --git a/skimage/measure/__init__.py b/skimage/measure/__init__.py index 7bcb166a2e5..d6852107165 100755 --- a/skimage/measure/__init__.py +++ b/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'] diff --git a/skimage/measure/profile.py b/skimage/measure/profile.py new file mode 100644 index 00000000000..98038e06237 --- /dev/null +++ b/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 + diff --git a/skimage/measure/tests/test_profile.py b/skimage/measure/tests/test_profile.py new file mode 100644 index 00000000000..e911673bf5e --- /dev/null +++ b/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() + diff --git a/skimage/viewer/plugins/lineprofile.py b/skimage/viewer/plugins/lineprofile.py index 0d555eb59bd..ea44981a029 100644 --- a/skimage/viewer/plugins/lineprofile.py +++ b/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