1// Package blend finds how two sources' pictures of the same ground line up.2//3// Two captures of one game map draw the same places at different scales and4// offsets: a scan of the official map and a rasterized in-game rendering do5// not share a pixel grid, and neither says how to reach the other. What they6// do share is places -- fast travel points, shops, landmarks -- pinned by7// name in both. Every name both sources pin, unambiguously, is a calibration8// point, and enough of them determine the affine transformation between the9// two pixel spaces deterministically: the same captures always align the10// same way, with a residual report that says how well.11package blend1213import (14 "fmt"15 "math"16 "regexp"17 "sort"18 "strings"19)2021// Anchor is one named point in a source's own world pixels.22type Anchor struct {23 Title string24 X float6425 Y float6426}2728// Affine carries base = A·donor + t, spelled per axis:29// x' = AX·x + BX·y + CX, and y' likewise. The cross terms carry any small30// rotation or shear between the renderings; north-up sources leave them31// near zero.32type Affine struct {33 AX, BX, CX float6434 AY, BY, CY float6435}3637// Apply sends a donor-space point into base space.38func (a Affine) Apply(x, y float64) (float64, float64) {39 return a.AX*x + a.BX*y + a.CX, a.AY*x + a.BY*y + a.CY40}4142// Invert solves the transformation the other way.43func (a Affine) Invert() (Affine, bool) {44 det := a.AX*a.BY - a.BX*a.AY45 if math.Abs(det) < 1e-12 {46 return Affine{}, false47 }48 inv := Affine{49 AX: a.BY / det, BX: -a.BX / det,50 AY: -a.AY / det, BY: a.AX / det,51 }52 inv.CX = -(inv.AX*a.CX + inv.BX*a.CY)53 inv.CY = -(inv.AY*a.CX + inv.BY*a.CY)54 return inv, true55}5657// Scale reports how many base pixels one donor pixel spans, averaged over58// the axes.59func (a Affine) Scale() float64 {60 sx := math.Hypot(a.AX, a.AY)61 sy := math.Hypot(a.BX, a.BY)62 return (sx + sy) / 263}6465// Report says what the fit stood on and how well it closed, in base pixels.66type Report struct {67 Anchors int68 MedianPx float6469 P90Px float6470 WorstPx float6471}7273func (r Report) String() string {74 return fmt.Sprintf("%d anchors, median %.1fpx, p90 %.1fpx, worst %.1fpx",75 r.Anchors, r.MedianPx, r.P90Px, r.WorstPx)76}7778// The fit stands only on names that are unambiguous in both sources, needs79// enough of them to mean something, trims its worst tenth twice to shed the80// pins one source simply placed sloppily, and refuses to answer at all when81// what remains still does not close: a bad alignment drawn confidently is82// worse than none.83const (84 minAnchors = 1285 trimRounds = 286 trimKeep = 0.987 maxMedianPx = 9688 maxConditioned = 1e-1289)9091// Fit finds the affine sending donor pixels onto base pixels from the places92// both name. It is pure: the same anchors always yield the same answer.93func Fit(donor, base []Anchor) (Affine, Report, error) {94 matchedDonor, matchedBase := matchByTitle(donor, base)95 if len(matchedDonor) < minAnchors {96 return Affine{}, Report{}, fmt.Errorf(97 "only %d unambiguous shared names; %d needed", len(matchedDonor), minAnchors)98 }99100 keep := make([]int, len(matchedDonor))101 for index := range keep {102 keep[index] = index103 }104 var affine Affine105 for round := 0; ; round++ {106 var ok bool107 affine, ok = solve(matchedDonor, matchedBase, keep)108 if !ok {109 return Affine{}, Report{}, fmt.Errorf("anchors are degenerate; no transformation stands")110 }111 if round == trimRounds {112 break113 }114 sort.Slice(keep, func(a, b int) bool {115 return residual(affine, matchedDonor[keep[a]], matchedBase[keep[a]]) <116 residual(affine, matchedDonor[keep[b]], matchedBase[keep[b]])117 })118 trimmed := int(float64(len(keep)) * trimKeep)119 if trimmed < minAnchors {120 trimmed = min(minAnchors, len(keep))121 }122 keep = keep[:trimmed]123 }124125 residuals := make([]float64, len(keep))126 for at, index := range keep {127 residuals[at] = residual(affine, matchedDonor[index], matchedBase[index])128 }129 sort.Float64s(residuals)130 report := Report{131 Anchors: len(keep),132 MedianPx: residuals[len(residuals)/2],133 P90Px: residuals[int(float64(len(residuals))*0.9)],134 WorstPx: residuals[len(residuals)-1],135 }136 if report.MedianPx > maxMedianPx {137 return Affine{}, report, fmt.Errorf(138 "alignment does not close: %s", report)139 }140 return affine, report, nil141}142143func residual(a Affine, donor, base Anchor) float64 {144 x, y := a.Apply(donor.X, donor.Y)145 return math.Hypot(x-base.X, y-base.Y)146}147148// matchByTitle pairs anchors by normalized name, using only names that149// appear exactly once in each source: a name two pins share could pair150// either way, and a guessed pair is a corrupted measurement.151func matchByTitle(donor, base []Anchor) ([]Anchor, []Anchor) {152 donorByName := uniqueByName(donor)153 baseByName := uniqueByName(base)154 names := make([]string, 0, len(donorByName))155 for name := range donorByName {156 if _, shared := baseByName[name]; shared {157 names = append(names, name)158 }159 }160 sort.Strings(names)161 matchedDonor := make([]Anchor, len(names))162 matchedBase := make([]Anchor, len(names))163 for at, name := range names {164 matchedDonor[at] = donorByName[name]165 matchedBase[at] = baseByName[name]166 }167 return matchedDonor, matchedBase168}169170func uniqueByName(anchors []Anchor) map[string]Anchor {171 counts := make(map[string]int, len(anchors))172 held := make(map[string]Anchor, len(anchors))173 for _, anchor := range anchors {174 name := NormalizeTitle(anchor.Title)175 if name == "" {176 continue177 }178 counts[name]++179 held[name] = anchor180 }181 for name, count := range counts {182 if count > 1 {183 delete(held, name)184 }185 }186 return held187}188189var (190 fastTravelSuffix = regexp.MustCompile(`\s*[-–—]?\s*fast travel( point)?$`)191 nonAlphanumeric = regexp.MustCompile(`[^a-z0-9]+`)192)193194// NormalizeTitle folds the spellings apart sources give one place. The fast195// travel suffix goes because one source titles the place and another titles196// the service at it, and they are the same anchor. Anything comparing pins197// across sources should compare through this, so a name means the same thing198// to the fit and to whoever consumes it.199func NormalizeTitle(title string) string {200 title = strings.ToLower(title)201 title = fastTravelSuffix.ReplaceAllString(title, "")202 return strings.TrimSpace(nonAlphanumeric.ReplaceAllString(title, " "))203}204205// solve fits one axis at a time: three unknowns against every kept anchor,206// by the normal equations, eliminated in place.207func solve(donor, base []Anchor, keep []int) (Affine, bool) {208 xRow, ok := solveAxis(donor, base, keep, func(a Anchor) float64 { return a.X })209 if !ok {210 return Affine{}, false211 }212 yRow, ok := solveAxis(donor, base, keep, func(a Anchor) float64 { return a.Y })213 if !ok {214 return Affine{}, false215 }216 return Affine{217 AX: xRow[0], BX: xRow[1], CX: xRow[2],218 AY: yRow[0], BY: yRow[1], CY: yRow[2],219 }, true220}221222func solveAxis(donor, base []Anchor, keep []int, target func(Anchor) float64) ([3]float64, bool) {223 var sxx, sxy, syy, sx, sy, n float64224 var stx, sty, st float64225 for _, index := range keep {226 x, y := donor[index].X, donor[index].Y227 t := target(base[index])228 sxx += x * x229 sxy += x * y230 syy += y * y231 sx += x232 sy += y233 n++234 stx += t * x235 sty += t * y236 st += t237 }238 matrix := [3][3]float64{{sxx, sxy, sx}, {sxy, syy, sy}, {sx, sy, n}}239 vector := [3]float64{stx, sty, st}240 for i := range 3 {241 pivot := matrix[i][i]242 if math.Abs(pivot) < maxConditioned {243 return [3]float64{}, false244 }245 for j := i + 1; j < 3; j++ {246 factor := matrix[j][i] / pivot247 for k := range 3 {248 matrix[j][k] -= factor * matrix[i][k]249 }250 vector[j] -= factor * vector[i]251 }252 }253 var out [3]float64254 for i := 2; i >= 0; i-- {255 out[i] = vector[i]256 for k := i + 1; k < 3; k++ {257 out[i] -= matrix[i][k] * out[k]258 }259 out[i] /= matrix[i][i]260 }261 return out, true262}