import numpy as np
rng = np.random.default_rng(4)
H, W = 14, 24
img = np.full((H, W), 60.0)
img[:, 13:] = 190.0 # a strong vertical edge
img[7:, :13] += 26.0 # a FAINT horizontal edge, touching it
img[11, 17] += 45.0 # one isolated bright speck
img += rng.normal(0, 3.5, (H, W))
img = np.clip(img, 0, 255)
FLAT = (slice(9, 14), slice(3, 12)) # a genuinely featureless patch
def show(a, label, hi):
ramp = " .:-=+*#%@"
print(" %s" % label)
for row in a:
print(" " + "".join(
ramp[min(9, int(9 * max(0.0, v) / (hi + 1e-9)))] for v in row))
show(img, "the image", 255)
print(" a strong vertical step at column 13 (60 -> 190), a faint")
print(" horizontal step at row 7 (60 -> 86) that runs into it, and one")
print(" isolated bright speck in the right-hand region. the faint step")
print(" is only about twice the noise, deliberately.")
print()
print("STEP 1 -- the crudest derivative: the difference with the neighbour.")
dx = np.abs(np.diff(img, axis=1))
dy = np.abs(np.diff(img, axis=0))
strong_raw = dx[3, 12]
faint_raw = dy[6, 7]
noise_raw = max(dx[FLAT].max(), dy[FLAT].max())
print(" across the strong edge: %.1f" % strong_raw)
print(" across the faint edge: %.1f" % faint_raw)
print(" worst noise in the flat: %.1f" % noise_raw)
print(" the strong edge is %.1fx the noise floor and easy. the faint"
% (strong_raw / noise_raw))
print(" edge is %.2fx it -- indistinguishable. a single-pixel"
% (faint_raw / noise_raw))
print(" difference amplifies noise as eagerly as it finds edges,")
print(" because differentiation is a high-pass operation and noise is")
print(" entirely high frequency.")
print()
print("STEP 2 -- SOBEL. smooth ALONG the edge while differentiating ACROSS")
print("it. that is the only reason the kernel is 3x3 and not 1x2:")
KX = np.array([[-1., 0, 1], [-2, 0, 2], [-1, 0, 1]])
KY = KX.T
print(" Kx = [-1 0 1] Ky = [-1 -2 -1]")
print(" [-2 0 2] [ 0 0 0]")
print(" [-1 0 1] [ 1 2 1]")
print(" each is an outer product of [1 2 1] and [-1 0 1]. the [1 2 1]")
print(" half averages 3 pixels along the edge; the [-1 0 1] half")
print(" differentiates across it.")
def conv(a, k):
p = np.pad(a, 1, mode="edge")
return np.array([[float((p[i:i + 3, j:j + 3] * k).sum())
for j in range(a.shape[1])] for i in range(a.shape[0])])
gx, gy = conv(img, KX), conv(img, KY)
mag = np.hypot(gx, gy)
s_sob, f_sob = mag[3, 13], mag[7, 6]
n_sob = mag[FLAT].max()
print(" the same three measurements, now with Sobel:")
print("%26s %12s %12s" % ("", "raw diff", "Sobel"))
print("%26s %12.1f %12.1f" % ("strong edge", strong_raw, s_sob))
print("%26s %12.1f %12.1f" % ("faint edge", faint_raw, f_sob))
print("%26s %12.1f %12.1f" % ("worst noise in the flat", noise_raw, n_sob))
print("%26s %11.2fx %11.2fx"
% ("faint-edge margin", faint_raw / noise_raw, f_sob / n_sob))
print(" Sobel has a gain of 8, so compare the RATIOS, not the values.")
print(" the faint edge went from %.2fx the noise -- lost -- to %.2fx"
% (faint_raw / noise_raw, f_sob / n_sob))
print(" it. averaging 3 pixels along a line cuts uncorrelated noise")
print(" by about sqrt(3) while leaving a straight edge untouched,")
print(" and a real edge is straight for far more than 3 pixels.")
print()
print("STEP 3 -- DIRECTION. gx and gy together give an angle, not just a")
print("strength:")
ang = np.degrees(np.arctan2(gy, gx))
print("%22s %10s %10s %12s" % ("location", "gx", "gy", "angle"))
for label, (r, c) in (("the strong edge", (3, 13)),
("the faint edge", (7, 6)),
("flat interior", (11, 6))):
print("%22s %10.1f %10.1f %11.0f deg" % (label, gx[r, c], gy[r, c],
ang[r, c]))
print(" the vertical edge has a big gx and almost no gy; the horizontal")
print(" edge is the reverse. the gradient points ACROSS an edge, never")
print(" along it, so an edge at 0 degrees of gradient is a vertical line.")
print(" over the flat patch the angle is meaningless -- it is the")
print(" direction of the noise, and it changes every pixel.")
print()
print("STEP 4 -- NON-MAXIMUM SUPPRESSION. Sobel answers with a ridge two")
print("or three pixels wide; an edge is one pixel wide. keep a pixel only")
print("if it beats its two neighbours ALONG ITS OWN GRADIENT:")
print(" row 7, columns 10-16 before: "
+ " ".join("%4.0f" % v for v in mag[7][10:17]))
nms = np.zeros_like(mag)
for i in range(1, H - 1):
for j in range(1, W - 1):
a = ang[i, j] % 180
if a < 22.5 or a >= 157.5:
n1, n2 = mag[i, j - 1], mag[i, j + 1]
elif a < 67.5:
n1, n2 = mag[i - 1, j + 1], mag[i + 1, j - 1]
elif a < 112.5:
n1, n2 = mag[i - 1, j], mag[i + 1, j]
else:
n1, n2 = mag[i - 1, j - 1], mag[i + 1, j + 1]
if mag[i, j] >= n1 and mag[i, j] >= n2:
nms[i, j] = mag[i, j]
print(" row 7, columns 10-16 after: "
+ " ".join("%4.0f" % v for v in nms[7][10:17]))
print(" the ridge collapsed onto its crest. this is the first stage")
print(" that needs the ANGLE rather than the strength, which is why")
print(" Sobel bothers to compute two kernels instead of one.")
print(" pixels above 80: %d before, %d after."
% (int((mag > 80).sum()), int((nms > 80).sum())))
print()
print("STEP 5 -- HYSTERESIS. one threshold cannot work here. try three:")
for t in (80, 200, 400):
kept = nms > t
print(" threshold %3d: %3d pixels: %2d on the faint edge, %d on the speck"
% (t, int(kept.sum()), int(kept[6:9, 0:12].sum()),
int(kept[10:13, 16:19].sum())))
faint_peak = nms[6:9, 0:12].max()
speck_peak = nms[10:13, 16:19].max()
print(" and here is why no single number works:")
print(" strongest response on the faint edge: %6.1f" % faint_peak)
print(" strongest response on the speck: %6.1f" % speck_peak)
print(" they are %.0f apart out of %.0f -- the same neighbourhood of"
% (abs(faint_peak - speck_peak), max(faint_peak, speck_peak)))
print(" values. a threshold that keeps one keeps the other, and a")
print(" threshold that drops one drops the other. BRIGHTNESS ALONE")
print(" CANNOT TELL THEM APART. so use TWO thresholds, and let strong")
print(" edges vouch for weak pixels connected to them:")
LO, HI = 80.0, 400.0
strong = nms >= HI
weak = (nms >= LO) & (nms < HI)
keep = strong.copy()
for _ in range(H * W):
grown = keep.copy()
for i in range(1, H - 1):
for j in range(1, W - 1):
if weak[i, j] and not keep[i, j] and keep[i - 1:i + 2, j - 1:j + 2].any():
grown[i, j] = True
if (grown == keep).all():
break
keep = grown
promoted = int(keep.sum() - strong.sum())
print(" strong (>= %.0f): %2d weak (%.0f to %.0f): %2d"
% (HI, int(strong.sum()), LO, HI, int(weak.sum())))
print(" after linking: %d kept. %d weak pixels were promoted because"
% (int(keep.sum()), promoted))
print(" they connect to a strong one, and %d were dropped because"
% int(weak.sum() - promoted))
print(" they connect to nothing.")
show(keep.astype(float), "the final edge map", 1.0)
print(" the faint edge came back along most of its length, with a")
print(" couple of gaps where the noise happened to flatten its crest.")
print(" and it came back from ONE end: not one of its pixels ever")
print(" cleared the high threshold -- all %d strong pixels are on the"
% int(strong.sum()))
print(" vertical line (%d of them on the faint edge). the chain starts"
% int(strong[6:9, 0:12].sum()))
print(" at the junction and walks left. linking is transitive: each")
print(" promoted pixel then vouches for its own neighbour.")
print(" the speck, whose peak response is within %.0f of the faint"
% abs(faint_peak - speck_peak))
print(" edge's, is gone -- %s in the final map."
% ("still present" if keep[10:13, 16:19].any() else "not one pixel of it"))
print(" it had no chain leading back to any evidence. that is the")
print(" whole idea: hysteresis judges a pixel by what it CONNECTS to,")
print(" not by how bright it is.")
print()
print("WHAT EACH STAGE REMOVED, counting pixels marked as edge:")
print("%-32s %8s %14s" % ("stage", "pixels", "speck pixels"))
speck_mask = np.zeros((H, W), bool)
speck_mask[10:13, 16:19] = True
for label, m in (("Sobel magnitude > 80", mag > LO),
("+ non-maximum suppression", nms > LO),
("+ hysteresis linking", keep)):
print("%-32s %8d %14d" % (label, int(m.sum()), int((m & speck_mask).sum())))
print(" thinning cuts the count without discriminating -- it shrinks the")
print(" speck to a crest just as neatly as it shrinks a real edge.")
print(" linking is the only stage that asks whether the evidence is")
print(" connected to anything, and it is the stage that kills the speck.")
print()
print("blur, differentiate, thin, threshold twice: that sequence is Canny,")
print("published in 1986 and still the default in OpenCV. each of the five")
print("stages exists to repair a specific failure of the stage before it,")
print("and none of them does anything a network could not learn -- only")
print("far cheaper, deterministically, with no training data at all.")