package main

import (
	"fmt"
	"image"
	"image/color"
	_ "image/jpeg"
	_ "image/png"
	"math"
	"os"
	"sort"
)

type ScanProcessResult struct {
	Faces  []TracedFace
	DPI    float64
	Errors []string
}

// ProcessFaceScans processes scanned face template images and extracts traced outlines.
// Each imagePath corresponds to a page (1-indexed by order).
func ProcessFaceScans(imagePaths []string, cfg FaceTemplateConfig) (*ScanProcessResult, error) {
	if len(imagePaths) == 0 {
		return nil, fmt.Errorf("no images provided")
	}

	if cfg.PageWidth <= 0 {
		cfg.PageWidth = 215.9
	}
	if cfg.PageHeight <= 0 {
		cfg.PageHeight = 279.4
	}
	if cfg.NumFaces < 1 {
		cfg.NumFaces = 6
	}
	if cfg.LongestSide <= 0 {
		cfg.LongestSide = 50
	}
	if cfg.GridSpacing <= 0 {
		cfg.GridSpacing = pickGridSpacing(cfg.LongestSide)
	}

	layout := computePageLayout(cfg)

	result := &ScanProcessResult{}

	for i, imgPath := range imagePaths {
		pageNum := i + 1
		if pageNum > len(layout) {
			result.Errors = append(result.Errors, fmt.Sprintf("page %d: no layout (only %d pages expected)", pageNum, len(layout)))
			continue
		}

		debugLog("ProcessFaceScans: page %d from %s", pageNum, imgPath)

		faces, dpi, err := processOnePage(imgPath, cfg, layout[i])
		if err != nil {
			result.Errors = append(result.Errors, fmt.Sprintf("page %d: %v", pageNum, err))
			continue
		}

		if result.DPI == 0 {
			result.DPI = dpi
		}
		result.Faces = append(result.Faces, faces...)
	}

	debugLog("ProcessFaceScans: extracted %d faces total", len(result.Faces))
	return result, nil
}

func processOnePage(imgPath string, cfg FaceTemplateConfig, page pageLayout) ([]TracedFace, float64, error) {
	img, err := loadScanImage(imgPath)
	if err != nil {
		return nil, 0, err
	}

	bounds := img.Bounds()
	w, h := bounds.Dx(), bounds.Dy()
	dpi := estimateDPI(w, h, cfg.PageWidth, cfg.PageHeight)
	debugLog("  image: %dx%d, estimated DPI: %.0f", w, h, dpi)

	// Phase 1: detect bullseye fiducial markers
	markers := DetectMarkers(img, dpi)
	debugLog("  detected %d markers", len(markers))

	expectedPos := markerPositionsMM(cfg.PageWidth, cfg.PageHeight)

	// Build alignment from detected markers → expected mm positions
	var srcPts, dstPts [4][2]float64
	nAlignPts := 0

	// Match detected markers to expected positions by corner ID
	for _, m := range markers {
		if m.Data.CornerID >= 0 && m.Data.CornerID <= 3 && nAlignPts < 4 {
			srcPts[nAlignPts] = m.PixelCenter
			dstPts[nAlignPts] = expectedPos[m.Data.CornerID]
			nAlignPts++
		}
	}

	// Fallback: if insufficient markers detected, use the old registration mark approach
	if nAlignPts < 3 {
		debugLog("  WARNING: only %d markers detected, falling back to registration circles", nAlignPts)
		margin := 15.0
		oldExpected := [4][2]float64{
			{margin, margin},
			{cfg.PageWidth - margin, margin},
			{margin, cfg.PageHeight - margin},
			{cfg.PageWidth - margin, cfg.PageHeight - margin},
		}
		detectedCorners, err := detectRegistrationMarks(img, dpi, oldExpected)
		if err != nil {
			return nil, dpi, fmt.Errorf("marker detection failed and registration fallback failed: %w", err)
		}
		srcPts = detectedCorners
		dstPts = oldExpected
		nAlignPts = 4
	}

	xform := computeAffine(srcPts, dstPts)
	debugLog("  affine: scale~%.4f, residual~%.2fpx", xform.scale(), xform.residual(srcPts, dstPts, dpi))

	// Phase 2: calibrate DPI from data bars (if markers were detected)
	if len(markers) >= 3 {
		barX := expectedPos[0][0] + float64(markerN)*markerCellMM/2 + 3
		barDPI := MeasureCalibBarDPI(img, otsuThreshold(img), xform, barX, cfg.PageHeight-15.0, calibBarSpecs())
		if barDPI > 100 {
			debugLog("  data bar calibrated DPI: %.1f (was %.1f)", barDPI, dpi)
			dpi = barDPI
		}
	}

	// Phase 3: extract face outlines using template-aware erasure
	isTemplate := TemplateElementMap(cfg, page)

	var faces []TracedFace
	for _, cell := range page.Cells {
		outline, err := extractFaceOutlineClean(img, xform, cell, dpi, isTemplate)
		if err != nil {
			debugLog("  face %d: %v", cell.FaceNum, err)
			continue
		}
		if len(outline) >= 3 {
			faces = append(faces, TracedFace{FaceNum: cell.FaceNum, Outline: outline})
			debugLog("  face %d: %d vertices", cell.FaceNum, len(outline))
		}
	}

	return faces, dpi, nil
}

