# Stereo calibration GUI -- split out of stereo_DPIV.py.
# UI only; all numerical work lives in stereo_DPIV.
import tkinter as tk
from tkinter import filedialog
from tkinter import messagebox
from PIL import Image, ImageTk, ImageDraw
import numpy as np
import cv2
import matplotlib.pyplot as plt
from scipy.spatial import KDTree
from dpivsoft.stereo_DPIV import stereo_calibration
[docs]
def read_grayscale_image(path):
"""Load an image as a uint8 grayscale array, robust to >8-bit TIFFs.
OpenCV's grayscale reader returns None on 12-bit TIFFs, so fall back to a
depth-aware read and then scikit-image, rescaling deeper formats to uint8.
"""
img = cv2.imread(path, cv2.IMREAD_GRAYSCALE)
if img is None:
img = cv2.imread(path, cv2.IMREAD_ANYDEPTH | cv2.IMREAD_ANYCOLOR)
if img is None:
from skimage import io as skio
img = skio.imread(path)
if img is None:
raise IOError(f"Could not read image: {path}")
if img.ndim == 3:
img = img[..., :3].mean(axis=2)
if img.dtype != np.uint8:
img = img.astype(np.float64)
peak = img.max()
if peak > 0:
img = img * (255.0 / peak)
img = img.astype(np.uint8)
return img
[docs]
def find_points_TSI(image):
"""Detect the TSI target dots on a binary image.
Returns (cX, cY, center, None): sub-pixel float centroids plus the
centre marker (the biggest blob). Blobs clipped by the image border are
skipped -- their centroids are biased by up to half a dot radius and
poison the calibration fit.
"""
num_labels, labels, stats, centroids = (
cv2.connectedComponentsWithStats(image)
)
H, W = image.shape
cX = []
cY = []
old_area = 0
center = None
# Keep the sub-pixel float centroids: connectedComponentsWithStats
# returns them as floats and truncating to int throws away precision
# that directly limits the calibration fit RMS.
for i in range(1, num_labels):
bx, by, bw, bh, area = stats[i]
if bx <= 2 or by <= 2 or bx + bw >= W - 2 or by + bh >= H - 2:
continue
cX.append(float(centroids[i][0]))
cY.append(float(centroids[i][1]))
# The centre marker is the biggest blob
if area > old_area:
old_area = area
center = [cX[-1], cY[-1]]
return cX, cY, center, None
[docs]
def find_points_LaVision(image):
"""Detect LaVision calibration dots plus the two orientation fiducials.
The square and triangle markers are located by shape, removed from the
returned dot list, and handed back together in ``fiducials``
([square, triangle]). ``center`` is the world origin: the dot next to
the square and on top of the triangle, chosen so both cameras pick the
same physical dot.
"""
# Find connected components: every blob is a dot candidate.
num_labels, labels, stats, centroids = (
cv2.connectedComponentsWithStats(image)
)
H, W = image.shape
cX, cY, extent, border = [], [], [], []
for i in range(1, num_labels):
bx, by, bw, bh, area = stats[i]
# Keep sub-pixel float centroids (see find_points_TSI).
cX.append(float(centroids[i][0]))
cY.append(float(centroids[i][1]))
# bbox extent = area / bounding-box area distinguishes the
# shapes: a filled square fills its box most, a triangle
# least, and a round dot lands in between.
extent.append(area / float(bw * bh))
border.append(bx <= 2 or by <= 2 or
bx + bw >= W - 2 or by + bh >= H - 2)
cX = np.array(cX)
cY = np.array(cY)
extent = np.array(extent)
border = np.array(border)
# Identify the two fiducials by extent, searched only among
# non-border blobs (border-clipped dots have unreliable extents):
# square = largest extent (fills its bounding box most)
# triangle = smallest extent (fills its bounding box least)
interior = np.where(~border)[0]
square_idx = interior[np.argmax(extent[interior])]
triangle_idx = interior[np.argmin(extent[interior])]
# Keep sub-pixel floats: square -> origin is the +X direction
# hint for the lattice sort.
square = [float(cX[square_idx]), float(cY[square_idx])]
triangle = [float(cX[triangle_idx]), float(cY[triangle_idx])]
fiducials = [square, triangle]
# Remove both fiducials from the dot list, and drop border-clipped
# dots: their centroids are biased by up to half a dot radius and
# poison the calibration fit.
keep = ~border
keep[[square_idx, triangle_idx]] = False
cX, cY = cX[keep], cY[keep]
# World origin: the dot next to the square and on top of the
# triangle -- the corner of the L the two markers form -- i.e.
# the dot nearest (triangle_x, square_y).
j = int(np.argmin(np.hypot(cX - triangle[0], cY - square[1])))
center = [float(cX[j]), float(cY[j])]
return list(cX), list(cY), center, fiducials
_UNASSIGNED = np.iinfo(np.int32).min
def _lattice_design(rc, order):
"""Design matrix for the smooth (row, col) -> pixel polynomial map."""
r, c = rc[:, 0], rc[:, 1]
cols = [np.ones_like(r), c, r]
if order >= 2:
cols += [c * c, c * r, r * r]
return np.column_stack(cols)
def _lattice_sort(pts, origin, xdir_hint, pattern):
"""Assign integer lattice coordinates (row, col) to detected target dots.
Predictive grid walk: seed an affine (row, col) -> pixel map from the
origin dot and its unambiguous nearest neighbours, then alternate between
refitting a smooth polynomial map over the assigned dots and re-matching
every dot against the predicted lattice-node positions. Because each dot
is matched by *predicted position* (tolerance ~0.3 of the local pitch)
instead of by counting sorted coordinates, the sort survives target
rotation, perspective/shear (non-perpendicular axes), the dual-plane
parallax stagger, and missing dots -- a hole in the grid stays a hole
instead of shifting every label after it.
pts: (N, 2) float pixel coordinates of the candidate dots
origin: pixel coordinates of the origin dot (the world (0, 0) dot)
xdir_hint: approximate pixel direction of world +X (need not be exact:
it only classifies the origin's nearest neighbours, so any
vector within ~25 deg of the true in-row direction works)
pattern: 'TSI' -- rectangular lattice, all (r, c) valid
'LaVision' -- staggered rows: col in half-pitch units and only
nodes with (r - c) even exist (in-row step is 2)
Returns (rc, origin_idx, report): rc is (N, 2) int with _UNASSIGNED for
dots that could not be labelled; report carries fit diagnostics.
"""
pts = np.asarray(pts, float)
n = len(pts)
lavision = pattern == "LaVision"
step_c = 2 if lavision else 1
origin_idx = int(np.argmin(np.hypot(pts[:, 0] - origin[0],
pts[:, 1] - origin[1])))
xhat = np.asarray(xdir_hint, float)
xhat = xhat / np.linalg.norm(xhat)
# --- seed: classify the origin's nearest neighbours by angle to xhat.
# The angle bands leave gaps on purpose: dots at ambiguous angles (e.g.
# the vertical two-row neighbour, which sits at the same distance as the
# in-row one) are simply not used as seeds -- the growth loop labels them
# later. The diagonal/in-column bands are wide enough to survive the
# dual-plane parallax skew (up to ~20% of a pitch).
d = pts - pts[origin_idx]
dist = np.hypot(d[:, 0], d[:, 1])
order_near = np.argsort(dist)[1:9]
inrow, diag = [], []
for k in order_near:
u = d[k] / dist[k]
cosang = abs(u @ xhat)
if cosang > 0.9:
inrow.append(k)
elif lavision and 0.35 < cosang < 0.85:
diag.append(k) # staggered-row diagonal (~45 deg)
elif not lavision and cosang < 0.35:
diag.append(k) # TSI in-column neighbour (~90 deg)
if not inrow or not diag:
raise RuntimeError(
"Lattice sort: could not identify the origin dot's in-row and "
"row-to-row neighbours. Check the threshold/minimum-size options "
"and that the target is not rotated more than ~25 degrees.")
# Seed one dot per lattice offset (nearest wins): a missing neighbour on
# one side must not let the 2nd dot on the other side steal its label.
seed = {origin_idx: (0, 0)}
taken = {(0, 0)}
for k in inrow[:2]:
sgn_c = 1 if d[k] @ xhat > 0 else -1
lab = (0, step_c * sgn_c)
if lab not in taken:
seed[k] = lab
taken.add(lab)
for k in diag[:2]:
sgn_r = 1 if d[k][1] > 0 else -1 # +row is down-image
if lavision:
sgn_dc = 1 if d[k] @ xhat > 0 else -1
lab = (sgn_r, sgn_dc)
else:
lab = (sgn_r, 0)
if lab not in taken:
seed[k] = lab
taken.add(lab)
assigned = dict(seed) # dot idx -> (r, c)
# --- iterative refit + re-match
for it in range(60):
idxs = np.array(list(assigned.keys()))
rc = np.array([assigned[i] for i in idxs], float)
use_order = 2 if len(idxs) >= 12 else 1
A = _lattice_design(rc, use_order)
coef_x, *_ = np.linalg.lstsq(A, pts[idxs, 0], rcond=None)
coef_y, *_ = np.linalg.lstsq(A, pts[idxs, 1], rcond=None)
# candidate lattice nodes: bounding box of the assigned dots plus a
# 2-node ring, so a missing dot cannot stop growth in its direction
rmin, rmax = int(rc[:, 0].min()) - 2, int(rc[:, 0].max()) + 2
cmin, cmax = (int(rc[:, 1].min()) - 2 * step_c,
int(rc[:, 1].max()) + 2 * step_c)
nodes = [(r, c)
for r in range(rmin, rmax + 1)
for c in range(cmin, cmax + 1)
if not (lavision and (r - c) % 2 != 0)]
nodes = np.array(nodes, float)
An = _lattice_design(nodes, use_order)
pred = np.column_stack([An @ coef_x, An @ coef_y])
# matching tolerance: 0.3 of the local in-row pitch at each node
nodes_p = nodes.copy()
nodes_p[:, 1] += step_c
Anp = _lattice_design(nodes_p, use_order)
pred_p = np.column_stack([Anp @ coef_x, Anp @ coef_y])
tol = 0.3 * np.hypot(*(pred_p - pred).T)
# match nodes to dots, greedy nearest-first, one-to-one
new = {}
D = np.hypot(pred[:, None, 0] - pts[None, :, 0],
pred[:, None, 1] - pts[None, :, 1])
cand = np.argwhere(D < tol[:, None])
cand = cand[np.argsort(D[cand[:, 0], cand[:, 1]])]
used_nodes, used_dots = set(), set()
for ni, di in cand:
if ni in used_nodes or di in used_dots:
continue
used_nodes.add(ni)
used_dots.add(di)
new[di] = (int(nodes[ni, 0]), int(nodes[ni, 1]))
new[origin_idx] = (0, 0) # keep the origin pinned
if new == assigned:
break
assigned = new
rc_out = np.full((n, 2), _UNASSIGNED, dtype=int)
for di, (r, c) in assigned.items():
rc_out[di] = (r, c)
idxs = np.array(list(assigned.keys()))
rc = np.array([assigned[i] for i in idxs], float)
A = _lattice_design(rc, 2 if len(idxs) >= 12 else 1)
coef_x, *_ = np.linalg.lstsq(A, pts[idxs, 0], rcond=None)
coef_y, *_ = np.linalg.lstsq(A, pts[idxs, 1], rcond=None)
res = np.hypot(A @ coef_x - pts[idxs, 0], A @ coef_y - pts[idxs, 1])
report = {'n_assigned': len(assigned), 'n_dots': n,
'fit_rms_px': float(np.sqrt((res ** 2).mean())),
'fit_max_px': float(res.max()), 'iters': it + 1}
return rc_out, origin_idx, report
[docs]
def sort_points(x, y, center, fiducial, pattern, org_height, showImages=1,
camera=None, camera_config=None):
# =========================================================
# Algorithm to sort all points in the image
# =========================================================
x = np.asarray(x)
y = np.asarray(y)
if pattern in ("TSI", "LaVision"):
# Predictive lattice walk (see _lattice_sort): fiducial-anchored,
# rotation/perspective tolerant, and missing dots stay holes instead
# of shifting the labels of every dot after them.
fiducial_idx = None
pts = np.column_stack([x, y]).astype(float)
if pattern == "LaVision":
# World +X direction from the fiducial pair: the origin dot sits
# to the right of the square marker, so square -> origin points
# along +X in the image whatever the target rotation.
square = fiducial[0]
xdir = (center[0] - square[0], center[1] - square[1])
else:
# TSI has no directional fiducial: use the origin's neighbour
# most aligned with image +x (assumes rotation below ~45 deg).
d = pts - np.asarray(center, float)
dd = np.hypot(d[:, 0], d[:, 1])
near = np.argsort(dd)[1:5]
k = near[np.argmax(np.abs(d[near, 0]) / dd[near])]
xdir = d[k] if d[k][0] > 0 else -d[k]
rc, origin_idx, report = _lattice_sort(pts, center, xdir, pattern)
# drop dots the walk could not label (stray blobs, dirt, reflections)
ok = rc[:, 0] > _UNASSIGNED
x = x[ok]
y = y[ok]
rows = rc[ok, 0]
cols = rc[ok, 1]
# TSI keeps its historical column-major order (index = column
# breaks); LaVision is row-major (index = row breaks).
if pattern == "TSI":
order = np.lexsort((rows, cols))
major = cols[order]
else:
order = np.lexsort((cols, rows))
major = rows[order]
x = x[order]
y = y[order]
rows = rows[order]
cols = cols[order]
breaks = np.where(np.diff(major) != 0)[0]
index = [-1] + list(breaks) + [len(x) - 1]
center_idx = int(np.where((rows == 0) & (cols == 0))[0][0])
q_rows = rows
q_cols = cols
if showImages:
fig, ax = plt.subplots()
ax.scatter(x, y, s=10, zorder=3)
for xi, yi, r, c in zip(x, y, rows, cols):
ax.annotate(f"({r},{c})", (xi, yi),
fontsize=5, ha='center', va='bottom')
ax.scatter(x[center_idx], y[center_idx], c='red', s=100,
marker='x', zorder=4, label='Origin dot')
ax.set_title(f'{pattern}: lattice assignment '
f'({report["n_assigned"]}/{report["n_dots"]} dots, '
f'fit RMS {report["fit_rms_px"]:.1f} px)')
ax.set_xlabel('x (px)')
ax.set_ylabel('y (px)')
ax.invert_yaxis()
ax.legend()
plt.show()
elif pattern == "Quadricle":
fiducial_idx = None
# All coordinates in image space (y increases downward).
x_img = x.astype(float)
y_img = y.astype(float)
pts = np.column_stack([x_img, y_img])
n = len(pts)
# Estimate grid pitch from median nearest-neighbour distance
tree = KDTree(pts)
nn_dists = tree.query(pts, k=2)[0][:, 1]
pitch = np.median(nn_dists)
tol = 0.45 * pitch
# Find the 4 intersections that form the corners of the centre square
# (the square whose centre is marked by the larger hole).
hole = np.array([center[0], center[1]], dtype=float)
corner_dists = np.hypot(pts[:, 0] - hole[0], pts[:, 1] - hole[1])
n_corners = min(4, len(pts))
corner_idxs = np.argsort(corner_dists)[:n_corners]
# Select the origin corner based on camera orientation.
# In image space: y_img increases downward, so "top" = min y_img.
# For the back camera in Front-Back config the image is x-mirrored
# w.r.t. the physical plane, so the physical top-right corner appears
# at the top-LEFT of the image. In every other case we want the
# top-right corner of the image.
cx = pts[corner_idxs, 0] # x in image space
cy = pts[corner_idxs, 1] # y in image space (y increases downward)
if camera == "left" and camera_config == "Front-Back":
# Back camera: top-left in image → score = -x - y_img
scores = -cx - cy
else:
# Front camera or side-by-side: top-right in image → score = x - y_img
scores = cx - cy
origin_idx = int(corner_idxs[int(np.argmax(scores))])
# BFS: propagate (row, col) grid coordinates from the origin
grid_rc = {origin_idx: (0, 0)}
queue = [origin_idx]
while queue:
curr = queue.pop(0)
curr_r, curr_c = grid_rc[curr]
curr_pt = pts[curr]
for nb in tree.query_ball_point(curr_pt, pitch + tol):
if nb in grid_rc:
continue
nb_pt = pts[nb]
dx_nb = nb_pt[0] - curr_pt[0]
dy_nb = nb_pt[1] - curr_pt[1]
dist = np.hypot(dx_nb, dy_nb)
if dist < pitch - tol:
continue # too close — not a grid neighbour
# Classify direction (cardinal only, 45° threshold)
if abs(dx_nb) >= abs(dy_nb):
dr, dc = 0, (1 if dx_nb > 0 else -1)
else:
dr, dc = (1 if dy_nb > 0 else -1), 0
grid_rc[nb] = (curr_r + dr, curr_c + dc)
queue.append(nb)
# Sort row-major and reorder the coordinate arrays
assigned = sorted(grid_rc.items(),
key=lambda item: (item[1][0], item[1][1]))
order = [a[0] for a in assigned]
rows = np.array([a[1][0] for a in assigned])
cols = np.array([a[1][1] for a in assigned])
x = x[order]
y = y[order]
# Remove partial boundary rows/columns.
# An outermost row or column is partial if it has fewer detected
# points than its immediate interior neighbour. Repeat until all
# boundaries are as complete as the next row/column inward.
changed = True
while changed:
changed = False
for axis, arr in [('row', rows), ('col', cols)]:
unique = sorted(set(arr.tolist()))
if len(unique) < 2:
continue
for boundary, neighbour in [
(unique[0], unique[1]),
(unique[-1], unique[-2])]:
if (np.sum(arr == boundary) <
np.sum(arr == neighbour)):
keep = arr != boundary
x = x[keep]
y = y[keep]
rows = rows[keep]
cols = cols[keep]
order = [o for o, k in zip(order, keep) if k]
changed = True
break # restart — arr reference is stale
if changed:
break
# Build index array: marks boundaries between rows
row_breaks = np.where(np.diff(rows) != 0)[0]
index = [-1] + list(row_breaks) + [len(x) - 1]
center_idx = (order.index(origin_idx) if origin_idx in order
else int(np.argmin(np.abs(rows) + np.abs(cols))))
# Store absolute grid coords for calibrationMesh
q_rows = rows
q_cols = cols
if showImages:
fig, ax = plt.subplots()
x_plot = x.astype(float)
y_plot = y.astype(float)
ax.scatter(x_plot, y_plot, s=10, zorder=3)
for xi, yi, r, c in zip(x_plot, y_plot, rows, cols):
ax.annotate(f"({r},{c})", (xi, yi),
fontsize=5, ha='center', va='bottom')
ax.scatter(hole[0], hole[1], c='red', s=100,
marker='x', zorder=4, label='Hole (origin)')
ax.set_title('Quadricle: BFS grid assignment')
ax.set_xlabel('x (px)')
ax.set_ylabel('y (px)')
ax.invert_yaxis()
ax.legend()
plt.show()
else:
raise ValueError(f"Unknown calibration pattern: '{pattern}'")
return x, y, index, center_idx, fiducial_idx, q_rows, q_cols
[docs]
def calibration_mesh(x, y, index, idx_o, camera, pattern, dx, dy, dz, thickness,
camera_config, q_rows=None, q_cols=None, z_plane=0.0):
# =========================================================
# Algorithm to associate each point to a real position
# =========================================================
if q_rows is None or q_cols is None:
raise ValueError(
"calibration_mesh needs the absolute lattice coordinates "
"(q_rows/q_cols) produced by sort_points.")
def get_real_TSI(x, y, index, idx_o, camera):
# Absolute lattice coordinates from the predictive sort: q_cols is
# the column index, q_rows the row index (down-image positive), both
# zero at the origin dot. The TSI dual-plane levels form a
# checkerboard, so the Z bit is the (row + col) parity.
real_x = (q_cols - q_cols[idx_o]).astype(float)
real_y = (q_rows - q_rows[idx_o]).astype(float)
real_z = np.abs((q_rows - q_rows[idx_o]) +
(q_cols - q_cols[idx_o])).astype(float) % 2
# Apply scale
# real_y is positive going down in image-convention y (small y_img =
# top of image = physically high). Negate so that physically-higher
# points get positive world Y.
# real_z is the dual-plane bit {0, 1}; *dz -> planes at {0, dz}.
real_x *= dx
real_y *= -dy
real_z *= dz
# Front-Back: cameras face opposite sides (left = -right). thickness is
# split half each side of z=0; dz extends outward into the target depth.
# The mirrored (left) camera has BOTH its in-plane X and its near/far
# dual-plane ordering flipped, so the dual-plane bit (real_z) enters with
# a - sign on the left and a + sign on the right. Applying the bit flip
# to the wrong camera makes the two sub-levels land on the wrong side and
# blows up the Soloff reprojection RMS (~14 px vs ~1.6 px when correct).
if camera_config == "Front-Back":
if camera == "left":
real_z = -real_z - thickness / 2
real_x = -real_x
elif camera == "right":
real_z = thickness / 2 + real_z
else:
# Front-Front: +z toward cameras, so dz goes away (reference at
# z_plane, other dot at z_plane - dz).
real_z = -real_z
return real_x, real_y, real_z
def get_real_LaVision(x, y, index, idx_o, camera):
"""Assign world (X, Y, Z) to the staggered dual-plane dot lattice.
Absolute lattice coordinates come from the predictive sort: q_rows
is the image-row index (down-image positive, adjacent rows one Z
level apart), q_cols the column index in HALF-pitch units (the two
levels interleave with a half-pitch X stagger), both zero at the
origin dot. Rows with the origin's parity carry the Z=0 (reference)
level; the alternate rows carry the Z=dz level.
NB units: dx is the in-row dot pitch (one level), while dy is the
spacing between ADJACENT image rows -- half the per-level Y pitch,
because the two levels interleave. E.g. the LaVision Type #21 plate
(15 mm per-level pitch) needs dx = 15, dy = 7.5.
"""
real_x = (q_cols - q_cols[idx_o]).astype(float) / 2.0
real_y = (q_rows - q_rows[idx_o]).astype(float)
real_z = (np.abs(q_rows - q_rows[idx_o]) % 2).astype(float)
# Apply scale (negate real_y: image-row index increases downward, so
# rows above the origin must get positive world Y).
real_x *= dx
real_y *= -dy
real_z *= dz
if camera_config == "Front-Back":
raise Exception("Front-Back configuration is not available yet for LaVision target, as I have no example to test")
# NB: mirror the corrected TSI convention when enabling this — the
# left (mirrored) camera flips BOTH real_x and the dual-plane bit:
#if camera == "left":
# real_z = -real_z - thickness / 2
# real_x = - real_x
#elif camera == "right":
# real_z = thickness / 2 + real_z
else:
# Front-Front: +z toward cameras, so dz goes away (reference at
# z_plane, other dot at z_plane - dz).
real_z = -real_z
return real_x, real_y, real_z
def get_real_Quadricle(x, y, index, idx_o, camera):
# Use the absolute (row, col) grid coordinates from BFS so that
# partial and asymmetric grid visibility is handled correctly.
# NOTE: with a flat target (Z = 0 for all points) the out-of-plane
# polynomial coefficient a6 is unconstrained. For better stereo
# accuracy image the target at two depths and merge the CPT files.
# q_rows increases downward (image-convention), so negate to get
# positive world Y above the origin.
real_x = (q_cols - q_cols[idx_o]).astype(float) * dx
real_y = (q_rows[idx_o] - q_rows).astype(float) * dy
real_z = np.zeros(len(x))
if camera_config == "Front-Back":
if camera == "left":
real_x = -real_x
return real_x, real_y, real_z
if pattern == "TSI":
real_x, real_y, real_z = get_real_TSI(x, y, index, idx_o, camera)
elif pattern == "LaVision":
real_x, real_y, real_z = get_real_LaVision(x, y, index, idx_o, camera)
elif pattern == "Quadricle":
real_x, real_y, real_z = get_real_Quadricle(x, y, index, idx_o, camera)
# Shift to this capture's absolute target depth so CPTs taken at different
# Z planes can be combined. z_plane is the reference (closer) plane's Z; the
# second dual-plane level sits at z_plane - dz. For the flat Quadricle this
# is the single plane depth -- capturing it at >=2 distinct z_plane values is
# what makes Soloff possible for that target.
real_z = real_z + z_plane
return real_x, real_y, real_z
[docs]
class StereoCalibrationApp(tk.Tk):
def __init__(self, show_diagnostics=False):
super().__init__()
self.title("Stereo Calibration")
self.geometry("800x500")
# Developer flag: when True, sort_points pops its diagnostic plots
# (TSI column-separation / Quadricle BFS grid assignment). Off by
# default so the normal workflow is not interrupted by blocking plots;
# enable via calibration_GUI(show_diagnostics=True).
self.show_diagnostics = show_diagnostics
self.left_image = None
self.right_image = None
self.left_image_label = None
self.right_image_label = None
# Initialize grid options
self.camera_config = "Front-Back"
self.pattern = "TSI"
self.dx = 10
self.dy = 10
self.dz = 1
self.thickness = 4.25 # gap between the two cameras' near planes (mm),
# Front-Back only
self.line_len = 80 # Quadricle: morphological line-detection kernel (px)
self.z_plane = 0 # absolute Z (mm) of this capture's target plane,
# added to every point's Z so CPTs from different
# depths can be combined into a multi-plane fit
# Initialize processing options
self.threshold = 10
self.min_size = 9
self.block_size = 51 # Quadricle: adaptive threshold neighbourhood (px, odd)
# PIL images currently displayed (before resize), for dynamic re-rendering
self._left_pil = None
self._right_pil = None
self._resize_job = None
# Soloff calibration model: 'reduced' (9 params, linear in Z, needs
# >=2 Z planes) or 'full' (adds the Z^2 curvature, needs >=3 Z planes).
self.soloff_mode = 'reduced'
# Popup window references — prevent duplicates
self._grid_menu = None
self._options_menu = None
self._soloff_menu = None
self.create_widgets()
# Variables to store the coordinates of the box
self.box_start = None
self.box_end = None
self.black_box = 0
self.bind('<Configure>', self._on_resize)
# Bind mouse events to the labels
self.left_image_label.bind("<Button-1>", self.start_box_l)
self.left_image_label.bind("<ButtonRelease-1>", self.end_box)
self.left_image_label.bind("<B1-Motion>", self._draw_box_preview)
self.right_image_label.bind("<Button-1>", self.start_box_r)
self.right_image_label.bind("<ButtonRelease-1>", self.end_box)
self.right_image_label.bind("<B1-Motion>", self._draw_box_preview)
def _adaptive_threshold(self, raw):
"""Adaptive threshold for Quadricle: handles uneven illumination.
Returns a binary image with lines as foreground (255)."""
bs = self.block_size
if bs % 2 == 0:
bs += 1 # blockSize must be odd
return cv2.adaptiveThreshold(
raw, 255,
cv2.ADAPTIVE_THRESH_GAUSSIAN_C,
cv2.THRESH_BINARY_INV,
bs, 3)
def _process_image(self, raw):
"""Return the display/processing image from a raw grayscale image.
For Quadricle use adaptive threshold so the user sees the same binary
the detection step works on (lines white, background black)."""
if self.pattern == "Quadricle":
return self._adaptive_threshold(raw)
img = self.apply_threshold(raw)
img = self.bwareaopen(img)
return img
[docs]
def load_images(self):
"""
Load the right and left images of the Stereo Calibration and display
them
"""
left_image_path = filedialog.askopenfilename(
title="Select Left Image")
if left_image_path:
self.left_image_o = read_grayscale_image(left_image_path)
self.left_image = self._process_image(self.left_image_o)
self._left_pil = Image.fromarray(self.left_image)
self.left_image_plot = self.resize_image(self._left_pil)
self.left_image_label.config(image=self.left_image_plot)
right_image_path = filedialog.askopenfilename(
title="Select Right Image")
if right_image_path:
self.right_image_o = read_grayscale_image(right_image_path)
self.right_image = self._process_image(self.right_image_o)
self._right_pil = Image.fromarray(self.right_image)
self.right_image_plot = self.resize_image(self._right_pil)
self.right_image_label.config(image=self.right_image_plot)
[docs]
def resize_image(self, image):
# Resize the images to fit in the window
self.org_width, self.org_height = image.size
width = max(int(self.winfo_width() / 2.5), 1)
height = int((width / self.org_width) * self.org_height)
resized_image = image.resize((width, height), Image.Resampling.LANCZOS)
photo_image = ImageTk.PhotoImage(resized_image)
return photo_image
def _on_resize(self, event):
if event.widget is not self:
return
if self._resize_job:
self.after_cancel(self._resize_job)
self._resize_job = self.after(100, self._update_images)
def _update_images(self):
if self._left_pil is not None:
self.left_image_plot = self.resize_image(self._left_pil)
self.left_image_label.config(image=self.left_image_plot)
if self._right_pil is not None:
self.right_image_plot = self.resize_image(self._right_pil)
self.right_image_label.config(image=self.right_image_plot)
def _set_busy(self, message="Processing..."):
self.config(cursor="watch")
self.results_label.config(text=message)
self.update_idletasks()
def _set_idle(self):
self.config(cursor="")
[docs]
def reload_images(self):
# reload images
if hasattr(self, 'left_image_o'):
self.left_image = self._process_image(self.left_image_o)
self._left_pil = Image.fromarray(self.left_image)
self.left_image_plot = self.resize_image(self._left_pil)
self.left_image_label.config(image=self.left_image_plot)
if hasattr(self, 'right_image_o'):
self.right_image = self._process_image(self.right_image_o)
self._right_pil = Image.fromarray(self.right_image)
self.right_image_plot = self.resize_image(self._right_pil)
self.right_image_label.config(image=self.right_image_plot)
[docs]
def load_cpt_file(self, file_path):
x = []
y = []
X = []
Y = []
Z = []
with open(file_path, 'r') as file:
for line in file:
values = line.strip().split()
x.append(float(values[0]))
y.append(float(values[1]))
X.append(float(values[2]))
Y.append(float(values[3]))
Z.append(float(values[4]))
final_matrix = np.column_stack([x, y, X, Y, Z])
return final_matrix
[docs]
def save_cpt_file(self, filename, final_matrix):
with open(filename, 'w') as fid:
for row in final_matrix:
fid.write(f'{row[0]:e} {row[1]:e} '
f'{row[2]:e} {row[3]:e} '
f'{row[4]:e}\n')
[docs]
def apply_threshold(self, image):
threshold = np.max(image) * (self.threshold / 100)
dummy, binary_image = cv2.threshold(
image, threshold, 255, cv2.THRESH_BINARY
)
return binary_image
[docs]
def bwareaopen(self, img):
"""Remove small objects from binary image (approximation of
bwareaopen in Matlab for 2D images).
Parameters
----------
img : ndarray
Binary image (dtype=uint8) to remove small objects from.
Returns
-------
ndarray
The binary image with small objects removed.
"""
min_size = int(self.min_size)
kernel = cv2.getStructuringElement(cv2.MORPH_RECT, (min_size, min_size))
img = cv2.morphologyEx(img, cv2.MORPH_OPEN, kernel)
return img
[docs]
def start_blackbox(self):
if hasattr(self, 'left_image_o') and hasattr(self, 'right_image_o'):
self.black_box = 1
# Change the cursor to crosshair when the button is clicked
self.config(cursor="crosshair")
else:
messagebox.showerror("Error", "Load the images first!")
[docs]
def start_box_l(self, event):
# Store the starting coordinates of the box
self.box_start = [event.x, event.y]
self._drag_start = (event.x, event.y)
self.right = 0
self.left = 1
[docs]
def start_box_r(self, event):
# Store the starting coordinates of the box
self.box_start = [event.x, event.y]
self._drag_start = (event.x, event.y)
self.left = 0
self.right = 1
def _draw_box_preview(self, event):
if not self.black_box or not hasattr(self, '_drag_start'):
return
pil = self._left_pil if self.left else self._right_pil
label = self.left_image_label if self.left else self.right_image_label
if pil is None:
return
disp_w = int(self.winfo_width() / 2.5)
disp_h = int((disp_w / pil.width) * pil.height)
preview = pil.convert('RGB').resize(
(disp_w, disp_h), Image.Resampling.LANCZOS
)
draw = ImageDraw.Draw(preview)
x0, y0 = self._drag_start
x1, y1 = event.x, event.y
draw.rectangle(
[min(x0, x1), min(y0, y1), max(x0, x1), max(y0, y1)],
outline='lime', width=2
)
photo = ImageTk.PhotoImage(preview)
label.config(image=photo)
label._preview = photo # prevent garbage collection
[docs]
def end_box(self, event):
if self.black_box:
# Store the ending coordinates of the box
self.box_end = [event.x, event.y]
factor = self.org_width / self.left_image_plot.width()
self.box_end[0] = int(self.box_end[0] * factor)
self.box_end[1] = int(self.box_end[1] * factor)
self.box_start[0] = int(self.box_start[0] * factor)
self.box_start[1] = int(self.box_start[1] * factor)
self.box_end = tuple(self.box_end)
self.box_start = tuple(self.box_start)
# Draw the blackbox — white for Quadricle (paper colour masks lines)
color = 255 if self.pattern == "Quadricle" else 0
if self.left:
self.left_image_o = cv2.rectangle(
self.left_image_o, self.box_start,
self.box_end, color, -1)
elif self.right:
self.right_image_o = cv2.rectangle(
self.right_image_o, self.box_start,
self.box_end, color, -1)
self.config(cursor="")
self.reload_images()
self.black_box = 0
# Print the coordinates
print("Box coordinates: {}, {}".format(self.box_start, self.box_end))
[docs]
def get_centroids(self):
if hasattr(self, 'left_image_o') and hasattr(self, 'right_image_o'):
self._set_busy("Finding centroids...")
try:
# Get the centroid of each target dot
self.get_coordinates()
# Update point counts under each image
n_l = len(self.cX1)
n_r = len(self.cX2)
self.left_count_label.config(text=f"{n_l} points detected")
self.right_count_label.config(text=f"{n_r} points detected")
# Rearrange the centroids
x_l, y_l, index_l, center_idx_l, fiducial_idx_l = (
self.sort_points(self.cX1, self.cY1, self.center_l,
self.fiducial_l, camera="left",
showImages=self.show_diagnostics)
)
self._q_rows_l = self._q_rows_tmp
self._q_cols_l = self._q_cols_tmp
if self.pattern == "Quadricle":
self.origin_l_px = (int(x_l[center_idx_l]),
int(y_l[center_idx_l]))
x_r, y_r, index_r, center_idx_r, fiducial_idx_r = (
self.sort_points(self.cX2, self.cY2, self.center_r,
self.fiducial_r, camera="right",
showImages=self.show_diagnostics)
)
self._q_rows_r = self._q_rows_tmp
self._q_cols_r = self._q_cols_tmp
if self.pattern == "Quadricle":
self.origin_r_px = (int(x_r[center_idx_r]),
int(y_r[center_idx_r]))
# Draw detected points — for Quadricle this also overlays the
# cyan origin circle, so it runs after origin_*_px are set.
self.draw_circles()
# Create a mesh with real positions
X_real_l, Y_real_l, Z_real_l = self.calibrationMesh(
x_l, y_l, index_l, center_idx_l, "left")
X_real_r, Y_real_r, Z_real_r = self.calibrationMesh(
x_r, y_r, index_r, center_idx_r, "right")
# Create matrix with all the values
self.final_matrix_l = np.column_stack(
[x_l, y_l, X_real_l, Y_real_l, Z_real_l])
self.final_matrix_r = np.column_stack(
[x_r, y_r, X_real_r, Y_real_r, Z_real_r])
self.results_label.config(
text=f" Centroids found — Left: {n_l} pts |"
f" Right: {n_r} pts | Ready to run calculation."
)
except Exception as e:
messagebox.showerror("Error", f"Failed to find centroids:\n{e}")
self.results_label.config(text="Centroid detection failed.")
finally:
self._set_idle()
else:
messagebox.showerror("Error", "Load the images first!")
[docs]
def calibrationMesh(self, x, y, index, idx_o, camera):
q_rows = getattr(self, '_q_rows_l' if camera == "left" else '_q_rows_r', None)
q_cols = getattr(self, '_q_cols_l' if camera == "left" else '_q_cols_r', None)
return calibration_mesh(
x, y, index, idx_o, camera, self.pattern,
self.dx, self.dy, self.dz, self.thickness, self.camera_config,
q_rows, q_cols, z_plane=self.z_plane)
[docs]
def sort_points(self, x, y, center, fiducial, showImages=False, camera=None):
x, y, index, center_idx, fiducial_idx, q_rows, q_cols = sort_points(
x, y, center, fiducial, self.pattern, self.org_height,
showImages=showImages, camera=camera,
camera_config=self.camera_config)
self._q_rows_tmp = q_rows
self._q_cols_tmp = q_cols
return x, y, index, center_idx, fiducial_idx
[docs]
def get_coordinates(self):
def find_points_quadricle(orig_image):
"""
Detect grid-line intersections and the origin hole.
Works directly on the original grayscale image so it is not
affected by the apply_threshold / bwareaopen pipeline (which is
designed for bright blobs on dark backgrounds).
Strategy:
1. THRESH_BINARY_INV → black lines become white foreground.
2. Morphological OPEN with a long horizontal kernel isolates
horizontal lines; same with a vertical kernel for vertical
lines. ANDing the two gives only the intersection pixels.
3. Dilate + connectedComponents → sub-pixel centroids.
4. HoughCircles on the blurred original → darkest circle = hole.
"""
# --- Intersections ---
lines_fg = self._adaptive_threshold(orig_image)
ll = max(self.line_len, 5)
h_ker = cv2.getStructuringElement(cv2.MORPH_RECT, (ll, 1))
v_ker = cv2.getStructuringElement(cv2.MORPH_RECT, (1, ll))
h_lines = cv2.morphologyEx(lines_fg, cv2.MORPH_OPEN, h_ker)
v_lines = cv2.morphologyEx(lines_fg, cv2.MORPH_OPEN, v_ker)
inter = cv2.bitwise_and(h_lines, v_lines)
dil = cv2.getStructuringElement(cv2.MORPH_ELLIPSE, (11, 11))
inter = cv2.dilate(inter, dil)
_, inter = cv2.threshold(inter, 127, 255, cv2.THRESH_BINARY)
num_labels, _, _, centroids = cv2.connectedComponentsWithStats(inter)
# Keep sub-pixel float centroids (see find_points_TSI).
cX = [float(centroids[i][0]) for i in range(1, num_labels)]
cY = [float(centroids[i][1]) for i in range(1, num_labels)]
# --- Origin hole: darkest circular region ---
blurred = cv2.GaussianBlur(orig_image, (7, 7), 0)
h_img, w_img = orig_image.shape
min_r = max(5, w_img // 200)
max_r = max(50, w_img // 40)
circles = cv2.HoughCircles(
blurred, cv2.HOUGH_GRADIENT, dp=1,
minDist=min_r * 4,
param1=50, param2=20,
minRadius=min_r, maxRadius=max_r)
if circles is not None:
circles = np.round(circles[0]).astype(int)
best = min(circles,
key=lambda c: int(orig_image[c[1], c[0]]))
center = [int(best[0]), int(best[1])]
else:
# Fallback: centroid of the darkest 1 % of pixels
p1 = int(np.percentile(orig_image, 1))
_, dark = cv2.threshold(
orig_image, p1, 255, cv2.THRESH_BINARY_INV)
nl, _, _, cents = cv2.connectedComponentsWithStats(dark)
center = ([int(cents[1][0]), int(cents[1][1])]
if nl > 1 else [cX[0], cY[0]])
return cX, cY, center, None # no fiducial for this target
if self.pattern == "TSI":
# Get coordinates of left image
cX1, cY1, center_l, fiducial_l = find_points_TSI(self.left_image)
# Get coordinates of right image
cX2, cY2, center_r, fiducial_r = find_points_TSI(self.right_image)
elif self.pattern == "LaVision":
# Both fiducials are returned together and excluded from the dot
# list; 'center' is the origin dot (next to the square, on top of
# the triangle). The [square, triangle] pair is passed on to
# sort_points, which uses square -> origin as the +X direction.
cX1, cY1, center_l, self.fiducials_l = find_points_LaVision(
self.left_image)
cX2, cY2, center_r, self.fiducials_r = find_points_LaVision(
self.right_image)
fiducial_l = self.fiducials_l
fiducial_r = self.fiducials_r
elif self.pattern == "Quadricle":
cX1, cY1, center_l, fiducial_l = find_points_quadricle(
self.left_image_o)
cX2, cY2, center_r, fiducial_r = find_points_quadricle(
self.right_image_o)
# Get the variables in the class
self.cX1 = cX1
self.cX2 = cX2
self.center_l = center_l
self.center_r = center_r
self.cY1 = np.asarray(cY1)
self.cY2 = np.asarray(cY2)
self.fiducial_l = fiducial_l
self.fiducial_r = fiducial_r
[docs]
def draw_circles(self):
thickness = 2
markerType = cv2.MARKER_CROSS
if self.pattern == "Quadricle":
# Draw on the original image; use smaller markers to avoid clutter
base_l = cv2.cvtColor(self.left_image_o, cv2.COLOR_GRAY2RGB)
base_r = cv2.cvtColor(self.right_image_o, cv2.COLOR_GRAY2RGB)
radius = 8
markerSize = 12
else:
base_l = cv2.cvtColor(self.left_image, cv2.COLOR_GRAY2RGB)
base_r = cv2.cvtColor(self.right_image, cv2.COLOR_GRAY2RGB)
radius = 40
markerSize = 15
color = (255, 0, 0)
cyan = (255, 255, 0) # BGR cyan
# Images here are RGB (GRAY2RGB + PIL), so (255, 255, 0) renders yellow.
yellow = (255, 255, 0)
left_image = base_l
for i in range(len(self.cX1)):
# cv2 needs integer pixel coordinates; the stored centroids are
# sub-pixel floats, so round only here for drawing.
center = (int(round(self.cX1[i])), int(round(self.cY1[i])))
left_image = cv2.circle(
left_image, center, radius, color, thickness)
cv2.drawMarker(
left_image, center, color, markerType, markerSize, thickness)
cv2.drawMarker(
left_image, (int(round(self.center_l[0])), int(round(self.center_l[1]))),
(0, 255, 0), markerType, 100, thickness)
# Mark the two LaVision fiducials (square, triangle) in yellow.
if self.pattern == "LaVision" and getattr(self, 'fiducials_l', None):
for fx, fy in self.fiducials_l:
cv2.circle(left_image, (int(fx), int(fy)), radius, yellow,
thickness)
cv2.drawMarker(left_image, (int(fx), int(fy)), yellow,
cv2.MARKER_TILTED_CROSS, markerSize, thickness)
if self.pattern == "Quadricle" and hasattr(self, 'origin_l_px'):
cv2.circle(left_image, self.origin_l_px, radius * 2, cyan, thickness + 1)
self._left_pil = Image.fromarray(left_image)
self.left_image_plot = self.resize_image(self._left_pil)
self.left_image_label.config(image=self.left_image_plot)
right_image = base_r
for i in range(len(self.cX2)):
# cv2 needs integer pixel coordinates; round the sub-pixel
# float centroids only here for drawing.
center = (int(round(self.cX2[i])), int(round(self.cY2[i])))
right_image = cv2.circle(
right_image, center, radius, color, thickness)
cv2.drawMarker(
right_image, center, color, markerType, markerSize, thickness)
cv2.drawMarker(
right_image, (int(round(self.center_r[0])), int(round(self.center_r[1]))),
(0, 255, 0), markerType, 100, thickness)
# Mark the two LaVision fiducials (square, triangle) in yellow.
if self.pattern == "LaVision" and getattr(self, 'fiducials_r', None):
for fx, fy in self.fiducials_r:
cv2.circle(right_image, (int(fx), int(fy)), radius, yellow,
thickness)
cv2.drawMarker(right_image, (int(fx), int(fy)), yellow,
cv2.MARKER_TILTED_CROSS, markerSize, thickness)
if self.pattern == "Quadricle" and hasattr(self, 'origin_r_px'):
cv2.circle(right_image, self.origin_r_px, radius * 2, cyan, thickness + 1)
self._right_pil = Image.fromarray(right_image)
self.right_image_plot = self.resize_image(self._right_pil)
self.right_image_label.config(image=self.right_image_plot)
[docs]
def run_calculation(self):
if not hasattr(self, 'final_matrix_l') or not hasattr(self, 'final_matrix_r'):
messagebox.showerror("Error", "Process Images or load CPT files")
else:
self._set_busy("Running calibration optimisation — please wait...")
try:
self.stereo_calibration()
except Exception as e:
messagebox.showerror("Error", f"Calibration failed:\n{e}")
self.results_label.config(text="Calibration failed.")
finally:
self._set_idle()
[docs]
def stereo_calibration(self):
(fwd_a_l, fwd_b_l, fwd_a_r, fwd_b_r,
map_a_l, map_b_l, map_a_r, map_b_r,
results_text) = stereo_calibration(
self.final_matrix_l, self.final_matrix_r,
self.org_width, self.org_height, mode=self.soloff_mode)
self.fwd_a_l = fwd_a_l
self.fwd_b_l = fwd_b_l
self.fwd_a_r = fwd_a_r
self.fwd_b_r = fwd_b_r
self.map_a_l = map_a_l
self.map_b_l = map_b_l
self.map_a_r = map_a_r
self.map_b_r = map_b_r
if self.camera_config == 'Front-Back':
sign_note = (
"\n\n Config: Front-Back"
"\n u_z sign convention: negative = toward left camera, positive = toward right camera"
)
else:
sign_note = (
"\n\n Config: Front-Front"
"\n u_z sign convention: positive = flow toward both cameras (out of light sheet)"
)
self.results_label.config(text=results_text + sign_note)
[docs]
def save_calibration(self):
if hasattr(self, 'fwd_a_l') and hasattr(self, 'map_a_l'):
file_path = filedialog.asksaveasfilename(
title="Calibration file", defaultextension=".npz",
filetypes=[("numpy files", "*.npz")]
)
if not file_path:
return
np.savez(file_path,
fwd_a_l=self.fwd_a_l, fwd_b_l=self.fwd_b_l,
fwd_a_r=self.fwd_a_r, fwd_b_r=self.fwd_b_r,
map_a_l=self.map_a_l, map_b_l=self.map_b_l,
map_a_r=self.map_a_r, map_b_r=self.map_b_r,
camera_config=self.camera_config)
else:
messagebox.showerror("Error", "Run calculation first!")
[docs]
def calibration_GUI(show_diagnostics=False):
app = StereoCalibrationApp(show_diagnostics=show_diagnostics)
app.mainloop()
if __name__ == "__main__":
calibration_GUI()