package main

import (
	"fmt"
	"image"
	"math"
	"math/rand"
	"testing"
)

// --- ground truth shapes (mm relative to center) ---

var squareMM = func(side float64) [][2]float64 {
	h := side / 2
	return [][2]float64{{-h, -h}, {h, -h}, {h, h}, {-h, h}}
}

var triangleMM = func(side float64) [][2]float64 {
	r := side / math.Sqrt(3)
	return [][2]float64{
		{0, -r},
		{r * math.Sqrt(3) / 2, r / 2},
		{-r * math.Sqrt(3) / 2, r / 2},
	}
}

var lShapeMM = func(size float64) [][2]float64 {
	s := size / 2
	return [][2]float64{
		{-s, -s}, {0, -s}, {0, 0}, {s, 0}, {s, s}, {-s, s},
	}
}

var starMM = func(outer float64) [][2]float64 {
	inner := outer * 0.4
	pts := make([][2]float64, 10)
	for i := 0; i < 5; i++ {
		a1 := float64(i)*2*math.Pi/5 - math.Pi/2
		a2 := a1 + math.Pi/5
		pts[2*i] = [2]float64{outer * math.Cos(a1), outer * math.Sin(a1)}
		pts[2*i+1] = [2]float64{inner * math.Cos(a2), inner * math.Sin(a2)}
	}
	return pts
}

var circleMM = func(r float64, n int) [][2]float64 {
	pts := make([][2]float64, n)
	for i := 0; i < n; i++ {
		a := 2 * math.Pi * float64(i) / float64(n)
		pts[i] = [2]float64{r * math.Cos(a), r * math.Sin(a)}
	}
	return pts
}

var irregularMM = func(size float64) [][2]float64 {
	s := size / 2
	return [][2]float64{
		{-s, -s * 0.3}, {-s * 0.4, -s}, {s * 0.6, -s * 0.8},
		{s, -s * 0.1}, {s * 0.7, s * 0.9}, {-s * 0.2, s},
		{-s * 0.8, s * 0.5},
	}
}

// --- accuracy metrics ---

func hausdorffDistance(a, b [][2]float64) float64 {
	directedMax := func(from, to [][2]float64) float64 {
		maxDist := 0.0
		for _, p := range from {
			minD := math.MaxFloat64
			for _, q := range to {
				d := math.Hypot(p[0]-q[0], p[1]-q[1])
				if d < minD {
					minD = d
				}
			}
			if minD > maxDist {
				maxDist = minD
			}
		}
		return maxDist
	}
	ab := directedMax(a, b)
	ba := directedMax(b, a)
	if ab > ba {
		return ab
	}
	return ba
}

func shoelaceArea(pts [][2]float64) float64 {
	n := len(pts)
	if n < 3 {
		return 0
	}
	sum := 0.0
	for i := 0; i < n; i++ {
		j := (i + 1) % n
		sum += pts[i][0]*pts[j][1] - pts[j][0]*pts[i][1]
	}
	return math.Abs(sum) / 2
}

func areaOverlapRatio(extracted, truth [][2]float64) float64 {
	ae := shoelaceArea(extracted)
	at := shoelaceArea(truth)
	if at == 0 {
		return 0
	}
	ratio := ae / at
	if ratio > 1 {
		return 1 / ratio
	}
	return ratio
}

// --- synthetic scan extensions for new marker system ---