// --- image loading ---

func loadScanImage(path string) (*image.Gray, error) {
	f, err := os.Open(path)
	if err != nil {
		return nil, err
	}
	defer f.Close()

	img, _, err := image.Decode(f)
	if err != nil {
		return nil, fmt.Errorf("decode %s: %w", path, err)
	}

	return toGrayscale(img), nil
}

func toGrayscale(img image.Image) *image.Gray {
	bounds := img.Bounds()
	gray := image.NewGray(bounds)
	for y := bounds.Min.Y; y < bounds.Max.Y; y++ {
		for x := bounds.Min.X; x < bounds.Max.X; x++ {
			r, g, b, _ := img.At(x, y).RGBA()
			lum := uint8((0.299*float64(r) + 0.587*float64(g) + 0.114*float64(b)) / 256)
			gray.SetGray(x, y, color.Gray{Y: lum})
		}
	}
	return gray
}

// --- DPI estimation ---

func estimateDPI(pixW, pixH int, pageWmm, pageHmm float64) float64 {
	dpiW := float64(pixW) / (pageWmm / 25.4)
	dpiH := float64(pixH) / (pageHmm / 25.4)
	return (dpiW + dpiH) / 2
}

func mmToPx(mm, dpi float64) float64 { return mm * dpi / 25.4 }
func pxToMM(px, dpi float64) float64 { return px * 25.4 / dpi }

// --- thresholding ---

func otsuThreshold(img *image.Gray) uint8 {
	bounds := img.Bounds()
	var hist [256]int
	for y := bounds.Min.Y; y < bounds.Max.Y; y++ {
		for x := bounds.Min.X; x < bounds.Max.X; x++ {
			hist[img.GrayAt(x, y).Y]++
		}
	}

	total := bounds.Dx() * bounds.Dy()
	var sumAll float64
	for i := 0; i < 256; i++ {
		sumAll += float64(i) * float64(hist[i])
	}

	var sumBg float64
	var wBg int
	maxVariance := 0.0
	bestT := uint8(0)

	for t := 0; t < 256; t++ {
		wBg += hist[t]
		if wBg == 0 {
			continue
		}
		wFg := total - wBg
		if wFg == 0 {
			break
		}

		sumBg += float64(t) * float64(hist[t])
		meanBg := sumBg / float64(wBg)
		meanFg := (sumAll - sumBg) / float64(wFg)

		variance := float64(wBg) * float64(wFg) * (meanBg - meanFg) * (meanBg - meanFg)
		if variance > maxVariance {
			maxVariance = variance
			bestT = uint8(t)
		}
	}
	return bestT
}

