From 011409f66a7e3736b6310e6b1cc12f52da9f9763 Mon Sep 17 00:00:00 2001 From: Juan Nunez-Iglesias Date: Fri, 24 Jan 2014 16:53:14 +1100 Subject: [PATCH 01/19] Add profile_line to measure package The profile_line function is currently part of the skimage LineProfile plugin. However, it's useful in non-interactive contexts, and importing it from the viewer is awkward, mostly hidden, and depends on PyQt for no good reason. By moving the function to `skimage.measure`, it is usable in many more contexts. --- skimage/measure/__init__.py | 4 +- skimage/measure/profile.py | 99 +++++++++++++++++++++++++++++++++++++ 2 files changed, 102 insertions(+), 1 deletion(-) create mode 100644 skimage/measure/profile.py diff --git a/skimage/measure/__init__.py b/skimage/measure/__init__.py index 7bcb166a..d6852107 100755 --- a/skimage/measure/__init__.py +++ b/skimage/measure/__init__.py @@ -4,6 +4,7 @@ from ._regionprops import regionprops, perimeter 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 @@ __all__ = ['find_contours', '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 00000000..c467cb94 --- /dev/null +++ b/skimage/measure/profile.py @@ -0,0 +1,99 @@ +import numpy as np +import scipy.ndimage as ndi + +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, mode='constant', cval=0.0): + """Return the intensity profile of an image measured along a scan line. + + 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, optional + Width of the scan, perpendicular to the line + 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, ((1, 2), (5, 2))) + array([[ 1.], + [ 1.], + [ 2.], + [ 2.]]) + """ + 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(x1, x2, 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, + mode=mode, cval=cval) for i in range(3)] + pixels = np.transpose(np.asarray(pixels), (1, 2, 0)) + else: + pixels = ndi.map_coordinates(img, perp_lines, mode=mode, cval=cval) + pixels = pixels[..., np.newaxis] + + intensities = pixels.mean(axis=1) + + if intensities.ndim == 1: + return intensities[..., np.newaxis] + else: + return intensities + From e24ae388147bd989d863bcb703dd3f1dbf7bfe72 Mon Sep 17 00:00:00 2001 From: Juan Nunez-Iglesias Date: Sun, 26 Jan 2014 00:55:52 +1100 Subject: [PATCH 02/19] Use the new profile_line function from `measure` Note that the indices were inverted relative to previous use, which I believe was incorrect. (See example in docstring, which works, where it used to be inverted.) --- skimage/viewer/plugins/lineprofile.py | 72 ++------------------------- 1 file changed, 5 insertions(+), 67 deletions(-) diff --git a/skimage/viewer/plugins/lineprofile.py b/skimage/viewer/plugins/lineprofile.py index 0d555eb5..30a0cf21 100644 --- a/skimage/viewer/plugins/lineprofile.py +++ b/skimage/viewer/plugins/lineprofile.py @@ -3,6 +3,7 @@ import warnings import numpy as np import scipy.ndimage as ndi from skimage.util.dtype import dtype_range +from skimage import measure from .plotplugin import PlotPlugin from ..canvastools import ThickLineTool @@ -70,7 +71,7 @@ class LineProfile(PlotPlugin): 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) self.reset_axes(scan_data) @@ -104,8 +105,9 @@ class LineProfile(PlotPlugin): 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, + linewidth=self.line_tool.linewidth) if scan.shape[1] != len(self.profile): self.reset_axes(scan) @@ -143,67 +145,3 @@ def _calc_vert(img, x1, x2, y1, y2, linewidth): 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. - - 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 From 7932db7832e62e2f62faa8686e32358d76959bbf Mon Sep 17 00:00:00 2001 From: Juan Nunez-Iglesias Date: Sun, 26 Jan 2014 13:17:31 +1100 Subject: [PATCH 03/19] Change profile line API to use src and dst points --- skimage/measure/profile.py | 17 +++++++++-------- 1 file changed, 9 insertions(+), 8 deletions(-) diff --git a/skimage/measure/profile.py b/skimage/measure/profile.py index c467cb94..8f302473 100644 --- a/skimage/measure/profile.py +++ b/skimage/measure/profile.py @@ -13,15 +13,17 @@ def _calc_vert(img, x1, x2, y1, y2, linewidth): return pixels.mean(axis=1)[:, np.newaxis] -def profile_line(img, end_points, linewidth=1, mode='constant', cval=0.0): +def profile_line(img, src, dst, linewidth=1, mode='constant', cval=0.0): """Return the intensity profile of an image measured along a scan line. 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. + 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 mode : string, one of {'constant', 'nearest', 'reflect', 'wrap'}, optional @@ -45,16 +47,15 @@ def profile_line(img, end_points, linewidth=1, mode='constant', cval=0.0): [1, 1, 1, 2, 2, 2], [1, 1, 1, 2, 2, 2], [0, 0, 0, 0, 0, 0]]) - >>> profile_line(img, ((1, 2), (5, 2))) + >>> profile_line(img, (1, 2), (5, 2)) array([[ 1.], [ 1.], [ 2.], [ 2.]]) """ - point1, point2 = end_points - x1, y1 = point1 = np.asarray(point1, dtype=float) - x2, y2 = point2 = np.asarray(point2, dtype=float) - dx, dy = point2 - point1 + x1, y1 = src = np.asarray(src, dtype=float) + x2, y2 = dst = np.asarray(dst, dtype=float) + dx, dy = dst - src channels = 1 if img.ndim == 3: channels = 3 From 0fd19e5708f09702d1feb50a66648d8d4d61689e Mon Sep 17 00:00:00 2001 From: Juan Nunez-Iglesias Date: Sun, 26 Jan 2014 15:08:51 +1100 Subject: [PATCH 04/19] Update lineprofile plugin to use new API --- skimage/viewer/plugins/lineprofile.py | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/skimage/viewer/plugins/lineprofile.py b/skimage/viewer/plugins/lineprofile.py index 30a0cf21..5cc70d0c 100644 --- a/skimage/viewer/plugins/lineprofile.py +++ b/skimage/viewer/plugins/lineprofile.py @@ -1,7 +1,6 @@ import warnings import numpy as np -import scipy.ndimage as ndi from skimage.util.dtype import dtype_range from skimage import measure @@ -71,7 +70,7 @@ class LineProfile(PlotPlugin): on_change=self.line_changed) self.line_tool.end_points = np.transpose([x, y]) - scan_data = measure.profile_line(image, self.line_tool.end_points) + scan_data = measure.profile_line(image, *self.line_tool.end_points) self.reset_axes(scan_data) @@ -106,7 +105,7 @@ class LineProfile(PlotPlugin): x, y = np.transpose(end_points) self.line_tool.end_points = end_points scan = measure.profile_line(self.image_viewer.original_image, - end_points, + *end_points, linewidth=self.line_tool.linewidth) if scan.shape[1] != len(self.profile): From a6ba745807a16cf451ef13346451b37043b2027c Mon Sep 17 00:00:00 2001 From: Juan Nunez-Iglesias Date: Mon, 27 Jan 2014 11:41:06 +1100 Subject: [PATCH 05/19] rename x-y to col-row to be less ambiguous --- skimage/measure/profile.py | 23 ++++++++++++----------- 1 file changed, 12 insertions(+), 11 deletions(-) diff --git a/skimage/measure/profile.py b/skimage/measure/profile.py index 8f302473..b19e1f74 100644 --- a/skimage/measure/profile.py +++ b/skimage/measure/profile.py @@ -47,35 +47,36 @@ def profile_line(img, src, dst, linewidth=1, mode='constant', cval=0.0): [1, 1, 1, 2, 2, 2], [1, 1, 1, 2, 2, 2], [0, 0, 0, 0, 0, 0]]) - >>> profile_line(img, (1, 2), (5, 2)) + >>> profile_line(img, (2, 1), (2, 5)) array([[ 1.], [ 1.], [ 2.], [ 2.]]) """ - x1, y1 = src = np.asarray(src, dtype=float) - x2, y2 = dst = np.asarray(dst, dtype=float) - dx, dy = dst - src + src_row, src_col = src = np.asarray(src, dtype=float) + dst_row, dst_col = dst = np.asarray(dst, dtype=float) + d_col, d_row = dst - src channels = 1 if img.ndim == 3: channels = 3 # Quick calculation if perfectly vertical; shortcuts div0 error - if x1 == x2: + if src_col == dst_col: if channels == 1: img = img[:, :, np.newaxis] img = np.rollaxis(img, -1) - intensities = np.hstack([_calc_vert(im, x1, x2, y1, y2, linewidth) + intensities = np.hstack([_calc_vert(im, src_col, dst_col, + src_row, dst_row, linewidth) for im in img]) return intensities - theta = np.arctan2(dy, dx) - a = dy / dx - b = y1 - a * x1 - length = np.hypot(dx, dy) + theta = np.arctan2(d_row, d_col) + a = d_row / d_col + b = src_row - a * src_col + length = np.hypot(d_row, d_col) - line_x = np.linspace(x1, x2, np.ceil(length)) + line_x = np.linspace(src_col, dst_col, 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, From 1170d87ba442acee7187f32b019725cbc261f0fa Mon Sep 17 00:00:00 2001 From: Juan Nunez-Iglesias Date: Mon, 27 Jan 2014 15:18:06 +1100 Subject: [PATCH 06/19] Update _calc_vert for simplicity --- skimage/measure/profile.py | 15 +++++---------- 1 file changed, 5 insertions(+), 10 deletions(-) diff --git a/skimage/measure/profile.py b/skimage/measure/profile.py index b19e1f74..3912fc4e 100644 --- a/skimage/measure/profile.py +++ b/skimage/measure/profile.py @@ -1,15 +1,10 @@ import numpy as np import scipy.ndimage as ndi -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, :] - +def _calc_vert(img, col, src_row, dst_row, linewidth): + # Quick calculation if perfectly vertical + pixels = img[src_row:dst_row:np.sign(dst_row - src_row), + col - linewidth / 2: col + linewidth / 2 + 1] return pixels.mean(axis=1)[:, np.newaxis] @@ -66,7 +61,7 @@ def profile_line(img, src, dst, linewidth=1, mode='constant', cval=0.0): img = img[:, :, np.newaxis] img = np.rollaxis(img, -1) - intensities = np.hstack([_calc_vert(im, src_col, dst_col, + intensities = np.hstack([_calc_vert(im, src_col, src_row, dst_row, linewidth) for im in img]) return intensities From ecb50571b10c0330ca42169b645ad4e6d809cb00 Mon Sep 17 00:00:00 2001 From: Juan Nunez-Iglesias Date: Mon, 27 Jan 2014 15:36:46 +1100 Subject: [PATCH 07/19] Return 1D array for single channel images --- skimage/measure/profile.py | 22 ++++++---------------- 1 file changed, 6 insertions(+), 16 deletions(-) diff --git a/skimage/measure/profile.py b/skimage/measure/profile.py index 3912fc4e..27724039 100644 --- a/skimage/measure/profile.py +++ b/skimage/measure/profile.py @@ -5,7 +5,7 @@ def _calc_vert(img, col, src_row, dst_row, linewidth): # Quick calculation if perfectly vertical pixels = img[src_row:dst_row:np.sign(dst_row - src_row), col - linewidth / 2: col + linewidth / 2 + 1] - return pixels.mean(axis=1)[:, np.newaxis] + return pixels.mean(axis=1)[..., np.newaxis] def profile_line(img, src, dst, linewidth=1, mode='constant', cval=0.0): @@ -43,28 +43,21 @@ def profile_line(img, src, dst, linewidth=1, mode='constant', cval=0.0): [1, 1, 1, 2, 2, 2], [0, 0, 0, 0, 0, 0]]) >>> profile_line(img, (2, 1), (2, 5)) - array([[ 1.], - [ 1.], - [ 2.], - [ 2.]]) + array([ 1., 1., 2., 2.]) """ src_row, src_col = src = np.asarray(src, dtype=float) dst_row, dst_col = dst = np.asarray(dst, dtype=float) - d_col, d_row = dst - src - channels = 1 - if img.ndim == 3: - channels = 3 + d_row, d_col = dst - src # Quick calculation if perfectly vertical; shortcuts div0 error if src_col == dst_col: - if channels == 1: + if img.ndim == 2: img = img[:, :, np.newaxis] - img = np.rollaxis(img, -1) intensities = np.hstack([_calc_vert(im, src_col, src_row, dst_row, linewidth) for im in img]) - return intensities + return np.squeeze(intensities) theta = np.arctan2(d_row, d_col) a = d_row / d_col @@ -89,8 +82,5 @@ def profile_line(img, src, dst, linewidth=1, mode='constant', cval=0.0): intensities = pixels.mean(axis=1) - if intensities.ndim == 1: - return intensities[..., np.newaxis] - else: - return intensities + return intensities From 5a12ef692d6ed7498fee4c6ac6886ba9010d0633 Mon Sep 17 00:00:00 2001 From: Juan Nunez-Iglesias Date: Mon, 27 Jan 2014 15:39:09 +1100 Subject: [PATCH 08/19] Allow multichannel images with >3 channels --- skimage/measure/profile.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/skimage/measure/profile.py b/skimage/measure/profile.py index 27724039..04916fe2 100644 --- a/skimage/measure/profile.py +++ b/skimage/measure/profile.py @@ -14,7 +14,7 @@ def profile_line(img, src, dst, linewidth=1, mode='constant', cval=0.0): Parameters ---------- img : 2d or 3d array - The image, in grayscale (2d) or RGB (3d) format. + The image, in grayscale (2d) or multichannel (2d + c) format. 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) @@ -73,8 +73,8 @@ def profile_line(img, src, dst, linewidth=1, mode='constant', cval=0.0): perp_lines = np.array([perp_ys, perp_xs]) if img.ndim == 3: - pixels = [ndi.map_coordinates(img[..., i], perp_lines, - mode=mode, cval=cval) for i in range(3)] + pixels = [ndi.map_coordinates(img[..., i], perp_lines, mode=mode, + cval=cval) for i in range(img.shape[2])] pixels = np.transpose(np.asarray(pixels), (1, 2, 0)) else: pixels = ndi.map_coordinates(img, perp_lines, mode=mode, cval=cval) From 2a6ec20136c6eb2c9c570f29035e3ac6df928e62 Mon Sep 17 00:00:00 2001 From: Juan Nunez-Iglesias Date: Tue, 28 Jan 2014 00:21:50 +1100 Subject: [PATCH 09/19] Fix doctest string and output shape --- skimage/measure/profile.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/skimage/measure/profile.py b/skimage/measure/profile.py index 04916fe2..c84b6329 100644 --- a/skimage/measure/profile.py +++ b/skimage/measure/profile.py @@ -43,7 +43,7 @@ def profile_line(img, src, dst, linewidth=1, mode='constant', cval=0.0): [1, 1, 1, 2, 2, 2], [0, 0, 0, 0, 0, 0]]) >>> profile_line(img, (2, 1), (2, 5)) - array([ 1., 1., 2., 2.]) + array([ 1., 1., 2., 2.]) """ src_row, src_col = src = np.asarray(src, dtype=float) dst_row, dst_col = dst = np.asarray(dst, dtype=float) @@ -78,7 +78,6 @@ def profile_line(img, src, dst, linewidth=1, mode='constant', cval=0.0): pixels = np.transpose(np.asarray(pixels), (1, 2, 0)) else: pixels = ndi.map_coordinates(img, perp_lines, mode=mode, cval=cval) - pixels = pixels[..., np.newaxis] intensities = pixels.mean(axis=1) From ae85e26ee3f81ba0b15d08c6d1d1c593756c5ad9 Mon Sep 17 00:00:00 2001 From: Juan Nunez-Iglesias Date: Wed, 29 Jan 2014 14:23:51 +1100 Subject: [PATCH 10/19] Reparameterize profile_line to avoid the _calc_vert hack --- skimage/measure/profile.py | 51 ++++++++++++++++++-------------------- 1 file changed, 24 insertions(+), 27 deletions(-) diff --git a/skimage/measure/profile.py b/skimage/measure/profile.py index c84b6329..f29cc512 100644 --- a/skimage/measure/profile.py +++ b/skimage/measure/profile.py @@ -1,12 +1,6 @@ import numpy as np import scipy.ndimage as ndi -def _calc_vert(img, col, src_row, dst_row, linewidth): - # Quick calculation if perfectly vertical - pixels = img[src_row:dst_row:np.sign(dst_row - src_row), - col - linewidth / 2: col + linewidth / 2 + 1] - return pixels.mean(axis=1)[..., np.newaxis] - def profile_line(img, src, dst, linewidth=1, mode='constant', cval=0.0): """Return the intensity profile of an image measured along a scan line. @@ -49,29 +43,32 @@ def profile_line(img, src, dst, linewidth=1, mode='constant', cval=0.0): dst_row, dst_col = dst = np.asarray(dst, dtype=float) d_row, d_col = dst - src - # Quick calculation if perfectly vertical; shortcuts div0 error - if src_col == dst_col: - if img.ndim == 2: - img = img[:, :, np.newaxis] - img = np.rollaxis(img, -1) - intensities = np.hstack([_calc_vert(im, src_col, - src_row, dst_row, linewidth) - for im in img]) - return np.squeeze(intensities) + if d_col == 0: + if d_row > 0: + theta = -np.pi / 2 + else: + theta = np.pi / 2 + else: + theta = np.arctan2(-d_row, d_col) - theta = np.arctan2(d_row, d_col) - a = d_row / d_col - b = src_row - a * src_col - length = np.hypot(d_row, d_col) + length = np.ceil(np.hypot(d_row, d_col)) + line_col = np.linspace(src_col, dst_col, length) + line_row = np.linspace(src_row, dst_row, length) - line_x = np.linspace(src_col, dst_col, 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]) + # this if clause is necessary to keep the line centered on the true + # source and destination points. Otherwise, the computed line has + # an offset of `linewidth/2` + if linewidth <= 1: + perp_lines = np.array([line_row[:, np.newaxis], + line_col[:, np.newaxis]]) + else: + col_width = linewidth * np.sin(theta) / 2 + row_width = linewidth * 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 = [ndi.map_coordinates(img[..., i], perp_lines, mode=mode, cval=cval) for i in range(img.shape[2])] From 47f6ddbce02955507c9e065092630d7a2775b447 Mon Sep 17 00:00:00 2001 From: Juan Nunez-Iglesias Date: Wed, 29 Jan 2014 15:06:00 +1100 Subject: [PATCH 11/19] Update interactive lineprofile tool to new API --- skimage/viewer/plugins/lineprofile.py | 21 +++++++-------------- 1 file changed, 7 insertions(+), 14 deletions(-) diff --git a/skimage/viewer/plugins/lineprofile.py b/skimage/viewer/plugins/lineprofile.py index 5cc70d0c..0eed732c 100644 --- a/skimage/viewer/plugins/lineprofile.py +++ b/skimage/viewer/plugins/lineprofile.py @@ -70,7 +70,10 @@ class LineProfile(PlotPlugin): on_change=self.line_changed) self.line_tool.end_points = np.transpose([x, y]) - scan_data = measure.profile_line(image, *self.line_tool.end_points) + scan_data = measure.profile_line(image, + *self.line_tool.end_points[:, ::-1]) + if scan_data.ndim == 1: + scan_data = scan_data[:, np.newaxis] self.reset_axes(scan_data) @@ -105,8 +108,10 @@ class LineProfile(PlotPlugin): x, y = np.transpose(end_points) self.line_tool.end_points = end_points scan = measure.profile_line(self.image_viewer.original_image, - *end_points, + *end_points[:, ::-1], linewidth=self.line_tool.linewidth) + if scan.ndim == 1: + scan = scan[:, np.newaxis] if scan.shape[1] != len(self.profile): self.reset_axes(scan) @@ -132,15 +137,3 @@ class LineProfile(PlotPlugin): scan_data[:, 1], 'g-', scan_data[:, 2], 'b-') - -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] - From 3d93582ae7599165468c621cac0fc01f70239ee8 Mon Sep 17 00:00:00 2001 From: Juan Nunez-Iglesias Date: Wed, 29 Jan 2014 15:13:36 +1100 Subject: [PATCH 12/19] Simplify parameters even more --- skimage/measure/profile.py | 11 ++--------- 1 file changed, 2 insertions(+), 9 deletions(-) diff --git a/skimage/measure/profile.py b/skimage/measure/profile.py index f29cc512..c2834368 100644 --- a/skimage/measure/profile.py +++ b/skimage/measure/profile.py @@ -42,14 +42,7 @@ def profile_line(img, src, dst, linewidth=1, mode='constant', cval=0.0): 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 - - if d_col == 0: - if d_row > 0: - theta = -np.pi / 2 - else: - theta = np.pi / 2 - else: - theta = np.arctan2(-d_row, d_col) + theta = np.arctan2(d_row, d_col) length = np.ceil(np.hypot(d_row, d_col)) line_col = np.linspace(src_col, dst_col, length) @@ -62,7 +55,7 @@ def profile_line(img, src, dst, linewidth=1, mode='constant', cval=0.0): perp_lines = np.array([line_row[:, np.newaxis], line_col[:, np.newaxis]]) else: - col_width = linewidth * np.sin(theta) / 2 + col_width = linewidth * np.sin(-theta) / 2 row_width = linewidth * 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]) From 80dcc0cd7b2244b1193fbff9b697c4f0736b58e7 Mon Sep 17 00:00:00 2001 From: Juan Nunez-Iglesias Date: Wed, 29 Jan 2014 19:08:37 +1100 Subject: [PATCH 13/19] Add 'order' kwarg for interpolation in profile_line() --- skimage/measure/profile.py | 16 +++++++++++----- 1 file changed, 11 insertions(+), 5 deletions(-) diff --git a/skimage/measure/profile.py b/skimage/measure/profile.py index c2834368..008948e2 100644 --- a/skimage/measure/profile.py +++ b/skimage/measure/profile.py @@ -1,8 +1,9 @@ import numpy as np -import scipy.ndimage as ndi +import scipy.ndimage as nd -def profile_line(img, src, dst, linewidth=1, mode='constant', cval=0.0): +def profile_line(img, src, dst, linewidth=1, + order=0, mode='constant', cval=0.0): """Return the intensity profile of an image measured along a scan line. Parameters @@ -15,6 +16,9 @@ def profile_line(img, src, dst, linewidth=1, mode='constant', cval=0.0): 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 @@ -63,11 +67,13 @@ def profile_line(img, src, dst, linewidth=1, mode='constant', cval=0.0): linewidth) for col_i in line_col]) perp_lines = np.array([perp_rows, perp_cols]) if img.ndim == 3: - pixels = [ndi.map_coordinates(img[..., i], perp_lines, mode=mode, - cval=cval) for i in range(img.shape[2])] + 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 = ndi.map_coordinates(img, perp_lines, mode=mode, cval=cval) + pixels = nd.map_coordinates(img, perp_lines, + order=order, mode=mode, cval=cval) intensities = pixels.mean(axis=1) From 00f4d38c2c743cfc3f5b22c286cb3fbce049fda0 Mon Sep 17 00:00:00 2001 From: Juan Nunez-Iglesias Date: Wed, 29 Jan 2014 22:47:52 +1100 Subject: [PATCH 14/19] Polish coordinate system It's important to distinguish between pixels and points. A square set of points belongs to one pixel. When users type `linewidth=3`, they usually mean a line 3 pixels wide. The distance between the pixel center points, then, is 2. --- skimage/measure/profile.py | 20 +++++++++++++++----- 1 file changed, 15 insertions(+), 5 deletions(-) diff --git a/skimage/measure/profile.py b/skimage/measure/profile.py index 008948e2..823bd9da 100644 --- a/skimage/measure/profile.py +++ b/skimage/measure/profile.py @@ -40,27 +40,37 @@ def profile_line(img, src, dst, linewidth=1, [1, 1, 1, 2, 2, 2], [1, 1, 1, 2, 2, 2], [0, 0, 0, 0, 0, 0]]) - >>> profile_line(img, (2, 1), (2, 5)) + >>> 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)) + 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) # this if clause is necessary to keep the line centered on the true # source and destination points. Otherwise, the computed line has - # an offset of `linewidth/2` + # an offset of `linewidth / 2` if linewidth <= 1: perp_lines = np.array([line_row[:, np.newaxis], line_col[:, np.newaxis]]) else: - col_width = linewidth * np.sin(-theta) / 2 - row_width = linewidth * np.cos(theta) / 2 + # 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, From 4bc5d66060b75e63e88a76c8de2f3511d8033546 Mon Sep 17 00:00:00 2001 From: Juan Nunez-Iglesias Date: Wed, 29 Jan 2014 22:50:19 +1100 Subject: [PATCH 15/19] Add numerous (14) tests for profile_line() --- skimage/measure/tests/test_profile.py | 110 ++++++++++++++++++++++++++ 1 file changed, 110 insertions(+) create mode 100644 skimage/measure/tests/test_profile.py diff --git a/skimage/measure/tests/test_profile.py b/skimage/measure/tests/test_profile.py new file mode 100644 index 00000000..e911673b --- /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() + From d71cd4d750600ab871977895fdf13da6a4609b79 Mon Sep 17 00:00:00 2001 From: Juan Nunez-Iglesias Date: Wed, 29 Jan 2014 22:52:49 +1100 Subject: [PATCH 16/19] Remove unnecessary if clause The `if linewidth <= 1` was necessary because the coordinate systems (pixel/point) were all muddled up. Now that it has been clarified, `linewidth=1` works transparently! --- skimage/measure/profile.py | 28 +++++++++++----------------- 1 file changed, 11 insertions(+), 17 deletions(-) diff --git a/skimage/measure/profile.py b/skimage/measure/profile.py index 823bd9da..5a754ee3 100644 --- a/skimage/measure/profile.py +++ b/skimage/measure/profile.py @@ -59,23 +59,17 @@ def profile_line(img, src, dst, linewidth=1, line_col = np.linspace(src_col, dst_col, length) line_row = np.linspace(src_row, dst_row, length) - # this if clause is necessary to keep the line centered on the true - # source and destination points. Otherwise, the computed line has - # an offset of `linewidth / 2` - if linewidth <= 1: - perp_lines = np.array([line_row[:, np.newaxis], - line_col[:, np.newaxis]]) - else: - # 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]) + # 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) From 040a21afb90d5e5c3e6ad43748142395259a680f Mon Sep 17 00:00:00 2001 From: Juan Nunez-Iglesias Date: Wed, 29 Jan 2014 22:55:18 +1100 Subject: [PATCH 17/19] Change default interpolation order to 1 --- skimage/measure/profile.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/skimage/measure/profile.py b/skimage/measure/profile.py index 5a754ee3..e49c710c 100644 --- a/skimage/measure/profile.py +++ b/skimage/measure/profile.py @@ -3,7 +3,7 @@ import scipy.ndimage as nd def profile_line(img, src, dst, linewidth=1, - order=0, mode='constant', cval=0.0): + order=1, mode='constant', cval=0.0): """Return the intensity profile of an image measured along a scan line. Parameters From ec0079dcba74c18afce2feff0e09eb4300e9f039 Mon Sep 17 00:00:00 2001 From: Juan Nunez-Iglesias Date: Wed, 29 Jan 2014 23:27:34 +1100 Subject: [PATCH 18/19] Output line drawing and scan data from plugin --- skimage/viewer/plugins/lineprofile.py | 20 ++++++++++++++++++++ 1 file changed, 20 insertions(+) diff --git a/skimage/viewer/plugins/lineprofile.py b/skimage/viewer/plugins/lineprofile.py index 0eed732c..ea44981a 100644 --- a/skimage/viewer/plugins/lineprofile.py +++ b/skimage/viewer/plugins/lineprofile.py @@ -2,6 +2,7 @@ import warnings import numpy as np from skimage.util.dtype import dtype_range +from skimage import draw from skimage import measure from .plotplugin import PlotPlugin @@ -72,6 +73,7 @@ class LineProfile(PlotPlugin): 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] @@ -110,6 +112,7 @@ class LineProfile(PlotPlugin): 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] @@ -137,3 +140,20 @@ class LineProfile(PlotPlugin): scan_data[:, 1], 'g-', scan_data[:, 2], 'b-') + def output(self): + """Return the drawn line and the resulting scan. + + 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 + From 01967e5060e93e2c91c092b199e0143213357c29 Mon Sep 17 00:00:00 2001 From: Juan Nunez-Iglesias Date: Thu, 30 Jan 2014 12:15:36 +1100 Subject: [PATCH 19/19] Update docstring according to @stefanv's comments --- skimage/measure/profile.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/skimage/measure/profile.py b/skimage/measure/profile.py index e49c710c..98038e06 100644 --- a/skimage/measure/profile.py +++ b/skimage/measure/profile.py @@ -8,8 +8,10 @@ def profile_line(img, src, dst, linewidth=1, Parameters ---------- - img : 2d or 3d array - The image, in grayscale (2d) or multichannel (2d + c) format. + 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)