func (s *syntheticScan) drawBullseyeMarkers(pageNum int) {
	positions := markerPositionsMM(s.cfg.PageWidth, s.cfg.PageHeight)
	cornerIDs := [4]int{CornerTL, CornerTR, CornerBL, CornerBR}

	for i, pos := range positions {
		grid := encodeMarkerGrid(MarkerData{
			PageNum:   pageNum,
			CornerID:  cornerIDs[i],
			NumFaces:  s.cfg.NumFaces,
			LongestMM: int(s.cfg.LongestSide),
		})
		cxPx := mmToPx(pos[0], s.dpi)
		cyPx := mmToPx(pos[1], s.dpi)
		cellPx := mmToPx(markerCellMM, s.dpi)
		halfN := float64(markerN) / 2.0

		// Compute cell boundaries cumulatively to avoid gaps/overlaps
		originX := cxPx - halfN*cellPx
		originY := cyPx - halfN*cellPx
		var colEdge [markerN + 1]int
		var rowEdge [markerN + 1]int
		for k := 0; k <= markerN; k++ {
			colEdge[k] = int(originX + float64(k)*cellPx)
			rowEdge[k] = int(originY + float64(k)*cellPx)
		}

		for r := 0; r < markerN; r++ {
			for c := 0; c < markerN; c++ {
				val := uint8(240)
				if grid[r][c] {
					val = 10
				}
				for py := rowEdge[r]; py < rowEdge[r+1]; py++ {
					for px := colEdge[c]; px < colEdge[c+1]; px++ {
						s.setPixelSafe(px, py, val)
					}
				}
			}
		}
	}
}

func (s *syntheticScan) drawDataBars() {
	bars := calibBarSpecs()
	mPos := markerPositionsMM(s.cfg.PageWidth, s.cfg.PageHeight)
	barY := s.cfg.PageHeight - 15.0
	bx := mPos[0][0] + float64(markerN)*markerCellMM/2 + 3

	for _, bar := range bars {
		cells := encodeCalibBar(bar.WidthMM)
		cellW := float64(bar.WidthMM) / float64(calibBarCells)
		y0 := int(mmToPx(barY, s.dpi))
		h := int(mmToPx(calibBarHeight, s.dpi))

		for c := 0; c < calibBarCells; c++ {
			cx0 := int(mmToPx(bx+float64(c)*cellW, s.dpi))
			cx1 := int(mmToPx(bx+float64(c+1)*cellW, s.dpi))
			val := uint8(240)
			if cells[c] {
				val = 10
			}
			for py := y0; py < y0+h; py++ {
				for px := cx0; px < cx1; px++ {
					s.setPixelSafe(px, py, val)
				}
			}
		}

		bx += float64(bar.WidthMM) + calibBarGap
	}
}

func (s *syntheticScan) fillRect(x, y, w, h int, value uint8) {
	for dy := 0; dy < h; dy++ {
		for dx := 0; dx < w; dx++ {
			s.setPixelSafe(x+dx, y+dy, value)
		}
	}
}

func (s *syntheticScan) drawTracedPolygonColor(cell faceCell, outlineMM [][2]float64, thickness int, intensity uint8) {
	cx := mmToPx(cell.X+cell.W/2, s.dpi)
	cy := mmToPx(cell.Y+cell.H/2, s.dpi)

	for i := 0; i < len(outlineMM); i++ {
		p0 := outlineMM[i]
		p1 := outlineMM[(i+1)%len(outlineMM)]
		x0 := cx + mmToPx(p0[0], s.dpi)
		y0 := cy + mmToPx(p0[1], s.dpi)
		x1 := cx + mmToPx(p1[0], s.dpi)
		y1 := cy + mmToPx(p1[1], s.dpi)
		s.drawThickLine(int(x0), int(y0), int(x1), int(y1), thickness, intensity)
	}
}

func (s *syntheticScan) drawCellBorderDashed(cell faceCell) {
	x0 := int(mmToPx(cell.X, s.dpi))
	y0 := int(mmToPx(cell.Y, s.dpi))
	x1 := int(mmToPx(cell.X+cell.W, s.dpi))
	y1 := int(mmToPx(cell.Y+cell.H, s.dpi))
	dashPx := int(mmToPx(1.2, s.dpi))
	if dashPx < 2 {
		dashPx = 2
	}

	for px := x0; px <= x1; px++ {
		if ((px-x0)/dashPx)%2 == 0 {
			s.setPixelSafe(px, y0, 170)
			s.setPixelSafe(px, y1, 170)
		}
	}
	for py := y0; py <= y1; py++ {
		if ((py-y0)/dashPx)%2 == 0 {
			s.setPixelSafe(x0, py, 170)
			s.setPixelSafe(x1, py, 170)
		}
	}
}