// adaptiveThreshold applies block-mean adaptive thresholding.
// Returns a binary mask: true = dark (ink), false = light (paper).
func adaptiveThreshold(img *image.Gray, blockSize, offset int) []bool {
	bounds := img.Bounds()
	w, h := bounds.Dx(), bounds.Dy()
	mask := make([]bool, w*h)

	// integral image for fast block mean
	integral := make([]int64, (w+1)*(h+1))
	stride := w + 1
	for y := 0; y < h; y++ {
		var rowSum int64
		for x := 0; x < w; x++ {
			rowSum += int64(img.GrayAt(bounds.Min.X+x, bounds.Min.Y+y).Y)
			integral[(y+1)*stride+(x+1)] = integral[y*stride+(x+1)] + rowSum
		}
	}

	half := blockSize / 2
	for y := 0; y < h; y++ {
		for x := 0; x < w; x++ {
			x0 := max(0, x-half)
			y0 := max(0, y-half)
			x1 := min(w, x+half+1)
			y1 := min(h, y+half+1)

			count := (x1 - x0) * (y1 - y0)
			sum := integral[y1*stride+x1] - integral[y0*stride+x1] - integral[y1*stride+x0] + integral[y0*stride+x0]
			mean := int(sum / int64(count))

			val := int(img.GrayAt(bounds.Min.X+x, bounds.Min.Y+y).Y)
			mask[y*w+x] = val < mean-offset
		}
	}
	return mask
}

// --- morphological operations ---

func circularKernel(radius int) []bool {
	size := 2*radius + 1
	k := make([]bool, size*size)
	r2 := float64(radius * radius)
	for y := 0; y < size; y++ {
		for x := 0; x < size; x++ {
			dy := float64(y - radius)
			dx := float64(x - radius)
			k[y*size+x] = dx*dx+dy*dy <= r2
		}
	}
	return k
}

func erode(mask []bool, w, h, radius int) []bool {
	kernel := circularKernel(radius)
	ksize := 2*radius + 1
	out := make([]bool, w*h)

	for y := radius; y < h-radius; y++ {
		for x := radius; x < w-radius; x++ {
			allSet := true
			for ky := 0; ky < ksize && allSet; ky++ {
				for kx := 0; kx < ksize && allSet; kx++ {
					if kernel[ky*ksize+kx] && !mask[(y-radius+ky)*w+(x-radius+kx)] {
						allSet = false
					}
				}
			}
			out[y*w+x] = allSet
		}
	}
	return out
}

func dilate(mask []bool, w, h, radius int) []bool {
	kernel := circularKernel(radius)
	ksize := 2*radius + 1
	out := make([]bool, w*h)

	for y := radius; y < h-radius; y++ {
		for x := radius; x < w-radius; x++ {
			if !mask[y*w+x] {
				continue
			}
			for ky := 0; ky < ksize; ky++ {
				for kx := 0; kx < ksize; kx++ {
					if kernel[ky*ksize+kx] {
						out[(y-radius+ky)*w+(x-radius+kx)] = true
					}
				}
			}
		}
	}
	return out
}

func morphOpen(mask []bool, w, h, radius int) []bool {
	return dilate(erode(mask, w, h, radius), w, h, radius)
}

// --- blob detection (connected components) ---

type blob struct {
	cx, cy float64
	area   int
}

func findBlobs(mask []bool, w, h int) []blob {
	labels := make([]int, w*h)
	nextLabel := 1
	blobPixels := map[int][][2]int{}

	for y := 0; y < h; y++ {
		for x := 0; x < w; x++ {
			if !mask[y*w+x] || labels[y*w+x] != 0 {
				continue
			}
			// BFS flood fill
			label := nextLabel
			nextLabel++
			queue := [][2]int{{x, y}}
			labels[y*w+x] = label
			var pixels [][2]int

			for len(queue) > 0 {
				p := queue[0]
				queue = queue[1:]
				pixels = append(pixels, p)

				for _, d := range [][2]int{{-1, 0}, {1, 0}, {0, -1}, {0, 1}} {
					nx, ny := p[0]+d[0], p[1]+d[1]
					if nx >= 0 && nx < w && ny >= 0 && ny < h && mask[ny*w+nx] && labels[ny*w+nx] == 0 {
						labels[ny*w+nx] = label
						queue = append(queue, [2]int{nx, ny})
					}
				}
			}
			blobPixels[label] = pixels
		}
	}

	var blobs []blob
	for _, pixels := range blobPixels {
		var sx, sy float64
		for _, p := range pixels {
			sx += float64(p[0])
			sy += float64(p[1])
		}
		n := float64(len(pixels))
		blobs = append(blobs, blob{cx: sx / n, cy: sy / n, area: len(pixels)})
	}
	return blobs
}

