"""ZIRULAB technical demo: run the Brinell measurement code on SYNTHETIC images.

This checks that the code does what the demo page says (pixel -> mm, ellipse fit, HBW, range check).
It is not a measurement of real indentations and says nothing about real-world accuracy
(real indentations have raised rims, scratches, uneven lighting and lens distortion).

usage: python3 brinell_synthetic_check.py   (needs Python 3, numpy, opencv-python)
"""
import math
import cv2
import numpy as np

# ---- measurement code ----
def measure_indentation(path_or_img, mm_per_px):
    """Return (d1, d2, d) in mm. Assumes the indentation is darker than its surroundings."""
    img = path_or_img if isinstance(path_or_img, np.ndarray) else cv2.imread(path_or_img, cv2.IMREAD_GRAYSCALE)
    if img is None:
        raise FileNotFoundError(path_or_img)
    blur = cv2.GaussianBlur(img, (5, 5), 0)
    _, mask = cv2.threshold(blur, 0, 255, cv2.THRESH_BINARY_INV + cv2.THRESH_OTSU)
    mask = cv2.morphologyEx(mask, cv2.MORPH_OPEN, np.ones((5, 5), np.uint8))
    contours, _ = cv2.findContours(mask, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_NONE)
    if not contours:
        raise ValueError("no indentation contour found")
    c = max(contours, key=cv2.contourArea)
    if len(c) < 5:
        raise ValueError("contour too small to fit an ellipse")
    (_, _), (axis1, axis2), _ = cv2.fitEllipse(c)  # full axis lengths in pixels
    problems = plausibility(c, mask.shape, axis1, axis2)
    if problems:
        raise ValueError("segmentation not plausible: " + "; ".join(problems))
    d1, d2 = axis1 * mm_per_px, axis2 * mm_per_px
    return d1, d2, (d1 + d2) / 2


def plausibility(contour, shape, axis1, axis2, max_axis_diff=0.05, min_fill=0.97):
    """Reasons to distrust a segmentation instead of silently returning a number."""
    problems = []
    x, y, w, h = cv2.boundingRect(contour)
    if x <= 1 or y <= 1 or x + w >= shape[1] - 1 or y + h >= shape[0] - 1:
        problems.append("contour touches the image border")
    if abs(axis1 - axis2) / max(axis1, axis2) > max_axis_diff:
        problems.append(f"axes differ by more than {max_axis_diff:.0%}")
    fill = cv2.contourArea(contour) / (np.pi * axis1 * axis2 / 4)
    if fill < min_fill:
        problems.append(f"contour fills only {fill:.0%} of the fitted ellipse")
    return problems


def brinell_hbw(force_n, ball_d_mm, d_mm):
    """HBW per ISO 6506-1. Returns (hbw, in_range); in_range is False if d is outside 0.24D-0.6D."""
    D, d = ball_d_mm, d_mm
    hbw = 0.102 * 2 * force_n / (np.pi * D * (D - np.sqrt(D**2 - d**2)))
    return hbw, 0.24 * D <= d <= 0.6 * D
# ----------------------------------------------------------------------------------------------


def synthetic(d1_px, d2_px, size=1200, noise=6.0, seed=0):
    """Light background, dark elliptical indentation with a blurred edge, Gaussian noise."""
    rng = np.random.default_rng(seed)
    img = np.full((size, size), 180.0)
    cv2.ellipse(img, (size // 2, size // 2), (int(round(d1_px / 2)), int(round(d2_px / 2))), 0, 0, 360, 70.0, -1, lineType=cv2.LINE_AA)
    img = cv2.GaussianBlur(img, (0, 0), 2.0) + rng.normal(0, noise, img.shape)
    return np.clip(img, 0, 255).astype(np.uint8)


if __name__ == "__main__":
    MM_PER_PX = 0.005  # assumed calibration: 5 um per pixel
    print(f"opencv {cv2.__version__}; assumed calibration {MM_PER_PX} mm/px; F = 29420 N, D = 10 mm")
    print("true d1/d2 (px) | measured d1/d2 (px) | d true -> measured (mm) | HBW true -> measured | in range")
    for d1, d2 in [(800, 800), (800, 780), (500, 500), (1000, 990), (400, 400)]:
        m1, m2, d = measure_indentation(synthetic(d1, d2), MM_PER_PX)
        true_d = (d1 + d2) / 2 * MM_PER_PX
        ht, _ = brinell_hbw(29420, 10, true_d)
        hm, ok = brinell_hbw(29420, 10, d)
        a, b = sorted((m1 / MM_PER_PX, m2 / MM_PER_PX), reverse=True)
        print(f"{d1}/{d2} | {a:.1f}/{b:.1f} | {true_d:.4f} -> {d:.4f} ({(d / true_d - 1) * 100:+.2f}%) | {ht:.1f} -> {hm:.1f} ({(hm / ht - 1) * 100:+.2f}%) | {ok}")

    print("\nuneven illumination (linear brightness gradient across the image, grey levels):")
    for g in (0, 40, 80, 120):
        img = synthetic(800, 800).astype(float) + np.linspace(-g / 2, g / 2, 1200)[None, :]
        img = np.clip(img, 0, 255).astype(np.uint8)
        try:
            _, _, d = measure_indentation(img, MM_PER_PX)
            print(f"  gradient {g:3d}: d = {d:.4f} mm ({(d / 4 - 1) * 100:+.2f}%)")
        except ValueError as e:
            print(f"  gradient {g:3d}: rejected ({e})")