initial commit

README.md

@@ -0,0 +1,75 @@
+# metrology demo
+
+this is a "quick and dirty" concept sketch.
+
+concepts/steps:
+- load picture
+- ensure "object" is white, background is black (maybe invert picture)
+- define a line to sample along
+- sample picture along line to get a 1-D signal
+- possibly lowpass filter against pixel noise
+- calculate gradient
+- find locations of maximum and minimum in gradient (argmax, argmin)
+- difference is distance
+
+
+## metrology
+
+can be taken a lot more seriously than I do here, in terms of precision, usage
+of units, calibration, ...
+
+it's basically "counting pixels". lots of low level image processing and
+geometric calculations.
+
+some industrially known libraries are Cognex "VisionPro" and MVTec "HALCON".
+they were made for "counting pixels" and user friendliness but they grow to
+include Deep Learning too now, inference at least.
+
+OpenCV has nearly no "easy to use" procedures for this purpose. It wasn't made
+for this. it was made for computer vision, a loftier goal. maybe it's a good
+idea to give opencv a metrology module that provides easy to use procedures for
+the things I do here (and more).
+
+## sampling along a line
+
+- two points define the line
+- calculate coordinates along that line
+  - equally spaced points (along X axis)
+  - affine transformation to map them onto the line
+- use some interpolation scheme (nearest neighbor, linear, cubic, higher order)
+  to sample the input image on those points
+
+## to find edges
+
+... there's one approach with two flavors, depending on how smooth your signal
+is: you basically have to find large changes in the signal's level.
+
+you could threshold the signal. then all changes are of the same magnitude:
+1 (binary). easy to find those.
+
+or you could deal with the signal as is. then you calculate a "gradient". then
+you have to find where the gradient is large:
+
+- either "relatively" large, so you find the largest gradient
+- or absolutely "large enough", where you may get multiple adjacent points
+  (hence ambiguity) that satisfy the condition
+
+## what does it look like
+
+there's this petri dish. it contains two plates with a translucent glob between
+them. we need to know how far apart these plates are.  the plates will move
+vertically.
+
+to determine the distance, we're interested in the top edge of the lower plate
+and the bottom edge of the upper plate. we define a line that spans both edges
+comfortably. the rest happens as described above.
+
+screenshot shows samples along line (top) and gradient (bottom):
+
+![screenshot](plot.png)
+
+input picture (not my own picture):
+
+![sample input](dish-1.jpg)
+
+

metrology-demo.py

@@ -0,0 +1,197 @@
+#!/usr/bin/env python3
+
+# written in 2020 by Christoph Rackwitz <christoph.rackwitz@gmail.com>
+# use only for war purposes and sarcasm
+
+import sys
+import numpy as np
+import cv2 as cv
+import scipy.interpolate # interp2d
+import scipy.ndimage # gaussian_filter1d
+
+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) # homography
+
+	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
+
+	# take affine part of homography
+	assert np.isclose(a=H[2], b=(0,0,1)).all() # we didn't touch those but let's better check
+	A = H[0:2, :]
+
+	return (nsamples, A)
+
+def sample_opencv(im, M, nsamples):
+
+	# 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]
+
+	# INTER_CUBIC seems to break down beyond 1/32 sampling (discretizes).
+	# there might be fixed point algorithms at work
+
+	return samples
+
+def sample_scipy(im, M, nsamples):
+	# https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.interp2d.html
+
+	coords = np.vstack([np.arange(nsamples), np.zeros(nsamples), np.ones(nsamples)])
+
+	coords_mapped = M.astype(np.float32) @ coords # @ = np.dot
+
+	# FIXME: interp2d() is an expensive operation if the image is large
+	#        maybe crop to bounding box of line (bbox of coords_mapped)?
+	sampler = scipy.interpolate.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!)
+	samples = sampler(*coords_mapped) # fairly fast compared to building the sampler (interp2d)
+
+	return samples
+
+
+if __name__ == '__main__':
+	do_display = True # see below
+	do_invert = True
+
+	# to remove pixel noise
+	smoothing_sigma = 2 # in pixels
+
+	# define a line segment to sample along
+	p0, p1 = (1320, 2500), (1320, 2100)
+	stride = 1/4 # sample stride in pixels
+
+	# the picture to work with
+	if len(sys.argv) >= 2:
+		imfname = sys.argv[1]
+	else:
+		imfname = "dish-1.jpg"
+
+	########## here be dragons ##########
+
+	decimals = max(0, int(np.ceil(-np.log10(stride))))
+
+	print("loading picture...", end=" ", flush=True)
+	im = cv.imread(imfname, cv.IMREAD_GRAYSCALE)
+	imh, imw = im.shape[:2]
+	if do_invert:
+		im = 255-im # invert
+	im = im.astype(np.float32)# * np.float32(1/255)
+	print("done")
+
+	# build transform
+	nsamples, M = build_transform(p0, p1, stride=stride)
+
+	print(f"taking {nsamples} samples along line {p0} -> {p1}...", end=" ", flush=True)
+
+	# pick one
+	samples = sample_opencv(im, M, nsamples) # does "normal" cubic (4 support points, continuous first derivative)
+	#samples = sample_scipy(im, M, nsamples) # does some fancy cubic with continuous higher derivatives
+
+	print("sampling done")
+
+	# smoothing to remove noise
+	if smoothing_sigma > 0:
+		samples = scipy.ndimage.gaussian_filter1d(samples, sigma=smoothing_sigma / stride)
+
+	# off-by-half in position 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
+
+	print(f"distance: {distance:.{decimals}f} pixels")
+
+	# this was the result. algorithm is done.
+	# now follows displaying code
+
+	if do_display:
+		gradient *= 255 / np.abs(gradient).max()
+
+		# plot signal
+		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 gradient
+		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()
+
+		# arrange vertically
+		canvas = np.vstack([canvas1, canvas2]) # 600 wide, 800 tall
+
+		# draw lines at edges (largest gradients)
+		# 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))
+
+		# some text to describe the picture
+		cv.putText(canvas, f"sampling {p0} -> {p1}",
+			(10, 350), cv.FONT_HERSHEY_SIMPLEX, 0.75, (255,255,255), thickness=1, lineType=cv.LINE_AA)
+
+		cv.putText(canvas, f"stride {stride} px, {nsamples} samples, sigma {smoothing_sigma}",
+			(10, 350+35), cv.FONT_HERSHEY_SIMPLEX, 0.75, (255,255,255), thickness=1, lineType=cv.LINE_AA)
+
+		cv.putText(canvas, f"distance: {distance:.{decimals}f} px",
+			(10, 350+70), cv.FONT_HERSHEY_SIMPLEX, 0.75, (255,255,255), thickness=1, lineType=cv.LINE_AA)
+
+		# save for posterity
+		cv.imwrite("plot.png", canvas)
+
+		cv.imshow("plot", canvas)
+
+		print("press Ctrl+C in the terminal, or press any key while the imshow() window is focused")
+
+		while True:
+			keycode = cv.waitKey(100)
+			if keycode == -1:
+				continue
+			else:
+				print(f"keycode: {keycode}")
+				break
+