// --- registration mark detection ---

func detectRegistrationMarks(img *image.Gray, dpi float64, expectedMM [4][2]float64) ([4][2]float64, error) {
	threshold := otsuThreshold(img)
	bounds := img.Bounds()
	w, h := bounds.Dx(), bounds.Dy()

	// Expected circle area: pi * r^2 where r = 0.8mm
	circRadPx := mmToPx(0.8, dpi)
	expectedArea := math.Pi * circRadPx * circRadPx
	minArea := int(expectedArea * 0.3)
	maxArea := int(expectedArea * 3.0)

	searchRadPx := int(mmToPx(20, dpi)) // 20mm search window

	var detected [4][2]float64

	for ci, mmPt := range expectedMM {
		ex := int(mmToPx(mmPt[0], dpi))
		ey := int(mmToPx(mmPt[1], dpi))

		x0 := max(0, ex-searchRadPx)
		y0 := max(0, ey-searchRadPx)
		x1 := min(w, ex+searchRadPx)
		y1 := min(h, ey+searchRadPx)

		roiW, roiH := x1-x0, y1-y0
		roiMask := make([]bool, roiW*roiH)
		for ry := 0; ry < roiH; ry++ {
			for rx := 0; rx < roiW; rx++ {
				roiMask[ry*roiW+rx] = img.GrayAt(bounds.Min.X+x0+rx, bounds.Min.Y+y0+ry).Y < threshold
			}
		}

		blobs := findBlobs(roiMask, roiW, roiH)

		// Find blob closest to expected position with area in range
		bestDist := math.MaxFloat64
		bestIdx := -1
		for bi, b := range blobs {
			if b.area < minArea || b.area > maxArea {
				continue
			}
			// blob coords are relative to ROI
			absX := float64(x0) + b.cx
			absY := float64(y0) + b.cy
			dist := math.Hypot(absX-float64(ex), absY-float64(ey))
			if dist < bestDist {
				bestDist = dist
				bestIdx = bi
			}
		}

		if bestIdx < 0 {
			// Fallback: use expected position
			debugLog("  WARNING: registration mark %d not found, using expected position", ci)
			detected[ci] = [2]float64{float64(ex), float64(ey)}
		} else {
			b := blobs[bestIdx]
			detected[ci] = [2]float64{float64(x0) + b.cx, float64(y0) + b.cy}
			debugLog("  mark %d: detected at (%.1f,%.1f), expected (%d,%d), dist=%.1fpx, area=%d",
				ci, detected[ci][0], detected[ci][1], ex, ey, bestDist, b.area)
		}
	}

	return detected, nil
}

// --- affine transform ---

type affineTransform struct {
	a, b, c float64 // xOut = a*xIn + b*yIn + c
	d, e, f float64 // yOut = d*xIn + e*yIn + f
}

