1 // Copyright 2016 The Go Authors. All rights reserved. 2 // Use of this source code is governed by a BSD-style 3 // license that can be found in the LICENSE file. 4 5 package vector 6 7 // This file contains a fixed point math implementation of the vector 8 // graphics rasterizer. 9 10 const ( 11 // ϕ is the number of binary digits after the fixed point. 12 // 13 // For example, if ϕ == 10 (and int1ϕ is based on the int32 type) then we 14 // are using 22.10 fixed point math. 15 // 16 // When changing this number, also change the assembly code (search for ϕ 17 // in the .s files). 18 ϕ = 9 19 20 fxOne int1ϕ = 1 << ϕ 21 fxOneAndAHalf int1ϕ = 1<<ϕ + 1<<(ϕ-1) 22 fxOneMinusIota int1ϕ = 1<<ϕ - 1 // Used for rounding up. 23 ) 24 25 // int1ϕ is a signed fixed-point number with 1*ϕ binary digits after the fixed 26 // point. 27 type int1ϕ int32 28 29 // int2ϕ is a signed fixed-point number with 2*ϕ binary digits after the fixed 30 // point. 31 // 32 // The Rasterizer's bufU32 field, nominally of type []uint32 (since that slice 33 // is also used by other code), can be thought of as a []int2ϕ during the 34 // fixedLineTo method. Lines of code that are actually like: 35 // 36 // buf[i] += uint32(etc) // buf has type []uint32. 37 // 38 // can be thought of as 39 // 40 // buf[i] += int2ϕ(etc) // buf has type []int2ϕ. 41 type int2ϕ int32 42 43 func fixedMax(x, y int1ϕ) int1ϕ { 44 if x > y { 45 return x 46 } 47 return y 48 } 49 50 func fixedMin(x, y int1ϕ) int1ϕ { 51 if x < y { 52 return x 53 } 54 return y 55 } 56 57 func fixedFloor(x int1ϕ) int32 { return int32(x >> ϕ) } 58 func fixedCeil(x int1ϕ) int32 { return int32((x + fxOneMinusIota) >> ϕ) } 59 60 func (z *Rasterizer) fixedLineTo(bx, by float32) { 61 ax, ay := z.penX, z.penY 62 z.penX, z.penY = bx, by 63 dir := int1ϕ(1) 64 if ay > by { 65 dir, ax, ay, bx, by = -1, bx, by, ax, ay 66 } 67 // Horizontal line segments yield no change in coverage. Almost horizontal 68 // segments would yield some change, in ideal math, but the computation 69 // further below, involving 1 / (by - ay), is unstable in fixed point math, 70 // so we treat the segment as if it was perfectly horizontal. 71 if by-ay <= 0.000001 { 72 return 73 } 74 dxdy := (bx - ax) / (by - ay) 75 76 ayϕ := int1ϕ(ay * float32(fxOne)) 77 byϕ := int1ϕ(by * float32(fxOne)) 78 79 x := int1ϕ(ax * float32(fxOne)) 80 y := fixedFloor(ayϕ) 81 yMax := fixedCeil(byϕ) 82 if yMax > int32(z.size.Y) { 83 yMax = int32(z.size.Y) 84 } 85 width := int32(z.size.X) 86 87 for ; y < yMax; y++ { 88 dy := fixedMin(int1ϕ(y+1)<<ϕ, byϕ) - fixedMax(int1ϕ(y)<<ϕ, ayϕ) 89 xNext := x + int1ϕ(float32(dy)*dxdy) 90 if y < 0 { 91 x = xNext 92 continue 93 } 94 buf := z.bufU32[y*width:] 95 d := dy * dir // d ranges up to ±1<<(1*ϕ). 96 x0, x1 := x, xNext 97 if x > xNext { 98 x0, x1 = x1, x0 99 } 100 x0i := fixedFloor(x0) 101 x0Floor := int1ϕ(x0i) << ϕ 102 x1i := fixedCeil(x1) 103 x1Ceil := int1ϕ(x1i) << ϕ 104 105 if x1i <= x0i+1 { 106 xmf := (x+xNext)>>1 - x0Floor 107 if i := clamp(x0i+0, width); i < uint(len(buf)) { 108 buf[i] += uint32(d * (fxOne - xmf)) 109 } 110 if i := clamp(x0i+1, width); i < uint(len(buf)) { 111 buf[i] += uint32(d * xmf) 112 } 113 } else { 114 oneOverS := x1 - x0 115 twoOverS := 2 * oneOverS 116 x0f := x0 - x0Floor 117 oneMinusX0f := fxOne - x0f 118 oneMinusX0fSquared := oneMinusX0f * oneMinusX0f 119 x1f := x1 - x1Ceil + fxOne 120 x1fSquared := x1f * x1f 121 122 // These next two variables are unused, as rounding errors are 123 // minimized when we delay the division by oneOverS for as long as 124 // possible. These lines of code (and the "In ideal math" comments 125 // below) are commented out instead of deleted in order to aid the 126 // comparison with the floating point version of the rasterizer. 127 // 128 // a0 := ((oneMinusX0f * oneMinusX0f) >> 1) / oneOverS 129 // am := ((x1f * x1f) >> 1) / oneOverS 130 131 if i := clamp(x0i, width); i < uint(len(buf)) { 132 // In ideal math: buf[i] += uint32(d * a0) 133 D := oneMinusX0fSquared // D ranges up to ±1<<(2*ϕ). 134 D *= d // D ranges up to ±1<<(3*ϕ). 135 D /= twoOverS 136 buf[i] += uint32(D) 137 } 138 139 if x1i == x0i+2 { 140 if i := clamp(x0i+1, width); i < uint(len(buf)) { 141 // In ideal math: buf[i] += uint32(d * (fxOne - a0 - am)) 142 // 143 // (x1i == x0i+2) and (twoOverS == 2 * (x1 - x0)) implies 144 // that twoOverS ranges up to +1<<(1*ϕ+2). 145 D := twoOverS<<ϕ - oneMinusX0fSquared - x1fSquared // D ranges up to ±1<<(2*ϕ+2). 146 D *= d // D ranges up to ±1<<(3*ϕ+2). 147 D /= twoOverS 148 buf[i] += uint32(D) 149 } 150 } else { 151 // This is commented out for the same reason as a0 and am. 152 // 153 // a1 := ((fxOneAndAHalf - x0f) << ϕ) / oneOverS 154 155 if i := clamp(x0i+1, width); i < uint(len(buf)) { 156 // In ideal math: 157 // buf[i] += uint32(d * (a1 - a0)) 158 // or equivalently (but better in non-ideal, integer math, 159 // with respect to rounding errors), 160 // buf[i] += uint32(A * d / twoOverS) 161 // where 162 // A = (a1 - a0) * twoOverS 163 // = a1*twoOverS - a0*twoOverS 164 // Noting that twoOverS/oneOverS equals 2, substituting for 165 // a0 and then a1, given above, yields: 166 // A = a1*twoOverS - oneMinusX0fSquared 167 // = (fxOneAndAHalf-x0f)<<(ϕ+1) - oneMinusX0fSquared 168 // = fxOneAndAHalf<<(ϕ+1) - x0f<<(ϕ+1) - oneMinusX0fSquared 169 // 170 // This is a positive number minus two non-negative 171 // numbers. For an upper bound on A, the positive number is 172 // P = fxOneAndAHalf<<(ϕ+1) 173 // < (2*fxOne)<<(ϕ+1) 174 // = fxOne<<(ϕ+2) 175 // = 1<<(2*ϕ+2) 176 // 177 // For a lower bound on A, the two non-negative numbers are 178 // N = x0f<<(ϕ+1) + oneMinusX0fSquared 179 // ≤ x0f<<(ϕ+1) + fxOne*fxOne 180 // = x0f<<(ϕ+1) + 1<<(2*ϕ) 181 // < x0f<<(ϕ+1) + 1<<(2*ϕ+1) 182 // ≤ fxOne<<(ϕ+1) + 1<<(2*ϕ+1) 183 // = 1<<(2*ϕ+1) + 1<<(2*ϕ+1) 184 // = 1<<(2*ϕ+2) 185 // 186 // Thus, A ranges up to ±1<<(2*ϕ+2). It is possible to 187 // derive a tighter bound, but this bound is sufficient to 188 // reason about overflow. 189 D := (fxOneAndAHalf-x0f)<<(ϕ+1) - oneMinusX0fSquared // D ranges up to ±1<<(2*ϕ+2). 190 D *= d // D ranges up to ±1<<(3*ϕ+2). 191 D /= twoOverS 192 buf[i] += uint32(D) 193 } 194 dTimesS := uint32((d << (2 * ϕ)) / oneOverS) 195 for xi := x0i + 2; xi < x1i-1; xi++ { 196 if i := clamp(xi, width); i < uint(len(buf)) { 197 buf[i] += dTimesS 198 } 199 } 200 201 // This is commented out for the same reason as a0 and am. 202 // 203 // a2 := a1 + (int1ϕ(x1i-x0i-3)<<(2*ϕ))/oneOverS 204 205 if i := clamp(x1i-1, width); i < uint(len(buf)) { 206 // In ideal math: 207 // buf[i] += uint32(d * (fxOne - a2 - am)) 208 // or equivalently (but better in non-ideal, integer math, 209 // with respect to rounding errors), 210 // buf[i] += uint32(A * d / twoOverS) 211 // where 212 // A = (fxOne - a2 - am) * twoOverS 213 // = twoOverS<<ϕ - a2*twoOverS - am*twoOverS 214 // Noting that twoOverS/oneOverS equals 2, substituting for 215 // am and then a2, given above, yields: 216 // A = twoOverS<<ϕ - a2*twoOverS - x1f*x1f 217 // = twoOverS<<ϕ - a1*twoOverS - (int1ϕ(x1i-x0i-3)<<(2*ϕ))*2 - x1f*x1f 218 // = twoOverS<<ϕ - a1*twoOverS - int1ϕ(x1i-x0i-3)<<(2*ϕ+1) - x1f*x1f 219 // Substituting for a1, given above, yields: 220 // A = twoOverS<<ϕ - ((fxOneAndAHalf-x0f)<<ϕ)*2 - int1ϕ(x1i-x0i-3)<<(2*ϕ+1) - x1f*x1f 221 // = twoOverS<<ϕ - (fxOneAndAHalf-x0f)<<(ϕ+1) - int1ϕ(x1i-x0i-3)<<(2*ϕ+1) - x1f*x1f 222 // = B<<ϕ - x1f*x1f 223 // where 224 // B = twoOverS - (fxOneAndAHalf-x0f)<<1 - int1ϕ(x1i-x0i-3)<<(ϕ+1) 225 // = (x1-x0)<<1 - (fxOneAndAHalf-x0f)<<1 - int1ϕ(x1i-x0i-3)<<(ϕ+1) 226 // 227 // Re-arranging the defintions given above: 228 // x0Floor := int1ϕ(x0i) << ϕ 229 // x0f := x0 - x0Floor 230 // x1Ceil := int1ϕ(x1i) << ϕ 231 // x1f := x1 - x1Ceil + fxOne 232 // combined with fxOne = 1<<ϕ yields: 233 // x0 = x0f + int1ϕ(x0i)<<ϕ 234 // x1 = x1f + int1ϕ(x1i-1)<<ϕ 235 // so that expanding (x1-x0) yields: 236 // B = (x1f-x0f + int1ϕ(x1i-x0i-1)<<ϕ)<<1 - (fxOneAndAHalf-x0f)<<1 - int1ϕ(x1i-x0i-3)<<(ϕ+1) 237 // = (x1f-x0f)<<1 + int1ϕ(x1i-x0i-1)<<(ϕ+1) - (fxOneAndAHalf-x0f)<<1 - int1ϕ(x1i-x0i-3)<<(ϕ+1) 238 // A large part of the second and fourth terms cancel: 239 // B = (x1f-x0f)<<1 - (fxOneAndAHalf-x0f)<<1 - int1ϕ(-2)<<(ϕ+1) 240 // = (x1f-x0f)<<1 - (fxOneAndAHalf-x0f)<<1 + 1<<(ϕ+2) 241 // = (x1f - fxOneAndAHalf)<<1 + 1<<(ϕ+2) 242 // The first term, (x1f - fxOneAndAHalf)<<1, is a negative 243 // number, bounded below by -fxOneAndAHalf<<1, which is 244 // greater than -fxOne<<2, or -1<<(ϕ+2). Thus, B ranges up 245 // to ±1<<(ϕ+2). One final simplification: 246 // B = x1f<<1 + (1<<(ϕ+2) - fxOneAndAHalf<<1) 247 const C = 1<<(ϕ+2) - fxOneAndAHalf<<1 248 D := x1f<<1 + C // D ranges up to ±1<<(1*ϕ+2). 249 D <<= ϕ // D ranges up to ±1<<(2*ϕ+2). 250 D -= x1fSquared // D ranges up to ±1<<(2*ϕ+3). 251 D *= d // D ranges up to ±1<<(3*ϕ+3). 252 D /= twoOverS 253 buf[i] += uint32(D) 254 } 255 } 256 257 if i := clamp(x1i, width); i < uint(len(buf)) { 258 // In ideal math: buf[i] += uint32(d * am) 259 D := x1fSquared // D ranges up to ±1<<(2*ϕ). 260 D *= d // D ranges up to ±1<<(3*ϕ). 261 D /= twoOverS 262 buf[i] += uint32(D) 263 } 264 } 265 266 x = xNext 267 } 268 } 269 270 func fixedAccumulateOpOver(dst []uint8, src []uint32) { 271 // Sanity check that len(dst) >= len(src). 272 if len(dst) < len(src) { 273 return 274 } 275 276 acc := int2ϕ(0) 277 for i, v := range src { 278 acc += int2ϕ(v) 279 a := acc 280 if a < 0 { 281 a = -a 282 } 283 a >>= 2*ϕ - 16 284 if a > 0xffff { 285 a = 0xffff 286 } 287 // This algorithm comes from the standard library's image/draw package. 288 dstA := uint32(dst[i]) * 0x101 289 maskA := uint32(a) 290 outA := dstA*(0xffff-maskA)/0xffff + maskA 291 dst[i] = uint8(outA >> 8) 292 } 293 } 294 295 func fixedAccumulateOpSrc(dst []uint8, src []uint32) { 296 // Sanity check that len(dst) >= len(src). 297 if len(dst) < len(src) { 298 return 299 } 300 301 acc := int2ϕ(0) 302 for i, v := range src { 303 acc += int2ϕ(v) 304 a := acc 305 if a < 0 { 306 a = -a 307 } 308 a >>= 2*ϕ - 8 309 if a > 0xff { 310 a = 0xff 311 } 312 dst[i] = uint8(a) 313 } 314 } 315 316 func fixedAccumulateMask(buf []uint32) { 317 acc := int2ϕ(0) 318 for i, v := range buf { 319 acc += int2ϕ(v) 320 a := acc 321 if a < 0 { 322 a = -a 323 } 324 a >>= 2*ϕ - 16 325 if a > 0xffff { 326 a = 0xffff 327 } 328 buf[i] = uint32(a) 329 } 330 } 331