func (s *syntheticScan) drawCenterCross(cell faceCell) {
	cx := int(mmToPx(cell.X+cell.W/2, s.dpi))
	cy := int(mmToPx(cell.Y+cell.H/2, s.dpi))
	arm := int(mmToPx(3.0, s.dpi))
	for dx := -arm; dx <= arm; dx++ {
		s.setPixelSafe(cx+dx, cy, 210)
	}
	for dy := -arm; dy <= arm; dy++ {
		s.setPixelSafe(cx, cy+dy, 210)
	}
}

func (s *syntheticScan) drawFaceLabel(cell faceCell) {
	x := int(mmToPx(cell.X+2, s.dpi))
	y := int(mmToPx(cell.Y+2, s.dpi))
	w := int(mmToPx(8, s.dpi))
	h := int(mmToPx(6, s.dpi))
	s.fillRect(x, y, w, h, 200)
}

// buildNewScan creates a synthetic scan using the new marker system.
func buildNewScan(t *testing.T, dpi float64, cfg FaceTemplateConfig, pageNum int) *syntheticScan {
	t.Helper()
	s := newSyntheticScan(dpi, cfg)
	s.drawBullseyeMarkers(pageNum)
	s.drawDataBars()
	for _, cell := range s.layout.Cells {
		s.drawGridLinesInCell(cell)
		s.drawCellBorderDashed(cell)
		s.drawCenterCross(cell)
		s.drawFaceLabel(cell)
	}
	return s
}

// --- Group A: Marker detection ---

func TestHarnessMarkerDetection(t *testing.T) {
	dpis := []float64{150, 200, 300, 400, 600}
	cfg := FaceTemplateConfig{NumFaces: 4, LongestSide: 50, PageWidth: 215.9, PageHeight: 279.4}

	for _, dpi := range dpis {
		t.Run(fmt.Sprintf("DPI_%.0f", dpi), func(t *testing.T) {
			s := buildNewScan(t, dpi, cfg, 1)
			markers := DetectMarkers(s.img, dpi)
			if len(markers) < 3 {
				t.Errorf("expected >=3 markers at %.0f DPI, got %d", dpi, len(markers))
				return
			}

			corners := map[int]bool{}
			for _, m := range markers {
				corners[m.Data.CornerID] = true
				if m.Data.PageNum != 1 {
					t.Errorf("marker corner %d: page=%d, want 1", m.Data.CornerID, m.Data.PageNum)
				}
				if m.Data.NumFaces != 4 {
					t.Errorf("marker corner %d: faces=%d, want 4", m.Data.CornerID, m.Data.NumFaces)
				}
				if m.Data.LongestMM != 50 {
					t.Errorf("marker corner %d: longest=%d, want 50", m.Data.CornerID, m.Data.LongestMM)
				}
			}

			for _, cid := range []int{CornerTL, CornerTR, CornerBL, CornerBR} {
				if !corners[cid] {
					t.Errorf("corner %d not detected", cid)
				}
			}
		})
	}
}

func TestHarnessMarkerWithNoise(t *testing.T) {
	cfg := FaceTemplateConfig{NumFaces: 4, LongestSide: 50, PageWidth: 215.9, PageHeight: 279.4}
	noiseLevels := []float64{0.005, 0.01, 0.02, 0.03}

	for _, noise := range noiseLevels {
		t.Run(fmt.Sprintf("noise_%.1f_pct", noise*100), func(t *testing.T) {
			s := buildNewScan(t, 300, cfg, 1)
			rng := rand.New(rand.NewSource(42))
			s.addSaltPepperNoise(noise, rng)

			markers := DetectMarkers(s.img, 300)
			if noise <= 0.02 && len(markers) < 3 {
				t.Errorf("%.1f%% noise: expected >=3 markers, got %d", noise*100, len(markers))
			}
			t.Logf("%.1f%% noise: %d markers detected", noise*100, len(markers))
		})
	}
}

