7.2 直线光栅化的误差界

#!/usr/bin/env python3
# file: lantern/sim/cg_raster_07_2_error.py
# tier: S
# verify: /opt/homebrew/bin/python3 sim/cg_raster_07_2_error.py
"""7.2 直线光栅化的误差界：把「半像素」写成一条纯整数判据并穷举验证。

第 i 步的误差 e_i = k_i - i*minor/major（像素单位）。两边同乘 2*major 消掉除法：

    2*|k_i*major - i*minor| <= major        <=>        |e_i| <= 1/2

左边是纯整数，故界可以逐组穷举判定，不需要任何浮点。
穷举 major in [1,63]、minor in [0,major]、i in [0,major] 共 89439 个 (i, 斜率) 采样点。
产物：out/cg_raster_07_2_error.csv、out/cg_raster_07_2_error.svg、out/cg_raster_07_2_error.json。
"""

import json
import os
import sys
import xml.etree.ElementTree as ET

sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))
import cglib  # noqa: E402

ROOT = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))
MAXM = 63


def kseq(major, minor, bias=None):
    """minor 偏移序列 k_i：bias=None 走整数 Bresenham（d > 0 进格，平局不动），
    否则走定点 DDA，取整为 (acc + bias) >> 16。"""
    out = []
    if bias is None:
        d, k = 2 * minor - major, 0
        for _ in range(major + 1):
            out.append(k)
            if d > 0:
                k += 1
                d += 2 * (minor - major)
            else:
                d += 2 * minor
    else:
        step, acc = cglib.fdiv(minor << 16, major << 16), 0
        for _ in range(major + 1):
            out.append((acc + bias) >> 16)
            acc += step
    return out


def residual(k, major, minor):
    """(max|2*(k*major - i*minor)|, 误差的 Q16.16 最小值, 最大值)。"""
    best, emin, emax = 0, 0, 0
    for i, ki in enumerate(k):
        r = ki * major - i * minor
        best = max(best, 2 * abs(r))
        e = cglib.fdiv(r << 16, major << 16)
        emin, emax = min(emin, e), max(emax, e)
    return best, emin, emax


