1package crawl23import (4 "image"5 "image/draw"6 "math"7)89// A deterministic, high-quality downsampler.10//11// The Blue Marble capture cuts a pinned source image down to the reference12// level the deriver folds, and the cut has to be the same bytes on every13// machine: the tiles' content hashes feed the derivation stamp, and the stamp14// feeds the bundle's name. The deriver's own 2×2 box fold cannot do this cut --15// the ratio is not a power of two -- so this file is the one other resampler in16// the lane, and it is built for the same property the fold has: integer17// arithmetic wherever a pixel is decided.18//19// The kernel is Catmull-Rom, widened by the scale factor so a downsample20// averages over the ground a destination pixel covers instead of aliasing21// across it. Tap weights are computed once per output column and row -- the22// only floating-point step -- with every operation explicitly rounded, so the23// compiler may not fuse two of them into an differently-rounded one, and are24// then frozen to 15-bit fixed point, normalized to sum exactly to one.25// Everything per-pixel after that is integer multiply-accumulate.2627// resampleShift is the fixed-point precision a frozen weight carries.28const resampleShift = 152930// resampleTaps is the weights one output index takes over the source axis.31type resampleTaps struct {32 first int33 weights []int3234}3536// resample scales src to width×height through the widened Catmull-Rom kernel,37// separably: columns first into a fixed-point intermediate, then rows.38func resample(src *image.NRGBA, width, height int) *image.NRGBA {39 bounds := src.Bounds()40 columns := resampleAxis(bounds.Dx(), width)41 rows := resampleAxis(bounds.Dy(), height)4243 out := image.NewNRGBA(image.Rect(0, 0, width, height))44 // Each output row needs only a window of horizontally-resampled source45 // rows, and output rows walk downward, so the intermediates live in a46 // small cache that trails the walk instead of a whole-image buffer.47 cache := make(map[int][]int32)48 for y := range height {49 taps := rows[y]50 for r := range taps.weights {51 row := taps.first + r52 if _, held := cache[row]; !held {53 cache[row] = resampleRow(src, row, columns, width)54 }55 }56 for column := range width {57 var acc [4]int6458 for r, weight := range taps.weights {59 held := cache[taps.first+r]60 for channel := range 4 {61 acc[channel] += int64(weight) * int64(held[column*4+channel])62 }63 }64 at := out.PixOffset(column, y)65 for channel := range 4 {66 out.Pix[at+channel] = clamp8(acc[channel])67 }68 }69 for row := range cache {70 if row < taps.first {71 delete(cache, row)72 }73 }74 }75 return out76}7778// resampleRow runs the horizontal pass over one source row, answering each79// output column's channels in 15-bit fixed point.80func resampleRow(src *image.NRGBA, row int, columns []resampleTaps, width int) []int32 {81 bounds := src.Bounds()82 out := make([]int32, width*4)83 base := src.PixOffset(bounds.Min.X, bounds.Min.Y+row)84 line := src.Pix[base : base+bounds.Dx()*4]85 for column := range width {86 taps := columns[column]87 var acc [4]int3288 for t, weight := range taps.weights {89 at := (taps.first + t) * 490 for channel := range 4 {91 acc[channel] += weight * int32(line[at+channel])92 }93 }94 copy(out[column*4:], acc[:])95 }96 return out97}9899// clamp8 rounds a doubly-scaled accumulator back to a byte, half away from100// zero, clamped to the channel's range.101func clamp8(value int64) uint8 {102 value = (value + 1<<(2*resampleShift-1)) >> (2 * resampleShift)103 if value < 0 {104 return 0105 }106 if value > 255 {107 return 255108 }109 return uint8(value)110}111112// resampleAxis freezes the tap weights for one axis. Taps are clamped inside113// the source, folding an edge overhang onto the edge sample, which extends the114// picture rather than darkening its border.115func resampleAxis(src, dst int) []resampleTaps {116 scale := float64(src) / float64(dst)117 if scale < 1 {118 scale = 1119 }120 support := float64(2 * scale)121 out := make([]resampleTaps, dst)122 for i := range dst {123 // The output pixel's centre on the source axis.124 center := float64(float64(2*i+1)*float64(src)) / float64(2*dst)125 first := int(math.Ceil(float64(center - support - 0.5)))126 last := int(math.Floor(float64(center + support - 0.5)))127 raw := make([]float64, 0, last-first+1)128 sum := 0.0129 for j := first; j <= last; j++ {130 t := float64(float64(float64(j)+0.5-center) / scale)131 weight := catmullRom(t)132 raw = append(raw, weight)133 sum = float64(sum + weight)134 }135 frozen := freezeWeights(raw, sum)136 // Clamp the window inside the source, folding out-of-range taps onto137 // the edge sample they would have read past.138 taps := resampleTaps{first: first, weights: frozen}139 out[i] = clampTaps(taps, src)140 }141 return out142}143144// catmullRom is the kernel, spelled so every floating-point operation rounds on145// its own: an explicit float64 conversion is a rounding point the compiler must146// honour, which keeps a fused multiply-add from computing a differently-rounded147// weight on one architecture than on another.148func catmullRom(t float64) float64 {149 if t < 0 {150 t = -t151 }152 switch {153 case t < 1:154 // 1.5t³ − 2.5t² + 1155 return float64(float64(float64(float64(1.5*t)-2.5)*float64(t*t)) + 1)156 case t < 2:157 // −0.5t³ + 2.5t² − 4t + 2158 return float64(float64(float64(float64(float64(-0.5*t)+2.5)*t)-4)*t + 2)159 default:160 return 0161 }162}163164// freezeWeights quantizes one window of weights to fixed point summing exactly165// to one, handing the rounding residual to the heaviest tap.166func freezeWeights(raw []float64, sum float64) []int32 {167 frozen := make([]int32, len(raw))168 total := int32(0)169 heaviest := 0170 for i, weight := range raw {171 frozen[i] = int32(math.RoundToEven(float64(weight/sum) * (1 << resampleShift)))172 total += frozen[i]173 if frozen[i] > frozen[heaviest] {174 heaviest = i175 }176 }177 frozen[heaviest] += 1<<resampleShift - total178 return frozen179}180181// clampTaps folds taps that fall outside the source onto its first and last182// samples.183func clampTaps(taps resampleTaps, src int) resampleTaps {184 out := resampleTaps{first: taps.first}185 if out.first < 0 {186 out.first = 0187 }188 last := min(taps.first+len(taps.weights)-1, src-1)189 if out.first > last {190 out.first, last = 0, 0191 }192 out.weights = make([]int32, last-out.first+1)193 for i, weight := range taps.weights {194 at := min(max(taps.first+i, out.first), last)195 out.weights[at-out.first] += weight196 }197 return out198}199200// toNRGBA lands a decoded image in the one layout the resampler reads. The201// conversion is the standard library's own integer arithmetic, so a JPEG202// decodes to the same bytes on every machine.203func toNRGBA(src image.Image) *image.NRGBA {204 if held, ok := src.(*image.NRGBA); ok {205 return held206 }207 out := image.NewNRGBA(image.Rect(0, 0, src.Bounds().Dx(), src.Bounds().Dy()))208 draw.Draw(out, out.Bounds(), src, src.Bounds().Min, draw.Src)209 return out210}