func TestHarnessMarkerRotated(t *testing.T) {
	cfg := FaceTemplateConfig{NumFaces: 4, LongestSide: 50, PageWidth: 215.9, PageHeight: 279.4}
	angles := []float64{0.5, 1.0, 1.5, 2.0, 3.0}

	for _, angle := range angles {
		t.Run(fmt.Sprintf("%.1f_deg", angle), func(t *testing.T) {
			s := buildNewScan(t, 300, cfg, 1)
			rotated := s.rotateSmall(angle)

			markers := DetectMarkers(rotated, 300)
			if angle <= 2.0 && len(markers) < 3 {
				t.Errorf("%.1f°: expected >=3 markers, got %d", angle, len(markers))
			}
			t.Logf("%.1f°: %d markers detected", angle, len(markers))
		})
	}
}

// --- Group B: Data bar DPI calibration ---

func TestHarnessDataBarCalibration(t *testing.T) {
	dpis := []float64{150, 300, 600}
	cfg := FaceTemplateConfig{NumFaces: 4, LongestSide: 50, PageWidth: 215.9, PageHeight: 279.4}

	for _, dpi := range dpis {
		t.Run(fmt.Sprintf("DPI_%.0f", dpi), func(t *testing.T) {
			s := buildNewScan(t, dpi, cfg, 1)

			// identity transform (synthetic scan has no distortion)
			srcPts := [4][2]float64{
				{mmToPx(15, dpi), mmToPx(15, dpi)},
				{mmToPx(cfg.PageWidth-15, dpi), mmToPx(15, dpi)},
				{mmToPx(15, dpi), mmToPx(cfg.PageHeight-15, dpi)},
				{mmToPx(cfg.PageWidth-15, dpi), mmToPx(cfg.PageHeight-15, dpi)},
			}
			dstPts := [4][2]float64{
				{15, 15},
				{cfg.PageWidth - 15, 15},
				{15, cfg.PageHeight - 15},
				{cfg.PageWidth - 15, cfg.PageHeight - 15},
			}
			xform := computeAffine(srcPts, dstPts)
			threshold := otsuThreshold(s.img)

			mPos := markerPositionsMM(cfg.PageWidth, cfg.PageHeight)
			barX := mPos[0][0] + float64(markerN)*markerCellMM/2 + 3
			measuredDPI := MeasureCalibBarDPI(s.img, threshold, xform, barX, cfg.PageHeight-15.0, calibBarSpecs())
			if measuredDPI < 50 {
				t.Errorf("measured DPI too low: %.1f (expected ~%.0f)", measuredDPI, dpi)
				return
			}

			errorPct := math.Abs(measuredDPI-dpi) / dpi * 100
			if errorPct > 15 {
				t.Errorf("DPI calibration error: %.1f%% (measured %.0f, expected %.0f)", errorPct, measuredDPI, dpi)
			}
			t.Logf("%.0f DPI: measured %.1f (%.1f%% error)", dpi, measuredDPI, errorPct)
		})
	}
}

// --- Group C: Template erasure ---