// computeAffine computes the least-squares affine transform mapping srcPts (pixels) to dstPts (mm).
func computeAffine(srcPts [4][2]float64, dstPts [4][2]float64) affineTransform {
	// Solve for [a,b,c] and [d,e,f] independently via normal equations.
	// For each: A * params = b where A is Nx3, params is 3x1, b is Nx1.
	// A = [[x0,y0,1],[x1,y1,1],...], bx = [dx0,dx1,...], by = [dy0,dy1,...]

	var ata [3][3]float64
	var atbx, atby [3]float64

	for i := 0; i < 4; i++ {
		sx, sy := srcPts[i][0], srcPts[i][1]
		dx, dy := dstPts[i][0], dstPts[i][1]
		row := [3]float64{sx, sy, 1}

		for r := 0; r < 3; r++ {
			for c := 0; c < 3; c++ {
				ata[r][c] += row[r] * row[c]
			}
			atbx[r] += row[r] * dx
			atby[r] += row[r] * dy
		}
	}

	inv := invert3x3(ata)
	var xform affineTransform
	for i := 0; i < 3; i++ {
		xform.a += inv[0][i] * atbx[i]
		xform.b += inv[1][i] * atbx[i]
		xform.c += inv[2][i] * atbx[i]
		xform.d += inv[0][i] * atby[i]
		xform.e += inv[1][i] * atby[i]
		xform.f += inv[2][i] * atby[i]
	}
	return xform
}

func invert3x3(m [3][3]float64) [3][3]float64 {
	det := m[0][0]*(m[1][1]*m[2][2]-m[1][2]*m[2][1]) -
		m[0][1]*(m[1][0]*m[2][2]-m[1][2]*m[2][0]) +
		m[0][2]*(m[1][0]*m[2][1]-m[1][1]*m[2][0])

	if math.Abs(det) < 1e-12 {
		return [3][3]float64{{1, 0, 0}, {0, 1, 0}, {0, 0, 1}}
	}

	invDet := 1.0 / det
	return [3][3]float64{
		{(m[1][1]*m[2][2] - m[1][2]*m[2][1]) * invDet, (m[0][2]*m[2][1] - m[0][1]*m[2][2]) * invDet, (m[0][1]*m[1][2] - m[0][2]*m[1][1]) * invDet},
		{(m[1][2]*m[2][0] - m[1][0]*m[2][2]) * invDet, (m[0][0]*m[2][2] - m[0][2]*m[2][0]) * invDet, (m[0][2]*m[1][0] - m[0][0]*m[1][2]) * invDet},
		{(m[1][0]*m[2][1] - m[1][1]*m[2][0]) * invDet, (m[0][1]*m[2][0] - m[0][0]*m[2][1]) * invDet, (m[0][0]*m[1][1] - m[0][1]*m[1][0]) * invDet},
	}
}

func (t affineTransform) transform(px, py float64) (float64, float64) {
	return t.a*px + t.b*py + t.c, t.d*px + t.e*py + t.f
}

func (t affineTransform) scale() float64 {
	return math.Sqrt(t.a*t.a + t.d*t.d)
}

func (t affineTransform) residual(srcPts [4][2]float64, dstPts [4][2]float64, dpi float64) float64 {
	var sum float64
	for i := 0; i < 4; i++ {
		ox, oy := t.transform(srcPts[i][0], srcPts[i][1])
		dx := ox - dstPts[i][0]
		dy := oy - dstPts[i][1]
		sum += mmToPx(math.Sqrt(dx*dx+dy*dy), dpi)
	}
	return sum / 4
}

// invertAffine returns the inverse transform (mm -> pixels).
func invertAffine(t affineTransform) affineTransform {
	det := t.a*t.e - t.b*t.d
	if math.Abs(det) < 1e-12 {
		return t
	}
	invDet := 1.0 / det
	return affineTransform{
		a: t.e * invDet,
		b: -t.b * invDet,
		c: (t.b*t.f - t.e*t.c) * invDet,
		d: -t.d * invDet,
		e: t.a * invDet,
		f: (t.d*t.c - t.a*t.f) * invDet,
	}
}

// --- face outline extraction (new: template-aware) ---

