How to extract a 1d profile (with integrated width) from a 2D array in an arbitrary direction

Question

i have the following problem: I would like to extract a 1D profile from a 2D array, which is relatively simple. And it is also easy to do this in an arbitrary direction (see here).

But i would like to give the profile a certain width, so that the values perpendicular to the profile are averaged. I managed to do this, but it is extremely slow. Does anyone have a good solution for that?

Thanks!

import numpy as np
import os
import math
import itertools
import matplotlib.pyplot as plt
from matplotlib.patches import Polygon

def closest_point(points, coords):
    min_distances = []
    coords = coords
    for point in points:
        distances = []
        for coord in coords:
            distances.append(np.sqrt((point[0]-coord[0])**2 + (point[1]-coord[1])**2))

        val, idx = min((val, idx) for (idx, val) in enumerate(distances))
        min_distances.append(coords[idx])
    return min_distances

def rect_profile(x0, y0, x1, y1, width):

    xd=x1-x0
    yd=y1-y0
    alpha = (np.angle(xd+1j*yd))

    y00 = y0 - np.cos(math.pi - alpha)*width
    x00 = x0 - np.sin(math.pi - alpha)*width

    y01 = y0 + np.cos(math.pi - alpha)*width
    x01 = x0 + np.sin(math.pi - alpha)*width

    y10 = y1 + np.cos(math.pi - alpha)*width
    x10 = x1 + np.sin(math.pi - alpha)*width

    y11 = y1 - np.cos(math.pi - alpha)*width
    x11 = x1 - np.sin(math.pi - alpha)*width

    vertices = ((y00, x00), (y01, x01), (y10, x10), (y11, x11))
    poly_points = [x00, x01, x10, x11], [y00, y01, y10, y11]
    poly = Polygon(((y00, x00), (y01, x01), (y10, x10), (y11, x11)))

    return poly, poly_points

def averaged_profile(image, x0, y0, x1, y1, width):
    num = np.sqrt((x1-x0)**2 + (y1-y0)**2)
    x, y = np.linspace(x0, x1, num), np.linspace(y0, y1, num)
    coords = list(zip(x, y))

    # Get all points that are in Rectangle
    poly, poly_points = rect_profile(x0, y0, x1, y1, width)
    points_in_poly = []
    for point in itertools.product(range(image.shape[0]), range(image.shape[1])):
        if poly.get_path().contains_point(point, radius=1) == True:
            points_in_poly.append((point[1], point[0]))

    # Finds closest point on line for each point in poly
    neighbour = closest_point(points_in_poly, coords)

    # Add all phase values corresponding to closest point on line
    data = []
    for i in range(len(coords)):
        data.append([])

    for idx in enumerate(points_in_poly):
        index = coords.index(neighbour[idx[0]])
        data[index].append(image[idx[1][1], idx[1][0]])

    # Average data perpendicular to profile
    for i in enumerate(data):
        data[i[0]] = np.nanmean(data[i[0]])

    # Plot
    fig, axes = plt.subplots(figsize=(10, 5), nrows=1, ncols=2)

    axes[0].imshow(image)
    axes[0].plot([poly_points[0][0], poly_points[0][1]], [poly_points[1][0], poly_points[1][1]], 'yellow')
    axes[0].plot([poly_points[0][1], poly_points[0][2]], [poly_points[1][1], poly_points[1][2]], 'yellow')
    axes[0].plot([poly_points[0][2], poly_points[0][3]], [poly_points[1][2], poly_points[1][3]], 'yellow')
    axes[0].plot([poly_points[0][3], poly_points[0][0]], [poly_points[1][3], poly_points[1][0]], 'yellow')
    axes[0].axis('image')
    axes[1].plot(data)
    return data

from scipy.misc import face
img = face(gray=True)
profile = averaged_profile(img, 10, 10, 500, 500, 10)