func TestHarnessTemplateErasure(t *testing.T) {
	cfg := FaceTemplateConfig{NumFaces: 4, LongestSide: 50, PageWidth: 215.9, PageHeight: 279.4, GridSpacing: 10}
	dpi := 300.0
	s := buildNewScan(t, dpi, cfg, 1)

	// Draw a thick diamond in cell 0
	cell := s.layout.Cells[0]
	shape := diamondMM
	thickPx := int(mmToPx(2.0, dpi))
	s.drawTracedPolygon(cell, shape, thickPx)

	// Build affine transform (identity for synthetic)
	xform := buildIdentityAffine(cfg, dpi)

	// Threshold and erase template
	blockSize := max(31, int(mmToPx(10, dpi))|1)
	if blockSize%2 == 0 {
		blockSize++
	}

	bounds := s.img.Bounds()
	w, h := bounds.Dx(), bounds.Dy()
	binaryMask := adaptiveThreshold(s.img, blockSize, 15)

	isTemplate := TemplateElementMap(cfg, s.layout)

	// Count black pixels before erasure
	blackBefore := 0
	for _, v := range binaryMask {
		if v {
			blackBefore++
		}
	}

	EraseTemplateFromMask(binaryMask, w, h, xform, isTemplate, 0, 0)

	// Count after
	blackAfter := 0
	for _, v := range binaryMask {
		if v {
			blackAfter++
		}
	}

	erased := blackBefore - blackAfter
	t.Logf("erasure: %d black before, %d after, %d erased (%.1f%%)",
		blackBefore, blackAfter, erased, float64(erased)/float64(blackBefore)*100)

	if blackAfter == 0 {
		t.Error("template erasure removed ALL black pixels including the traced shape")
	}
	if erased == 0 {
		t.Error("template erasure removed nothing")
	}

	// The traced shape pixels should survive — verify some remain in the cell area
	cellPxX0 := int(mmToPx(cell.X+5, dpi))
	cellPxY0 := int(mmToPx(cell.Y+5, dpi))
	cellPxX1 := int(mmToPx(cell.X+cell.W-5, dpi))
	cellPxY1 := int(mmToPx(cell.Y+cell.H-5, dpi))
	cellBlack := 0
	for py := cellPxY0; py < cellPxY1; py++ {
		for px := cellPxX0; px < cellPxX1; px++ {
			if px < w && py < h && binaryMask[py*w+px] {
				cellBlack++
			}
		}
	}
	if cellBlack < 100 {
		t.Errorf("only %d black pixels remaining in cell interior — traced shape may have been erased", cellBlack)
	}
	t.Logf("cell interior: %d black pixels remaining (traced shape)", cellBlack)
}

// --- Group D: End-to-end shape extraction ---

func TestHarnessE2EShapes(t *testing.T) {
	type shapeSpec struct {
		name   string
		points [][2]float64
		area   float64 // expected area in mm^2
	}

	shapes := []shapeSpec{
		{"diamond_20mm", diamondMM, 200},
		{"square_15mm", squareMM(15), 225},
		{"triangle_20mm", triangleMM(20), triangleArea(20)},
		{"L_shape_20mm", lShapeMM(20), 300},    // 3/4 of 20x20
		{"star_10mm", starMM(10), 0},           // complex, skip area check
		{"circle_12mm", circleMM(12, 24), 0},   // approximation
		{"irregular_18mm", irregularMM(18), 0}, // skip area check
	}

	for _, shape := range shapes {
		t.Run(shape.name, func(t *testing.T) {
			cfg := FaceTemplateConfig{NumFaces: 1, LongestSide: 50, PageWidth: 215.9, PageHeight: 279.4}
			dpi := 300.0
			s := buildNewScan(t, dpi, cfg, 1)
			cell := s.layout.Cells[0]
			thickPx := int(mmToPx(1.5, dpi))
			s.drawTracedPolygon(cell, shape.points, thickPx)

			dir := t.TempDir()
			path := saveTestImage(t, s.img, dir, "scan.png")

			result, err := ProcessFaceScans([]string{path}, cfg)
			if err != nil {
				t.Fatal(err)
			}
			if len(result.Faces) == 0 {
				t.Fatal("no faces extracted")
			}

			face := result.Faces[0]
			if len(face.Outline) < 3 {
				t.Fatalf("extracted outline has too few points: %d", len(face.Outline))
			}

			hd := hausdorffDistance(face.Outline, shape.points)
			areaRatio := areaOverlapRatio(face.Outline, shape.points)
			vertRatio := float64(len(face.Outline)) / float64(len(shape.points))

			t.Logf("shape=%s: hausdorff=%.2fmm, area_ratio=%.2f, verts=%d/%d (ratio=%.2f)",
				shape.name, hd, areaRatio, len(face.Outline), len(shape.points), vertRatio)

			if hd > 5.0 {
				t.Errorf("hausdorff distance too large: %.2fmm", hd)
			}
			if shape.area > 0 && areaRatio < 0.5 {
				t.Errorf("area ratio too low: %.2f", areaRatio)
			}
		})
	}
}

func triangleArea(side float64) float64 {
	return side * side * math.Sqrt(3) / 4
}

