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
5package vector
6
7// This file contains a fixed point math implementation of the vector
8// graphics rasterizer.
9
10const (
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.
27type 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ϕ.
41type int2ϕ int32
42
43func fixedMax(x, y int1ϕ) int1ϕ {
44 if x > y {
45 return x
46 }
47 return y
48}
49
50func fixedMin(x, y int1ϕ) int1ϕ {
51 if x < y {
52 return x
53 }
54 return y
55}
56
57func fixedFloor(x int1ϕ) int32 { return int32(x >> ϕ) }
58func fixedCeil(x int1ϕ) int32 { return int32((x + fxOneMinusIota) >> ϕ) }
59
60func (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
270func 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
295func 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
316func 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}