// extractFaceOutlineClean uses template element erasure instead of morphological hacks.
// It thresholds the cell area, erases all known template elements by position, then
// traces the remaining contour — which is the user's hand-drawn outline.
func extractFaceOutlineClean(img *image.Gray, xform affineTransform, cell faceCell, dpi float64, isTemplate func(float64, float64) bool) ([][2]float64, error) {
	inv := invertAffine(xform)

	pad := 2.0
	mmX0, mmY0 := cell.X-pad, cell.Y-pad
	mmX1, mmY1 := cell.X+cell.W+pad, cell.Y+cell.H+pad

	corners := [4][2]float64{
		{mmX0, mmY0}, {mmX1, mmY0}, {mmX0, mmY1}, {mmX1, mmY1},
	}
	var pxMinX, pxMinY float64 = math.MaxFloat64, math.MaxFloat64
	var pxMaxX, pxMaxY float64 = -math.MaxFloat64, -math.MaxFloat64
	for _, c := range corners {
		px, py := inv.transform(c[0], c[1])
		pxMinX = math.Min(pxMinX, px)
		pxMinY = math.Min(pxMinY, py)
		pxMaxX = math.Max(pxMaxX, px)
		pxMaxY = math.Max(pxMaxY, py)
	}

	bounds := img.Bounds()
	cropX0 := max(bounds.Min.X, int(pxMinX))
	cropY0 := max(bounds.Min.Y, int(pxMinY))
	cropX1 := min(bounds.Max.X, int(pxMaxX)+1)
	cropY1 := min(bounds.Max.Y, int(pxMaxY)+1)
	cropW := cropX1 - cropX0
	cropH := cropY1 - cropY0

	if cropW < 10 || cropH < 10 {
		return nil, fmt.Errorf("cell crop too small: %dx%d", cropW, cropH)
	}

	sub := image.NewGray(image.Rect(0, 0, cropW, cropH))
	for y := 0; y < cropH; y++ {
		for x := 0; x < cropW; x++ {
			sub.SetGray(x, y, img.GrayAt(cropX0+x, cropY0+y))
		}
	}

	blockSize := max(31, int(mmToPx(10, dpi))|1)
	if blockSize%2 == 0 {
		blockSize++
	}
	binaryMask := adaptiveThreshold(sub, blockSize, 15)

	// Template-aware erasure: remove all known template elements by their computed positions
	EraseTemplateFromMask(binaryMask, cropW, cropH, xform, isTemplate, cropX0, cropY0)

	// Light morphological opening to clean up noise (NOT for removing template elements)
	noiseRadius := max(1, int(math.Round(mmToPx(0.15, dpi))))
	cleaned := morphOpen(binaryMask, cropW, cropH, noiseRadius)

	contour := traceLargestContour(cleaned, cropW, cropH)
	if len(contour) < 3 {
		return nil, fmt.Errorf("no contour found")
	}

	fContour := make([][2]float64, len(contour))
	for i, pt := range contour {
		fContour[i] = [2]float64{float64(pt[0] + cropX0), float64(pt[1] + cropY0)}
	}

	tolerance := mmToPx(0.5, dpi)
	simplified := douglasPeucker(fContour, tolerance)
	if len(simplified) < 3 {
		return nil, fmt.Errorf("contour too simple after DP: %d points", len(simplified))
	}

	cellCenterX := cell.X + cell.W/2
	cellCenterY := cell.Y + cell.H/2
	mmOutline := make([][2]float64, len(simplified))
	for i, pt := range simplified {
		mx, my := xform.transform(pt[0], pt[1])
		mmOutline[i] = [2]float64{mx - cellCenterX, my - cellCenterY}
	}

	return mmOutline, nil
}

// --- face outline extraction (legacy fallback) ---

