Csdrhrt

I have a question that I originally posted to Stackoverflow, but got no answers, and was encouraged to re-post here. I am developing a method for “straightening” binary-masks of pictures of carrot roots. I have traced midlines through the masks. For example, I start with a binary mask:

binary-mask

I then use Voronoi diagrams to estimate the midline, as shown in red:

binary-mask with midline

I would now like to find the length of the “slices” that run perpendicular to the midline at each point, which would then constitute the appropriate width of the root at any given point along its length. I am using a method similar to that explained here to find the derivative of the fitted spline connecting points along the midline:

def get_spline_derivative(points):



x_max = points[:, 0].max()

# x_min = points[:, 0].min()

new_length = x_max



x = points[:, 0]

y = points[:, 1]



# new_x = np.linspace(x_min, x_max, new_length)



# new_x = np.linspace(0, 1, 2000)  # this is the one that works

new_x = np.linspace(0, 1, new_length)  # this is the one that works



tck, u = itp.splprep([x, y])

dx, dy = itp.splev(new_x, tck, der=1)



return list(zip(dx, dy))

This produces an enormous array, shortened here (full output here):

[(-559.8864500362733, -105.4032498463578), (-560.2848158832071, -109.774387868758), (-560.6711489905051, -113.82253391400316), (-561.0454493581669, -117.54768798209324), (-561.4077169861929, -120.94985007302827), ..., (-533.1106034956821, 74.87936144117594)]

I am unsure how to then use this output to find the vector perpendicular to this derivative, and then calculate the length of the “bin” of pixels that is within the polygon, running along this vector.

Does anyone have experience with performing this computations using Python?

I think the two big steps are 1. finding the perpendicular lines to the midline, and 2. finding the distance to the edge along these lines.

The first task is relatively straightforward. Suppose the midline is defined by the points $(x_1,y_1), (x_2,y_2), ldots, (x_k,y_k)$. You can find the tangent vector $(Delta x,Delta y) = (x_2-x_1,y_2-y_1)$ at a point $(x_1,y_1)$, and then find the vector orthogonal to this: $ (-Delta y, Delta x) $ and normalize it.

Now, note that you will get $k-1$ perpendicular lines. So instead of using the original data for the midline, it would probably make sense to smooth it (using some spline interpolation) and then sample the spline to get one more point than the number of perpendicular lines you want.

The second task is a bit tricker. It looks like the carrot shape is more or less convex which means you could essentially define the boundary as a bunch of hyperplanes, and then see where the line from task 1 intersects them. I'm sure there are libraries for this, and if you search "intersect line with convex shape" there are lots of results.

Another thought is to find the points on the boundary such that the slope between the point on the midline and the boundary point is as close to possible as the slope of the normal at the point on the midline. Once you have those points, you can then compute the distance to them from the midline point.

If speed is a concern, searching for these points can probably be done using some kind of bisection.

Sample code

We store the boudary points in an array of points called boundary (black points) and the midline in an array of points called midline (red points). The following code identifies the two blue points which approximate the intersection of the normal line (in grey) with the boundary.

The first task is to find the unit tangent at a point. We can do this for all points simultaneously:

slopes = midline[1:] - midline[:-1]

tangents = [1,-1]*slopes[:,::-1]

tangents /= np.linalg.norm(tangents,axis=1)[:,None] 

Now suppose you have a known point p on and a known direction n.
We first compute the direction from p to each boundary point.

vec_to_boundary = (boundary - p)

vec_to_boundary /= np.linalg.norm(slope_to_boundary,axis=1)[:,None]

Now compute the inner product with each of these vectors.

inner_product_normal = slope_to_boundary@n

The largest (in magnitude) inner products will tell us which points on the boundary are closest to the normal line passing through p.

b0 = np.argmax(inner_product_normal)

b1 = np.argmin(inner_product_normal)

Finally, we can approximate the distance to the boundary from p along n as the distance to the boundary points we found.

d1 = np.linalg.norm(boundary[b0] - p)

d2 = np.linalg.norm(boundary[b1] - p)

Note that for this to work well you need that the boundary have enough points that the stored boundary points are "close" to the intersection of the normal line with the "true" boundary. It also assumes the boundary is convex. In practice you will probably have to smooth the midline and boundary before using this type of approach.