Hi there! Please post your code ([MVCE](http://stackoverflow.com/help/mcve)), otherwise it is virtually impossible to help with that. — MB-F, Apr 19 '16 at 08:34
Hi,i added the code, which might be a bit messy, but as i said, i am very unhappy with the solution i found. — F. Win, Apr 19 '16 at 08:50

score 3 · Answer 1 · answered Aug 29 '19 at 23:52

As another option, There's now a scipy measure function that does exactly this (get profile between arbitrary points in a 2d array, with optional width specified):skimage.measure.profile_line.

As a big plus, it also lets you specify the interpolation value to use for off-grid locations.

I'm not sure how it compares to the above code though - I know for orthogonal cases it is much much faster (ie order magnitude or more) to use simple array slicing/summing.

Like heltonbiker says, if you're really needing speed (large array and/or many times) it is faster to first rotate the matrix, then just use slicing. A technique I've used in the past is basically his approach, but also first essentially masking the original unrotated array, and then only rotating the portion of the array that is the size of the area of your profile (plus a bit).

The downside with that approach (for speed) is that for the rotation you need to use some form of interpolation, which is generally slow, and to get accurate results you need at least linear (order 1) interpolation. However most of the python library modules (there are at least 3) for array rotation seem fairly optimised.

...However for pure convenience, then profile_line is the way to go

score 2 · Accepted Answer · edited May 23 '17 at 11:52

The main performance hog is the function closest_point. Computing the distances between all points on the line with all points in the rectangle is really slow.

You can speed the function up considerably by projecting all rectangle points onto the line. The projected point is the closest point on the line, so there is no need for computing all distances. Further, by correctly normalizing and rounding the projection (distance from line start) it can be directly used as an index.

def closest_point(points, x0, y0, x1, y1):
    line_direction = np.array([x1 - x0, y1 - y0], dtype=float)
    line_length = np.sqrt(line_direction[0]**2 + line_direction[1]**2)
    line_direction /= line_length

    n_bins = int(np.ceil(line_length))

    # project points on line
    projections = np.array([(p[0] * line_direction[0] + p[1] * line_direction[1]) for p in points])

    # normalize projections so that they can be directly used as indices
    projections -= np.min(projections)
    projections *= (n_bins - 1) / np.max(projections)
    return np.floor(projections).astype(int), n_bins

If you wonder about the strange for inside brackets - these are list comprehensions.

Use the function like this inside averaged_profile:

#...

# Finds closest point on line for each point in poly
neighbours, n_bins = closest_point(points_in_poly, x0, y0, x1, y1)

# Add all phase values corresponding to closest point on line
data = [[] for _ in range(n_bins)]
for idx in enumerate(points_in_poly):
    index = neighbours[idx[0]]
    data[index].append(image[idx[1][1], idx[1][0]])

#...

This optimization will make the computation noticably faster. If it is still too slow for you, you can also optimize how you find the points inside the polygon. Instead of testing if each point in the image is inside the rectangle you can use a polygon rasterization algorithm to directly generate the coordinates. See here for details.

Finally, although it is not a performance issue, the use of complex numbers to compute an angle is very creative :) Instead of trigonometric functions you can use the fact that the normal vector of the line is [yd, -xd] devided by line length:

def rect_profile(x0, y0, x1, y1, width):

    xd = x1 - x0
    yd = y1 - y0
    length = np.sqrt(xd**2 + yd**2)

    y00 = y0 + xd * width / length
    x00 = x0 - xd * width / length

    y01 = y0 - xd * width / length
    x01 = x0 + xd * width / length

    y10 = y1 - xd * width / length
    x10 = x1 + xd * width / length

    y11 = y1 + xd * width / length
    x11 = x1 - xd * width / length

    poly_points = [x00, x01, x10, x11], [y00, y01, y10, y11]
    poly = Polygon(((y00, x00), (y01, x01), (y10, x10), (y11, x11)))

    return poly, poly_points

score 0 · Answer 3 · answered Apr 19 '16 at 12:13

I would do the following:

Find out the direction of the desired line, and thus its angle;
Rotate the 2D array via matrix multiplication, using a rotation matrix with the angle you found;
Do a simple bounding box filtering with a rectangle representing your selected area. This will by definition be aligned to one of the axes;
Discard the y coordinate of the points inside the bounding box;
Smooth the results, possibly via 1D spline interpolation (available in scipy).

How to extract a 1d profile (with integrated width) from a 2D array in an arbitrary direction

3 Answers3