func extractFaceOutline(img *image.Gray, xform affineTransform, cell faceCell, dpi float64) ([][2]float64, error) {
	inv := invertAffine(xform)

	// Cell bounds in mm (with padding)
	pad := 2.0 // mm extra around cell
	mmX0, mmY0 := cell.X-pad, cell.Y-pad
	mmX1, mmY1 := cell.X+cell.W+pad, cell.Y+cell.H+pad

	// Convert corners to pixel coords
	corners := [4][2]float64{
		{mmX0, mmY0}, {mmX1, mmY0}, {mmX0, mmY1}, {mmX1, mmY1},
	}
	var pxMinX, pxMinY float64 = math.MaxFloat64, math.MaxFloat64
	var pxMaxX, pxMaxY float64 = -math.MaxFloat64, -math.MaxFloat64
	for _, c := range corners {
		px, py := inv.transform(c[0], c[1])
		pxMinX = math.Min(pxMinX, px)
		pxMinY = math.Min(pxMinY, py)
		pxMaxX = math.Max(pxMaxX, px)
		pxMaxY = math.Max(pxMaxY, py)
	}

	bounds := img.Bounds()
	cropX0 := max(bounds.Min.X, int(pxMinX))
	cropY0 := max(bounds.Min.Y, int(pxMinY))
	cropX1 := min(bounds.Max.X, int(pxMaxX)+1)
	cropY1 := min(bounds.Max.Y, int(pxMaxY)+1)
	cropW := cropX1 - cropX0
	cropH := cropY1 - cropY0

	if cropW < 10 || cropH < 10 {
		return nil, fmt.Errorf("cell crop too small: %dx%d", cropW, cropH)
	}

	// Extract sub-image
	sub := image.NewGray(image.Rect(0, 0, cropW, cropH))
	for y := 0; y < cropH; y++ {
		for x := 0; x < cropW; x++ {
			sub.SetGray(x, y, img.GrayAt(cropX0+x, cropY0+y))
		}
	}

	// Adaptive threshold
	blockSize := max(31, int(mmToPx(10, dpi))|1) // ensure odd
	if blockSize%2 == 0 {
		blockSize++
	}
	binaryMask := adaptiveThreshold(sub, blockSize, 15)

	// Morphological opening to remove template elements:
	// Bracket corners are 1.0mm stroke, dashed border is 0.6mm stroke.
	// Erosion radius must exceed half the bracket width to fully remove it.
	// At 600 DPI: 0.6mm = 14px, 1.0mm = 24px → need > 12px erosion.
	// User's traced line (Sharpie/pen) is ~2-3mm = 47-71px → survives easily.
	erosionRadius := max(3, int(math.Round(mmToPx(0.6, dpi))))
	cleaned := morphOpen(binaryMask, cropW, cropH, erosionRadius)

	// Erase anything within 3mm of cell edges to remove residual bracket/border fragments.
	maskMarginMM := 3.0
	for py := 0; py < cropH; py++ {
		for px := 0; px < cropW; px++ {
			if !cleaned[py*cropW+px] {
				continue
			}
			mmX, mmY := xform.transform(float64(px+cropX0), float64(py+cropY0))
			if mmX-cell.X < maskMarginMM || cell.X+cell.W-mmX < maskMarginMM ||
				mmY-cell.Y < maskMarginMM || cell.Y+cell.H-mmY < maskMarginMM {
				cleaned[py*cropW+px] = false
			}
		}
	}

	// Find the largest contour
	contour := traceLargestContour(cleaned, cropW, cropH)
	if len(contour) < 3 {
		return nil, fmt.Errorf("no contour found")
	}

	// Convert pixel contour to float
	fContour := make([][2]float64, len(contour))
	for i, pt := range contour {
		fContour[i] = [2]float64{float64(pt[0] + cropX0), float64(pt[1] + cropY0)}
	}

	// Simplify
	tolerance := mmToPx(0.5, dpi)
	simplified := douglasPeucker(fContour, tolerance)
	if len(simplified) < 3 {
		return nil, fmt.Errorf("contour too simple after DP: %d points", len(simplified))
	}

	// Convert to mm relative to cell center
	cellCenterX := cell.X + cell.W/2
	cellCenterY := cell.Y + cell.H/2
	mmOutline := make([][2]float64, len(simplified))
	for i, pt := range simplified {
		mx, my := xform.transform(pt[0], pt[1])
		mmOutline[i] = [2]float64{mx - cellCenterX, my - cellCenterY}
	}

	return mmOutline, nil
}

// --- contour tracing (Moore neighborhood) ---

