How to extract line profile (ray trace line) through 3d matrix (ndarray 3-dim) by specifying initial & final point

Question

This seems like a fairly simple problem, and I haven't been able to find a canned solution. Essentially, this would be the 3d equivalent of skimage.measure profile_line

Consider a matrix A, with dimensions (i,j,k). Each matrix element is a number. In real terms, consider this a voxelized 3 dimensional temperature distribution.

I would like an efficient method to extract a line profile (a ray trace line) through this data from point (i_1,j_1,k_1) to (i_2,j_2,k_2). Or, similarly, defining an initial point (i_1,j_1,k_1) and a line trajectory using theta, phi and radius in polar coordinates.

I recognize an accurate result would consider boundary crossing and partial voxel volumes along the path, but I would be satisfied with a crude approximation which samples nearest voxel values with a regular step size (say, 0.1*voxel dimension) along the ray trace line.

Help is much appreciated. I am happy to describe further as needed,

It would be nice if this worked..

from skimage.measure import profile_line

line = profile_line(3dim_ndarray, (i1,j1,k1), (i2,j2,k2))
print(line)

Answer 1

Just fyi, we have a long-standing PR for 3D profile_line here. You might find some inspiration there.

However, the difficulty and reason why it hasn't been merged yet is that averaging over a cylinder is harder than averaging over a rectangle, as in profile_line (2D). If you don't need averaging, it is reasonably easy to implement with SciPy:

import numpy as np
from scipy import ndimage as ndi, spatial

def profile_line(image, start_coords, end_coords, *,
                 spacing=1, order=0, endpoint=True):
    coords = []
    n_points = int(np.ceil(spatial.distance.euclidean(start_coords, end_coords)
                           / spacing))
    for s, e in zip(start_coords, end_coords):
        coords.append(np.linspace(s, e, n_points, endpoint=endpoint))
    profile = ndi.map_coordinates(image, coords, order=order)
    return profile

order is the order of the interpolation: default 0 above, ie value of nearest voxel, but can be 1 (linear interpolation) or 3 (cubic interpolation). The best choice depends on your needs.

There's also various design choices to be made, e.g. instead of np.ceil you might want to instead calculate what the endpoint should be so that the spacing itself is exact. e.g. if the distance between start and end is 2.5, then your profile might include 3 points and will end 0.5 voxels short of the endpoint you specified. In the implementation above, you will get 3 points between the start and endpoint, but the real spacing will end up being 2.5/3 = 0.8333, instead of 1.