// --- Group E: Pen width variations ---

func TestHarnessPenWidths(t *testing.T) {
	widthsMM := []float64{0.5, 1.0, 1.5, 2.0, 3.0, 4.0, 5.0}
	cfg := FaceTemplateConfig{NumFaces: 1, LongestSide: 50, PageWidth: 215.9, PageHeight: 279.4}

	for _, wmm := range widthsMM {
		t.Run(fmt.Sprintf("%.1fmm", wmm), func(t *testing.T) {
			dpi := 300.0
			s := buildNewScan(t, dpi, cfg, 1)
			cell := s.layout.Cells[0]
			thickPx := max(1, int(mmToPx(wmm/2, dpi)))
			s.drawTracedPolygon(cell, squareMM(15), thickPx)

			dir := t.TempDir()
			path := saveTestImage(t, s.img, dir, "scan.png")

			result, err := ProcessFaceScans([]string{path}, cfg)
			if err != nil {
				t.Fatal(err)
			}

			if len(result.Faces) == 0 {
				t.Errorf("pen width %.1fmm: no faces extracted", wmm)
				return
			}
			face := result.Faces[0]
			hd := hausdorffDistance(face.Outline, squareMM(15))
			area := areaOverlapRatio(face.Outline, squareMM(15))
			t.Logf("pen=%.1fmm: hausdorff=%.2fmm, area=%.2f, verts=%d", wmm, hd, area, len(face.Outline))

			if wmm >= 1.0 && wmm <= 4.0 && hd > 6.0 {
				t.Errorf("pen %.1fmm should extract well, hausdorff=%.2fmm", wmm, hd)
			}
		})
	}
}

// --- Group F: Pen intensity variations ---

func TestHarnessPenIntensity(t *testing.T) {
	intensities := []uint8{10, 30, 50, 80, 100, 120, 140}
	cfg := FaceTemplateConfig{NumFaces: 1, LongestSide: 50, PageWidth: 215.9, PageHeight: 279.4}

	for _, intensity := range intensities {
		t.Run(fmt.Sprintf("gray_%d", intensity), func(t *testing.T) {
			dpi := 300.0
			s := buildNewScan(t, dpi, cfg, 1)
			cell := s.layout.Cells[0]
			thickPx := int(mmToPx(1.0, dpi))
			s.drawTracedPolygonColor(cell, squareMM(15), thickPx, intensity)

			dir := t.TempDir()
			path := saveTestImage(t, s.img, dir, "scan.png")

			result, err := ProcessFaceScans([]string{path}, cfg)
			if err != nil {
				t.Fatal(err)
			}

			extracted := len(result.Faces) > 0 && len(result.Faces[0].Outline) >= 3
			t.Logf("intensity=%d: extracted=%v", intensity, extracted)

			if intensity <= 100 && !extracted {
				t.Errorf("intensity %d should be dark enough to extract", intensity)
			}
		})
	}
}

// --- Group G: Noise levels with shape extraction ---

func TestHarnessNoiseE2E(t *testing.T) {
	noiseLevels := []float64{0.001, 0.005, 0.01, 0.02, 0.03, 0.05}
	cfg := FaceTemplateConfig{NumFaces: 1, LongestSide: 50, PageWidth: 215.9, PageHeight: 279.4}

	for _, noise := range noiseLevels {
		t.Run(fmt.Sprintf("%.1f_pct", noise*100), func(t *testing.T) {
			dpi := 300.0
			s := buildNewScan(t, dpi, cfg, 1)
			cell := s.layout.Cells[0]
			thickPx := int(mmToPx(1.5, dpi))
			s.drawTracedPolygon(cell, squareMM(15), thickPx)

			rng := rand.New(rand.NewSource(42))
			s.addSaltPepperNoise(noise, rng)

			dir := t.TempDir()
			path := saveTestImage(t, s.img, dir, "scan.png")

			result, err := ProcessFaceScans([]string{path}, cfg)
			if err != nil {
				t.Fatal(err)
			}

			if len(result.Faces) > 0 {
				face := result.Faces[0]
				hd := hausdorffDistance(face.Outline, squareMM(15))
				t.Logf("noise=%.1f%%: hausdorff=%.2fmm, verts=%d", noise*100, hd, len(face.Outline))
			} else {
				t.Logf("noise=%.1f%%: no faces extracted", noise*100)
				if noise <= 0.02 {
					t.Errorf("%.1f%% noise should not prevent extraction", noise*100)
				}
			}
		})
	}
}

