foo2.py

#!/usr/bin/env python3

import sys
import numpy as np
import cv2 as cv
from scipy.interpolate import interp2d

def build_transform(p0, p1, stride=None, nsamples=None):
	"builds an affine transform with x+ along defined line"
	# use one of stride (in pixels) or nsamples (absolute value)

	(x0, y0) = p0
	(x1, y1) = p1

	dx = x1 - x0
	dy = y1 - y0

	length = np.hypot(dx, dy)

	if nsamples is not None:
		#stride = length / nsamples
		factor = 1 / nsamples

	else:
		if stride is None:
			stride = 1.0

		factor = stride / length
		nsamples = int(round(length / stride))

	# map: src <- dst (use WARP_INVERSE_MAP flag for warpAffine)
	H = np.eye(3, dtype=np.float64)

	H[0:2, 0] = (dx, dy) # x unit vector
	H[0:2, 1] = (-dy, dx) # y unit vector is x rotated by 90 degrees

	H[0:2, 0:2] *= factor

	H[0:2, 2] = (x0, y0) # translate onto starting point

	A = H[0:2, :] # take affine part

	# number of samples along line
	return (nsamples, A)

print("loading picture...", end=" ", flush=True)
im = cv.imread(sys.argv[1] if len(sys.argv) >= 2 else "200801-2056.png", cv.IMREAD_GRAYSCALE)
imh, imw = im.shape[:2]
im = 255-im # invert
im = im.astype(np.float32)
print("done")

# define a line segment to sample along
p0, p1 = (81, 143), (81, 108)
stride = 1/64 # sample stride in pixels

# build transform
nsamples, M = build_transform(p0, p1, stride=stride)

print(f"taking {nsamples} samples along line {p0} -> {p1}")

if False:
	# INTER_CUBIC seems to break down beyond 1/32 sampling (discretizes). there might be fixed point algorithms at work

	# use transform to get samples
	samples = cv.warpAffine(im, M=M, dsize=(nsamples, 1), flags=cv.WARP_INVERSE_MAP | cv.INTER_CUBIC)
	# data is a row vector
	samples = samples[0]
	#samples = samples.astype(np.int16)

	# opencv seems to use fixed point arithmetic that produces significant error
	# either that or the implementation is broken

# https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.interp2d.html
sampler = interp2d(x=np.arange(imw), y=np.arange(imh), z=im, kind='cubic')
sampler = np.vectorize(sampler) # doesn't cake coordinate pairs as is, vectorize() handles that (!= execution speed!)
coords = np.vstack([np.arange(nsamples), np.zeros(nsamples), np.ones(nsamples)])
coords2 = M.astype(np.float32) @ coords # @ = np.dot
samples2 = sampler(*coords2)
samples = samples2

# two variants:
# (1) get gradient, find minimum (falling edge) and maximum (rising edge)
# (2) use levels, find first black sample from both ends

# off-by-half because for values [0,1,1,0] this returns [+1,0,-1]
gradient = np.diff(samples) / stride

i_falling = np.argmin(gradient) # in samples
i_rising = np.argmax(gradient) # in samples

distance = (i_rising - i_falling) * stride # in pixels

## np.argmin finds the smallest element, and if there are multiple smallest, the first one
## find first black pixel from one end
## find last white pixel from other end (first black + 1)
#thresholded = (samples >= 128)
#v1 = np.argmin(thresholded)
#v2 = len(thresholded) - np.argmin(thresholded[::-1])

print(f"distance: {distance} pixels")

plot = cv.plot.Plot2d_create(np.arange(nsamples, dtype=np.float64), samples.astype(np.float64))
plot.setMinY(256+32)
plot.setMaxY(-32)
plot.setMinX(0)
plot.setMaxX(nsamples)
plot.setGridLinesNumber(5)
plot.setShowText(False) # callout for specific point, setPointIdxToPrint(index)
plot.setPlotGridColor((64,)*3)
canvas1 = plot.render()

plot = cv.plot.Plot2d_create(np.arange(nsamples-1) + 0.5, gradient.astype(np.float64))
plot.setMinY(256+64)
plot.setMaxY(-256-64)
plot.setMinX(0)
plot.setMaxX(nsamples)
plot.setGridLinesNumber(5)
plot.setShowText(False) # callout for specific point, setPointIdxToPrint(index)
plot.setPlotGridColor((64,)*3)
canvas2 = plot.render()

canvas = np.vstack([canvas1, canvas2]) # 600 wide, 800 tall

# plots are 600x400 pixels... and there's no way to plot multiple or plot lines in "plot space"
px_falling = int(600 * (i_falling+0.5) / nsamples)
px_rising = int(600 * (i_rising+0.5) / nsamples)
cv.line(canvas, (px_falling, 0), (px_falling, 400*2), color=(255,0,0))
cv.line(canvas, (px_rising, 0), (px_rising, 400*2), color=(255,0,0))

cv.imwrite("plot.png", canvas)

cv.imshow("plot", canvas)
cv.waitKey(-1)


if 0:
	# let's draw something
	canvas = cv.cvtColor(im, cv.COLOR_GRAY2BGR)
	canvas = cv.resize(canvas, dsize=None, fx=4, fy=4, interpolation=cv.INTER_NEAREST)
	
	cv.line(canvas, (sx*4+2, sy0*4+2), (sx*4+2, sy1*4+2), color=(0,0,255))
	
	cv.line(canvas, ((sx-2)*4+2, (sy0+v1)*4+2), ((sx+2)*4+2, (sy0+v1)*4+2), color=(0, 255, 255))
	cv.line(canvas, ((sx-2)*4+2, (sy0+v2)*4+2), ((sx+2)*4+2, (sy0+v2)*4+2), color=(255, 0, 0))
	
	cv.imshow("canvas", canvas)
	cv.waitKey(-1)