func traceLargestContour(mask []bool, w, h int) [][2]int {
	// Find all outer contours, return the largest by area
	visited := make([]bool, w*h)
	var largest [][2]int

	for y := 1; y < h-1; y++ {
		for x := 1; x < w-1; x++ {
			if !mask[y*w+x] || visited[y*w+x] {
				continue
			}
			// Check if it's a boundary pixel (has at least one non-set neighbor)
			isBoundary := false
			for _, d := range [][2]int{{-1, 0}, {1, 0}, {0, -1}, {0, 1}} {
				nx, ny := x+d[0], y+d[1]
				if !mask[ny*w+nx] {
					isBoundary = true
					break
				}
			}
			if !isBoundary {
				continue
			}

			contour := mooreTrace(mask, w, h, x, y)
			for _, pt := range contour {
				visited[pt[1]*w+pt[0]] = true
			}

			if len(contour) > len(largest) {
				largest = contour
			}
		}
	}

	return largest
}

func mooreTrace(mask []bool, w, h, startX, startY int) [][2]int {
	// 8-connected Moore neighborhood tracing
	// Direction offsets: 0=E, 1=SE, 2=S, 3=SW, 4=W, 5=NW, 6=N, 7=NE
	dx := [8]int{1, 1, 0, -1, -1, -1, 0, 1}
	dy := [8]int{0, 1, 1, 1, 0, -1, -1, -1}

	var contour [][2]int
	contour = append(contour, [2]int{startX, startY})

	// Start direction: came from the west (direction 0 is east)
	dir := 7 // start looking NE (entry from west means backtrack = west = dir 4, so start at (4+1)%8 = 5... actually standard: start at (backtrack+1) mod 8)
	// Standard: if we found the pixel by scanning left-to-right, entry was from west, so backtrack direction is 4 (W), start searching at (4+1)%8 = 5 (NW)
	dir = 5

	cx, cy := startX, startY
	maxIter := w * h * 2

	for i := 0; i < maxIter; i++ {
		found := false
		for j := 0; j < 8; j++ {
			nd := (dir + j) % 8
			nx, ny := cx+dx[nd], cy+dy[nd]
			if nx >= 0 && nx < w && ny >= 0 && ny < h && mask[ny*w+nx] {
				cx, cy = nx, ny
				contour = append(contour, [2]int{cx, cy})
				// Backtrack direction
				dir = (nd + 5) % 8 // opposite of nd is (nd+4)%8, start at (nd+4+1)%8 = (nd+5)%8
				found = true
				break
			}
		}

		if !found || (cx == startX && cy == startY) {
			break
		}
	}

	return contour
}

// --- Douglas-Peucker simplification ---

func douglasPeucker(pts [][2]float64, epsilon float64) [][2]float64 {
	if len(pts) <= 2 {
		return pts
	}

	maxDist := 0.0
	maxIdx := 0
	end := len(pts) - 1

	for i := 1; i < end; i++ {
		d := pointToSegmentDist(pts[i], pts[0], pts[end])
		if d > maxDist {
			maxDist = d
			maxIdx = i
		}
	}

	if maxDist > epsilon {
		left := douglasPeucker(pts[:maxIdx+1], epsilon)
		right := douglasPeucker(pts[maxIdx:], epsilon)
		return append(left[:len(left)-1], right...)
	}

	return [][2]float64{pts[0], pts[end]}
}

func pointToSegmentDist(p, a, b [2]float64) float64 {
	abx, aby := b[0]-a[0], b[1]-a[1]
	apx, apy := p[0]-a[0], p[1]-a[1]
	ab2 := abx*abx + aby*aby
	if ab2 < 1e-12 {
		return math.Hypot(apx, apy)
	}
	t := (apx*abx + apy*aby) / ab2
	t = math.Max(0, math.Min(1, t))
	projX := a[0] + t*abx
	projY := a[1] + t*aby
	return math.Hypot(p[0]-projX, p[1]-projY)
}

// --- utility ---

// sortFaces sorts traced faces by face number.
func sortFaces(faces []TracedFace) {
	sort.Slice(faces, func(i, j int) bool {
		return faces[i].FaceNum < faces[j].FaceNum
	})
}