def build_svg(path):
    """误差曲线 + ±1/2 包络（纯 xml.etree，不用任何绘图库）。"""
    W, H, pad = 560, 240, 24
    svg = ET.Element("svg", {"xmlns": "http://www.w3.org/2000/svg", "width": str(W), "height": str(H)})
    ET.SubElement(svg, "rect", {"x": "0", "y": "0", "width": str(W), "height": str(H), "fill": "#101014"})
    y0, span = H // 2, (H - 2 * pad) // 2
    for sign in (-1, 1):
        ET.SubElement(svg, "line", {"x1": str(pad), "y1": str(y0 + sign * span // 2),
                                    "x2": str(W - pad), "y2": str(y0 + sign * span // 2),
                                    "stroke": "#ff9a3c", "stroke-dasharray": "4 4"})
    for (major, minor), col in zip([(8, 1), (6, 1), (42, 17)], ["#6ad4ff", "#9dff6a", "#ff6a9d"]):
        pts = []
        for i, ki in enumerate(kseq(major, minor)):
            e = (ki * major - i * minor) / major  # 仅用于绘图，判据本身是纯整数
            pts.append(f"{pad + (W - 2 * pad) * i // major},{y0 - int(e * span)}")
        ET.SubElement(svg, "polyline", {"points": " ".join(pts), "fill": "none",
                                        "stroke": col, "stroke-width": "2"})
    label = ET.SubElement(svg, "text", {"x": str(pad), "y": str(pad - 8),
                                        "fill": "#e8e8f0", "font-size": "12"})
    label.text = "e_i 曲线与 +-1/2 包络（斜率 8/1、6/1、42/17）"
    os.makedirs(os.path.dirname(path), exist_ok=True)
    ET.ElementTree(svg).write(path, encoding="utf-8", xml_declaration=True)
    ET.fromstring(open(path, encoding="utf-8").read())  # S 档断言：SVG 必须能被解析
    return os.path.getsize(path)


def main():
    rows, tight, samples = [], [], 0
    stat = {"bresenham": [0, 0], "dda_round": [0, 0], "dda_floor": [0, 0]}
    for major in range(1, MAXM + 1):
        for minor in range(0, major + 1):
            samples += major + 1
            vals = {}
            for name, k in (("bresenham", kseq(major, minor)),
                            ("dda_round", kseq(major, minor, 32768)),
                            ("dda_floor", kseq(major, minor, 0))):
                vals[name] = residual(k, major, minor)
                stat[name][0] += vals[name][0] > major   # 判据：2*|res| <= major 即 |e| <= 1/2
                stat[name][1] = max(stat[name][1], vals[name][0] / (2 * major))
            if vals["bresenham"][0] == major:            # 界被恰好达到 -> 界是紧的
                tight.append([major, minor, major // 2, vals["bresenham"][2]])
            rows.append((major, minor, vals["bresenham"][0], vals["dda_round"][0], vals["dda_floor"][0],
                         vals["dda_floor"][1], vals["dda_floor"][2]))

    # 逐像素比对：Bresenham 与 DDA(+32768) 的分歧只在平局处，且两边 |e| 都是 1/2
    div = [m for m in range(1, MAXM + 1) if kseq(m, 1) != kseq(m, 1, 32768)]
    kb, kd = kseq(8, 1), kseq(8, 1, 32768)
    res = lambda k: [2 * (k[i] * 8 - i) for i in range(9)]  # noqa: E731

    out = os.path.join(ROOT, "out")
    cglib.write_csv(os.path.join(out, "cg_raster_07_2_error.csv"),
                    ["major", "minor", "res_bresenham", "res_dda_round", "res_dda_floor",
                     "e_floor_min_q16", "e_floor_max_q16"], rows)
    svg = build_svg(os.path.join(out, "cg_raster_07_2_error.svg"))
    with open(os.path.join(out, "cg_raster_07_2_error.json"), "w", encoding="utf-8") as fh:
        json.dump({"max_major": MAXM, "samples": samples, "rows": len(rows),
                   "violations": {n: stat[n][0] for n in stat},
                   "max_abs_e_half": {n: stat[n][1] for n in stat},
                   "tight_cases": len(tight), "tight_examples": tight[:6],
                   "demo_major": 8, "demo_minor": 1, "demo_k_bresenham": kb, "demo_k_dda_round": kd,
                   "demo_residual_bresenham": res(kb), "demo_residual_dda_round": res(kd),
                   "diverge_minor1": div, "svg_bytes": svg}, fh, ensure_ascii=False, indent=2)

    print(f"采样点         : {samples} 个 (i, 斜率)，覆盖 {len(rows)} 组 (major, minor)")
    print("判据 2*|k*major - i*minor| <= major  (= |e| <= 1/2 像素)：")
    for n in ("bresenham", "dda_round", "dda_floor"):
        print(f"  {n:10s} 越界 {stat[n][0]:5d} 组   实测 max|e| = {stat[n][1]:.4f} 像素")
    print(f"界是紧的       : {len(tight)} 组恰好达到 |e| = 1/2，例 {tight[0][:3]}")
    print(f"分歧只在平局    : minor=1 时 Bresenham 与 DDA(+32768) 分歧的 major = {div}")
    print("  演示 major=8 minor=1 :")
    print(f"    k_bresenham = {kb}  残差 2*(k*major - i*minor) = {res(kb)}")
    print(f"    k_dda_round = {kd}  残差 = {res(kd)}")
    print("    平局处 i=4 两边误差都是 1/2，符号相反；两边都满足界，逐像素比对才分得开")
    print(f"产物           : error.csv {len(rows)} 行 / error.svg {svg} bytes")


if __name__ == "__main__":
    main()

# snippet: 复现命令（在 lantern 工作区根执行）
cd $HOME/.cgbook/lantern
/opt/homebrew/bin/python3 sim/cg_raster_07_2_error.py
node --test 'core/test/cg_raster_07_2_error.test.ts'

// file: lantern/core/test/cg_raster_07_2_error.test.ts
// tier: P
// verify: node --test 'core/test/cg_raster_07_2_error.test.ts'
// 7.2 的跨语言断言：误差界 2*|k*major - i*minor| <= major 的穷举验证必须与
// Python 侧 sim/cg_raster_07_2_error.py 逐数一致，且负例必须被同一判据拒绝。

import assert from "node:assert/strict";
import { readFileSync } from "node:fs";
import { dirname, join } from "node:path";
import { test } from "node:test";
import { fileURLToPath } from "node:url";

import { bresenhamLine, ddaLine } from "../src/lantern/cg_raster_07_line.ts";

type G = {
  max_major: number; samples: number; rows: number;
  violations: Record<string, number>; max_abs_e_half: Record<string, number>;
  tight_cases: number; demo_major: number; demo_minor: number;
  demo_k_bresenham: number[]; demo_k_dda_round: number[];
  demo_residual_bresenham: number[]; demo_residual_dda_round: number[];
  diverge_minor1: number[];
};
type Fn = (a: number, b: number, c: number, d: number) => [number, number][];

const root = join(dirname(fileURLToPath(import.meta.url)), "..", "..");
const g = JSON.parse(readFileSync(join(root, "out", "cg_raster_07_2_error.json"), "utf8")) as G;
/** 取 minor 偏移序列 k_i：斜率 0 <= minor <= major 时主轴是 x，故像素的 y 就是 k_i。 */
const ks = (fn: Fn, major: number, minor: number): number[] => fn(0, 0, major, minor).map((p) => p[1]);
/** 判据：2*|k*major - i*minor| <= num*major 即 |e| <= num/2。num = 1 是半像素界。 */
const holds = (k: number[], major: number, minor: number, num: number): boolean =>
  k.every((ki, i) => 2 * Math.abs(ki * major - i * minor) <= num * major);
const ddaR: Fn = (a, b, c, d) => ddaLine(a, b, c, d, 32768);
const ddaF: Fn = (a, b, c, d) => ddaLine(a, b, c, d, 0);

test("穷举：半像素界在 89439 个采样点上成立，越界组数与 Python 侧一致", () => {
  let samples = 0;
  let rows = 0;
  let vB = 0;
  let vR = 0;
  let vF = 0;
  let tight = 0;
  let maxE = 0;
  for (let major = 1; major <= g.max_major; major += 1) {
    for (let minor = 0; minor <= major; minor += 1) {
      rows += 1;
      samples += major + 1;
      const kb = ks(bresenhamLine, major, minor);
      const kr = ks(ddaR, major, minor);
      if (!holds(kb, major, minor, 1)) vB += 1;
      if (!holds(kr, major, minor, 1)) vR += 1;
      if (!holds(ks(ddaF, major, minor), major, minor, 1)) vF += 1;
      for (const k of [kb, kr]) {
        k.forEach((ki, i) => { maxE = Math.max(maxE, Math.abs(ki * major - i * minor) / major); });
      }
      if (kb.some((ki, i) => 2 * Math.abs(ki * major - i * minor) === major)) tight += 1;
    }
  }
  assert.equal(samples, g.samples);
  assert.equal(rows, g.rows);
  assert.equal(vB, g.violations.bresenham);
  assert.equal(vR, g.violations.dda_round);
  assert.equal(vF, g.violations.dda_floor);
  assert.equal(maxE, g.max_abs_e_half.bresenham);
  assert.equal(tight, g.tight_cases);
  // 界是紧的：恰有 651 组（= 7.1 的平局组数）达到 |e| = 1/2
  assert.ok(tight > 0 && maxE === 0.5, "界必须被达到，否则不是紧的");
});

test("平局处的对称分叉：i=4 两边误差都是 1/2，符号相反", () => {
  const { demo_major: M, demo_minor: m } = g;
  const kb = ks(bresenhamLine, M, m);
  const kr = ks(ddaR, M, m);
  assert.deepEqual(kb, g.demo_k_bresenham);
  assert.deepEqual(kr, g.demo_k_dda_round);
  assert.deepEqual(kb.map((k, i) => 2 * (k * M - i * m)), g.demo_residual_bresenham);
  assert.deepEqual(kr.map((k, i) => 2 * (k * M - i * m)), g.demo_residual_dda_round);
  // 两个残差数组只在平局点分叉，且互为相反数
  assert.deepEqual(kb.map((_, i) => i).filter((i) => kb[i] !== kr[i]), [4]);
  assert.equal(g.demo_residual_bresenham[4], -M);
  assert.equal(g.demo_residual_dda_round[4], M);
});

test("负例 A：把界放宽 1 倍（|e| <= 1）-> 判据不再拒绝 floor DDA", () => {
  let rejectedTight = 0;
  let rejectedLoose = 0;
  for (let major = 1; major <= g.max_major; major += 1) {
    for (let minor = 0; minor <= major; minor += 1) {
      const kf = ks(ddaF, major, minor);
      if (!holds(kf, major, minor, 1)) rejectedTight += 1;   // 1/2 界必须拒绝 floor DDA
      if (!holds(kf, major, minor, 2)) rejectedLoose += 1;   // 放宽到 1 像素后应全部通过
    }
  }
  assert.equal(rejectedTight, g.violations.dda_floor);
  assert.ok(rejectedTight > 0, "1/2 界必须拒绝 floor DDA");
  assert.equal(rejectedLoose, 0, "放宽到 1 像素后判据失效：floor DDA 全部通过");
});

test("负例 B：误差界不能替代逐像素比对（分歧只在平局，两边都满足界）", () => {
  // minor = 1 时两者分歧的 major 全是 2 的幂：平局恰好落在整步长上
  const div: number[] = [];
  for (let major = 1; major <= g.max_major; major += 1) {
    if (JSON.stringify(ks(bresenhamLine, major, 1)) !== JSON.stringify(ks(ddaR, major, 1))) div.push(major);
  }
  assert.deepEqual(div, g.diverge_minor1);
  assert.ok(div.length > 0);
  // 两边都满足 1/2 界，却给出不同像素 —— 故 sha256 比对不可省略
  for (const major of div) {
    assert.ok(holds(ks(bresenhamLine, major, 1), major, 1, 1));
    assert.ok(holds(ks(ddaR, major, 1), major, 1, 1));
  }
});