// --- Group H: Geometric distortions ---

func TestHarnessRotationE2E(t *testing.T) {
	angles := []float64{0, 0.5, 1.0, 1.5, 2.0, 3.0}
	cfg := FaceTemplateConfig{NumFaces: 1, LongestSide: 50, PageWidth: 215.9, PageHeight: 279.4}

	for _, angle := range angles {
		t.Run(fmt.Sprintf("%.1f_deg", angle), func(t *testing.T) {
			dpi := 300.0
			s := buildNewScan(t, dpi, cfg, 1)
			cell := s.layout.Cells[0]
			thickPx := int(mmToPx(1.5, dpi))
			s.drawTracedPolygon(cell, squareMM(15), thickPx)

			var scanImg *image.Gray
			if angle == 0 {
				scanImg = s.img
			} else {
				scanImg = s.rotateSmall(angle)
			}

			dir := t.TempDir()
			path := saveTestImage(t, scanImg, dir, "scan.png")

			result, err := ProcessFaceScans([]string{path}, cfg)
			if err != nil {
				t.Fatal(err)
			}

			if len(result.Faces) > 0 {
				face := result.Faces[0]
				hd := hausdorffDistance(face.Outline, squareMM(15))
				t.Logf("angle=%.1f°: hausdorff=%.2fmm, verts=%d", angle, hd, len(face.Outline))
			} else {
				t.Logf("angle=%.1f°: no faces extracted", angle)
				if angle <= 2.0 {
					t.Errorf("%.1f° rotation should not prevent extraction", angle)
				}
			}
		})
	}
}

// --- Group I: Multi-cell extraction ---

func TestHarnessMultiCell(t *testing.T) {
	cfg := FaceTemplateConfig{NumFaces: 4, LongestSide: 50, PageWidth: 215.9, PageHeight: 279.4}
	dpi := 300.0
	s := buildNewScan(t, dpi, cfg, 1)

	shapes := [][][2]float64{
		squareMM(15),
		diamondMM,
		triangleMM(18),
		lShapeMM(16),
	}
	for i, cell := range s.layout.Cells {
		if i < len(shapes) {
			thickPx := int(mmToPx(1.5, dpi))
			s.drawTracedPolygon(cell, shapes[i], thickPx)
		}
	}

	dir := t.TempDir()
	path := saveTestImage(t, s.img, dir, "scan.png")

	result, err := ProcessFaceScans([]string{path}, cfg)
	if err != nil {
		t.Fatal(err)
	}

	t.Logf("multi-cell: %d faces extracted from %d cells", len(result.Faces), len(s.layout.Cells))
	for _, face := range result.Faces {
		t.Logf("  face %d: %d vertices", face.FaceNum, len(face.Outline))
	}

	if len(result.Faces) < 2 {
		t.Errorf("expected at least 2 faces from 4 cells, got %d", len(result.Faces))
	}
}

// --- Group J: Edge cases ---

func TestHarnessEmptyCell(t *testing.T) {
	cfg := FaceTemplateConfig{NumFaces: 2, LongestSide: 50, PageWidth: 215.9, PageHeight: 279.4}
	dpi := 300.0
	s := buildNewScan(t, dpi, cfg, 1)

	// Only trace in cell 0, leave cell 1 empty
	thickPx := int(mmToPx(1.5, dpi))
	s.drawTracedPolygon(s.layout.Cells[0], squareMM(15), thickPx)

	dir := t.TempDir()
	path := saveTestImage(t, s.img, dir, "scan.png")

	result, err := ProcessFaceScans([]string{path}, cfg)
	if err != nil {
		t.Fatal(err)
	}

	// Should extract face 1 but not face 2
	hasFace1 := false
	for _, f := range result.Faces {
		if f.FaceNum == 1 {
			hasFace1 = true
		}
	}
	if !hasFace1 {
		t.Error("face 1 (traced) should be extracted")
	}
	t.Logf("empty cell: %d faces extracted (expected 1)", len(result.Faces))
}

