i don't think there is an easy way out, but
this isn't particularly hard provided that you know python well.
here is what i was able to whip up in about 2 hours or so (in jupyter, cells are separated with #========):
import numpy as np
import matplotlib.pyplot as plt
# step 1: load image
img = plt.imread("./your_image.png")
img = img[..., 0] # get rid of rgb channels that are present in the sample image
plt.imshow(img)
# ======================
# step 2, identify relevant area, crop the image to it
img_bin = img > 0.05 # threshold image to get where the roi is
# btw, plt.imread gets pixel values to 0-1 range
where_y, where_x = np.where(img_bin)
xmi, xma = np.quantile(where_x, (0.01, 0.99)) # better for images noise
ymi, yma = np.quantile(where_y, (0.01, 0.99))
# show visualize the region of interest
plt.imshow(img_bin)
plt.gca().add_artist(plt.Rectangle((xmi, ymi), xma-xmi, yma-ymi, facecolor='none', edgecolor='red'))
img2 = img[int(ymi):int(yma), int(xmi):int(xma)].copy() # crop the image
# ========================
# step 3: find some sort of starting point for the algorithm
ci = img2.ravel().argmax() # get brightest point
width = img2.shape[1]
cy, cx = ci//width, ci%width
plt.imshow(img2)
plt.plot([cx, ], [cy, ], marker="o", color='red', markerfacecolor="none")
def get_line_ends(cx, cy, len_pixels, rads):
l = len_pixels/2
cos = np.cos(rads)
sin = np.sin(rads)
y0 = cy-l*sin
y1 = cy+l*sin
x0 = cx-l*cos
x1 = cx+l*cos
return x0, x1, y0, y1
x0, x1, y0, y1 = get_line_ends(cx, cy, 100, -np.pi/11)
# notice that because y axis is inverted, line will be rotating clockwise, instead of counter-clockwise
print(x0, x1, y0, y1)
plt.plot([x0, x1] , [y0, y1], c='red')
# ===========================
# step 4. line sampling prototype
x0, x1, y0, y1 = get_line_ends(cx, cy, 100, -np.pi/11)
plt.imshow(img2)
print(x0, x1, y0, y1)
# plt.plot([x0, x1] , [y0, y1], c='red')
xs = np.linspace(x0, x1, 100).astype(int)
ys = np.linspace(y0, y1, 100).astype(int)
plt.plot(xs, ys, c="red", ls="none", marker="s", markersize=1)
plt.xlim(x0-5, x1+5)
plt.ylim(y0+5, y1-5) # y is still inverted
# ===============================
# step 5 sample pixels along the line at a bunch of angles, and find correct angle,
# to find direction of your lines
def sample_coordinates_along_a_line(img, cx, cy, len_pixels, rads): # same variable names
x0, x1, y0, y1 = get_line_ends(cx, cy, len_pixels, rads)
xs = np.linspace(x0, x1, int(len_pixels)).astype(int)
ys = np.linspace(y0, y1, int(len_pixels)).astype(int)
return img[ys, xs]
rs = np.linspace(-np.pi, 0, 100)
ms = np.empty_like(rs)
for i, r in enumerate(rs):
sample = sample_coordinates_along_a_line(img2, cx, cy, 100, r)
ms[i] = sample.mean()
r_est = r_estimated = rs[ms.argmax()] # will be nearly identical to our guess, so i don't plot it
plt.plot(rs, ms)
plt.axvline(-np.pi/11, c='red') # our guess, should match a maximum
# =================================
# step 6: sample along perpendicular direction, to identify lines
r_90 = r_est + np.pi/2 # perpendicular,
# since r_est is between -pi and 0, r_90 is always between -pi/2, pi/2
r_90 = r_est + np.pi/2 # perpendicular,
# since r_est is between -pi and 0, r_90 is always between -pi/2, pi/2
def get_line_in_a_box(cx, cy, r, xmi, xma, ymi, yma):
"get line that is inside of rectangular box, that goes through point (cx, cy), at angle r"
is_steep = np.abs(r) > np.pi/4 # if angle is > 45 deg, then line grows faster vertically
if is_steep:
y0 = ymi; y1 = yma
x0 = cx - cy/np.tan(r)
x1 = cx+(yma-cy)/np.tan(r)
else:
x0 = xmi; x1 = xma
y0 = cy - cx*np.tan(r)
y1 = cy + (xma-cx)*np.tan(r)
return x0, x1, y0, y1
plt.imshow(img2)
x0, x1, y0, y1 = get_line_in_a_box(cx, cy, r_est, xmi=0, xma=img2.shape[1]-1, ymi=0, yma=img2.shape[0]-1)
plt.plot([x0, x1] , [y0, y1], c='red') # along lines
x0, x1, y0, y1 = get_line_in_a_box(cx, cy, r_90, xmi=0, xma=img2.shape[1]-1, ymi=0, yma=img2.shape[0]-1)
plt.plot([x0, x1] , [y0, y1], c='red') # perpendicular to lines
# ================================
# now we figure out where peaks are from sampling along perpendicular
plt.figure()
def sample_coordinates_along_a_line2(img, x0, x1, y0, y1):
len_pixels = np.sqrt((x1-x0)**2+(y1-y0)**2)
print(len_pixels)
xs = np.linspace(x0, x1, int(len_pixels)).astype(int)
ys = np.linspace(y0, y1, int(len_pixels)).astype(int)
return img[ys, xs]
plt.figure()
x0, x1, y0, y1 = get_line_in_a_box(cx, cy, r_90, xmi=0, xma=img2.shape[1]-1, ymi=0, yma=img2.shape[0]-1)
sampl = sample_coordinates_along_a_line2(img2, x0, x1, y0, y1)
trend = np.convolve(sampl, [1/100]*100, mode="same")
sampl_detrended = sampl-trend
plt.plot(sampl_detrended)
# ==============================
# step 7: find maxima in detrended sample
# i've found this function somewhere in my processing scripts, it's more generic, but good one
def bool_to_regions2(xs: np.ndarray[bool], which=None) -> np.ndarray[int]:
"""return (start, end) pairs of each continious region as (n, 2) int array (end not inclusive, as usual). example
```
a = np.array([1,0,0,0,0,1,1,1,0,1,1,], dtype=bool)
for b, e in bool_to_regions2(a):
print(b, e)
print("".join(np.where(a,'1','0')))
print(' '*b + '^'+' '*(e-b-1)+'^')
```
set which to True or False to return only regions for these values
"""
heads = np.diff(xs, prepend=~xs[0])
tails = np.diff(xs, append=~xs[-1])
nh = abs(heads.sum()) # if array is 0 and 1 instead of proper boolean
nt = abs(tails.sum())
assert nh == nt, f"A: function `bool_to_regions` {nh=}, {nt=}, nh!=nt"
r = np.stack([np.where(heads)[0], np.where(tails)[0]+1], axis=-1)
if which is None: return r
elif which is True: return r[::2] if bool(xs[0]) else r[1::2]
elif which is False: return r[::2] if not xs[0] else r[1::2]
else: raise Exception("`which` should be True, False or None")
plt.plot(sampl_detrended)
maxima = bool_to_regions2(sampl_detrended>0.05, which=True).mean(axis=-1)
# maxima are positions of your lines along a perpendicular to them
for m in maxima:
plt.axvline(m, color='red', alpha=0.5)
# =======================================
# step 8: project maxima back to image space, by using linear interpolation along the line
plt.imshow(img2)
# remember, x0, y0, x1, y1 are from the perpendicular
line_x_coords = x0 + (x1-x0) * maxima / len(sampl_detrended)
line_y_coords = y0 + (y1-y0) * maxima / len(sampl_detrended)
plt.plot(line_x_coords, line_y_coords, ls="none", c="red", marker="+")
#================================================
# step 9: sample all lines
line_sampls = []
plt.imshow(img2)
x0, x1, y0, y1 = get_line_in_a_box(cx, cy, r_90, xmi=0, xma=img2.shape[1]-1, ymi=0, yma=img2.shape[0]-1)
plt.plot([x0, x1,], [y0, y1], c="red") # perpendicular
for line_number in range(len(line_x_coords)):
lcx = line_x_coords[line_number]
lcy = line_y_coords[line_number]
x0, x1, y0, y1 = get_line_in_a_box(lcx, lcy, r_est, xmi=0, xma=img2.shape[1]-1, ymi=0, yma=img2.shape[0]-1)
if x0 < 0 or x1 > img2.shape[1]-1 or y0 > img2.shape[0]-1 or y1 < 0:
continue # todo: clip line instead of skipping it
else:
plt.plot([x0, x1,], [y0, y1], c="red") # should cover the lines in the img2
sample = sample_coordinates_along_a_line2(img2, x0, x1, y0, y1)
line_sampls.append(sample)
line_sampls= np.stack(line_sampls)
# ===============================================
# this is how intensity samples look along the lines
# it should be easy to compute center of each line
for sampl in line_sampls:
plt.plot(sampl)
note: jupyter starts new plot in each cell. i'm not sure how to get nice visual feedback in a plain .py script. maybe call plt.show() at the end of each cell?