func TestHarnessTinyShape(t *testing.T) {
	cfg := FaceTemplateConfig{NumFaces: 1, LongestSide: 50, PageWidth: 215.9, PageHeight: 279.4}
	dpi := 300.0
	s := buildNewScan(t, dpi, cfg, 1)

	// Very small 3mm square
	thickPx := int(mmToPx(0.5, dpi))
	s.drawTracedPolygon(s.layout.Cells[0], squareMM(3), thickPx)

	dir := t.TempDir()
	path := saveTestImage(t, s.img, dir, "scan.png")

	result, err := ProcessFaceScans([]string{path}, cfg)
	if err != nil {
		t.Fatal(err)
	}
	t.Logf("tiny shape: %d faces extracted", len(result.Faces))
}

func TestHarnessLargeShape(t *testing.T) {
	cfg := FaceTemplateConfig{NumFaces: 1, LongestSide: 80, PageWidth: 215.9, PageHeight: 279.4}
	dpi := 300.0
	s := buildNewScan(t, dpi, cfg, 1)

	// Shape that nearly fills the cell
	thickPx := int(mmToPx(2.0, dpi))
	s.drawTracedPolygon(s.layout.Cells[0], squareMM(70), thickPx)

	dir := t.TempDir()
	path := saveTestImage(t, s.img, dir, "scan.png")

	result, err := ProcessFaceScans([]string{path}, cfg)
	if err != nil {
		t.Fatal(err)
	}
	if len(result.Faces) > 0 {
		face := result.Faces[0]
		hd := hausdorffDistance(face.Outline, squareMM(70))
		t.Logf("large shape: hausdorff=%.2fmm, verts=%d", hd, len(face.Outline))
	} else {
		t.Error("large shape should be extractable")
	}
}

func TestHarnessDPIConsistency(t *testing.T) {
	cfg := FaceTemplateConfig{NumFaces: 1, LongestSide: 50, PageWidth: 215.9, PageHeight: 279.4}
	shape := squareMM(15)
	dpis := []float64{150, 300, 600}
	var areas []float64

	for _, dpi := range dpis {
		s := buildNewScan(t, dpi, cfg, 1)
		thickPx := int(mmToPx(1.5, dpi))
		s.drawTracedPolygon(s.layout.Cells[0], shape, thickPx)

		dir := t.TempDir()
		path := saveTestImage(t, s.img, dir, "scan.png")

		result, err := ProcessFaceScans([]string{path}, cfg)
		if err != nil {
			t.Errorf("%.0f DPI: %v", dpi, err)
			continue
		}
		if len(result.Faces) == 0 {
			t.Errorf("%.0f DPI: no faces extracted", dpi)
			continue
		}
		a := shoelaceArea(result.Faces[0].Outline)
		areas = append(areas, a)
		t.Logf("%.0f DPI: area=%.1f mm^2, verts=%d", dpi, a, len(result.Faces[0].Outline))
	}

	if len(areas) >= 2 {
		for i := 1; i < len(areas); i++ {
			ratio := areas[i] / areas[0]
			if ratio < 0.5 || ratio > 2.0 {
				t.Errorf("area at %.0f DPI (%.1f) vs %.0f DPI (%.1f): ratio=%.2f, too inconsistent",
					dpis[i], areas[i], dpis[0], areas[0], ratio)
			}
		}
	}
}

// --- helpers ---

func buildIdentityAffine(cfg FaceTemplateConfig, dpi float64) affineTransform {
	// Maps pixel coords to mm coords via simple scale (no rotation/skew)
	scale := 25.4 / dpi
	return affineTransform{
		a: scale, b: 0, c: 0,
		d: 0, e: scale, f: 0,
	}
}