1/*
2 * Copyright (c) 2020 - 2023 the ThorVG project. All rights reserved.
3
4 * Permission is hereby granted, free of charge, to any person obtaining a copy
5 * of this software and associated documentation files (the "Software"), to deal
6 * in the Software without restriction, including without limitation the rights
7 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
8 * copies of the Software, and to permit persons to whom the Software is
9 * furnished to do so, subject to the following conditions:
10
11 * The above copyright notice and this permission notice shall be included in all
12 * copies or substantial portions of the Software.
13
14 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
15 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
16 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
17 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
18 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
19 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
20 * SOFTWARE.
21 */
22
23#include <math.h>
24#include "tvgSwCommon.h"
25
26
27/************************************************************************/
28/* Internal Class Implementation */
29/************************************************************************/
30
31//clz: count leading zero’s
32#if defined(_MSC_VER) && !defined(__clang__)
33 #include <intrin.h>
34 static uint32_t __inline _clz(uint32_t value)
35 {
36 unsigned long leadingZero = 0;
37 if (_BitScanReverse(&leadingZero, value)) return 31 - leadingZero;
38 else return 32;
39 }
40#else
41 #define _clz(x) __builtin_clz((x))
42#endif
43
44
45constexpr SwFixed CORDIC_FACTOR = 0xDBD95B16UL; //the Cordic shrink factor 0.858785336480436 * 2^32
46
47//this table was generated for SW_FT_PI = 180L << 16, i.e. degrees
48constexpr static auto ATAN_MAX = 23;
49constexpr static SwFixed ATAN_TBL[] = {
50 1740967L, 919879L, 466945L, 234379L, 117304L, 58666L, 29335L,
51 14668L, 7334L, 3667L, 1833L, 917L, 458L, 229L, 115L,
52 57L, 29L, 14L, 7L, 4L, 2L, 1L};
53
54static inline SwCoord SATURATE(const SwCoord x)
55{
56 return (x >> (sizeof(SwCoord) * 8 - 1));
57}
58
59
60static inline SwFixed PAD_ROUND(const SwFixed x, int32_t n)
61{
62 return (((x) + ((n)/2)) & ~((n)-1));
63}
64
65
66static SwCoord _downscale(SwFixed x)
67{
68 //multiply a give value by the CORDIC shrink factor
69 auto s = abs(x);
70 int64_t t = (s * static_cast<int64_t>(CORDIC_FACTOR)) + 0x100000000UL;
71 s = static_cast<SwFixed>(t >> 32);
72 if (x < 0) s = -s;
73 return s;
74}
75
76
77static int32_t _normalize(SwPoint& pt)
78{
79 /* the highest bit in overflow-safe vector components
80 MSB of 0.858785336480436 * sqrt(0.5) * 2^30 */
81 constexpr auto SAFE_MSB = 29;
82
83 auto v = pt;
84
85 //High order bit(MSB)
86 int32_t shift = 31 - _clz(abs(v.x) | abs(v.y));
87
88 if (shift <= SAFE_MSB) {
89 shift = SAFE_MSB - shift;
90 pt.x = static_cast<SwCoord>((unsigned long)v.x << shift);
91 pt.y = static_cast<SwCoord>((unsigned long)v.y << shift);
92 } else {
93 shift -= SAFE_MSB;
94 pt.x = v.x >> shift;
95 pt.y = v.y >> shift;
96 shift = -shift;
97 }
98 return shift;
99}
100
101
102static void _polarize(SwPoint& pt)
103{
104 auto v = pt;
105 SwFixed theta;
106
107 //Get the vector into [-PI/4, PI/4] sector
108 if (v.y > v.x) {
109 if (v.y > -v.x) {
110 auto tmp = v.y;
111 v.y = -v.x;
112 v.x = tmp;
113 theta = SW_ANGLE_PI2;
114 } else {
115 theta = v.y > 0 ? SW_ANGLE_PI : -SW_ANGLE_PI;
116 v.x = -v.x;
117 v.y = -v.y;
118 }
119 } else {
120 if (v.y < -v.x) {
121 theta = -SW_ANGLE_PI2;
122 auto tmp = -v.y;
123 v.y = v.x;
124 v.x = tmp;
125 } else {
126 theta = 0;
127 }
128 }
129
130 auto atan = ATAN_TBL;
131 uint32_t i;
132 SwFixed j;
133
134 //Pseudorotations. with right shifts
135 for (i = 1, j = 1; i < ATAN_MAX; j <<= 1, ++i) {
136 if (v.y > 0) {
137 auto tmp = v.x + ((v.y + j) >> i);
138 v.y = v.y - ((v.x + j) >> i);
139 v.x = tmp;
140 theta += *atan++;
141 } else {
142 auto tmp = v.x - ((v.y + j) >> i);
143 v.y = v.y + ((v.x + j) >> i);
144 v.x = tmp;
145 theta -= *atan++;
146 }
147 }
148
149 //round theta
150 if (theta >= 0) theta = PAD_ROUND(theta, 32);
151 else theta = -PAD_ROUND(-theta, 32);
152
153 pt.x = v.x;
154 pt.y = theta;
155}
156
157
158static void _rotate(SwPoint& pt, SwFixed theta)
159{
160 SwFixed x = pt.x;
161 SwFixed y = pt.y;
162
163 //Rotate inside [-PI/4, PI/4] sector
164 while (theta < -SW_ANGLE_PI4) {
165 auto tmp = y;
166 y = -x;
167 x = tmp;
168 theta += SW_ANGLE_PI2;
169 }
170
171 while (theta > SW_ANGLE_PI4) {
172 auto tmp = -y;
173 y = x;
174 x = tmp;
175 theta -= SW_ANGLE_PI2;
176 }
177
178 auto atan = ATAN_TBL;
179 uint32_t i;
180 SwFixed j;
181
182 for (i = 1, j = 1; i < ATAN_MAX; j <<= 1, ++i) {
183 if (theta < 0) {
184 auto tmp = x + ((y + j) >> i);
185 y = y - ((x + j) >> i);
186 x = tmp;
187 theta += *atan++;
188 } else {
189 auto tmp = x - ((y + j) >> i);
190 y = y + ((x + j) >> i);
191 x = tmp;
192 theta -= *atan++;
193 }
194 }
195
196 pt.x = static_cast<SwCoord>(x);
197 pt.y = static_cast<SwCoord>(y);
198}
199
200
201/************************************************************************/
202/* External Class Implementation */
203/************************************************************************/
204
205SwFixed mathMean(SwFixed angle1, SwFixed angle2)
206{
207 return angle1 + mathDiff(angle1, angle2) / 2;
208}
209
210
211bool mathSmallCubic(const SwPoint* base, SwFixed& angleIn, SwFixed& angleMid, SwFixed& angleOut)
212{
213 auto d1 = base[2] - base[3];
214 auto d2 = base[1] - base[2];
215 auto d3 = base[0] - base[1];
216
217 if (d1.small()) {
218 if (d2.small()) {
219 if (d3.small()) {
220 //basically a point.
221 //do nothing to retain original direction
222 } else {
223 angleIn = angleMid = angleOut = mathAtan(d3);
224 }
225 } else {
226 if (d3.small()) {
227 angleIn = angleMid = angleOut = mathAtan(d2);
228 } else {
229 angleIn = angleMid = mathAtan(d2);
230 angleOut = mathAtan(d3);
231 }
232 }
233 } else {
234 if (d2.small()) {
235 if (d3.small()) {
236 angleIn = angleMid = angleOut = mathAtan(d1);
237 } else {
238 angleIn = mathAtan(d1);
239 angleOut = mathAtan(d3);
240 angleMid = mathMean(angleIn, angleOut);
241 }
242 } else {
243 if (d3.small()) {
244 angleIn = mathAtan(d1);
245 angleMid = angleOut = mathAtan(d2);
246 } else {
247 angleIn = mathAtan(d1);
248 angleMid = mathAtan(d2);
249 angleOut = mathAtan(d3);
250 }
251 }
252 }
253
254 auto theta1 = abs(mathDiff(angleIn, angleMid));
255 auto theta2 = abs(mathDiff(angleMid, angleOut));
256
257 if ((theta1 < (SW_ANGLE_PI / 8)) && (theta2 < (SW_ANGLE_PI / 8))) return true;
258 return false;
259}
260
261
262int64_t mathMultiply(int64_t a, int64_t b)
263{
264 int32_t s = 1;
265
266 //move sign
267 if (a < 0) {
268 a = -a;
269 s = -s;
270 }
271 if (b < 0) {
272 b = -b;
273 s = -s;
274 }
275 int64_t c = (a * b + 0x8000L) >> 16;
276 return (s > 0) ? c : -c;
277}
278
279
280int64_t mathDivide(int64_t a, int64_t b)
281{
282 int32_t s = 1;
283
284 //move sign
285 if (a < 0) {
286 a = -a;
287 s = -s;
288 }
289 if (b < 0) {
290 b = -b;
291 s = -s;
292 }
293 int64_t q = b > 0 ? ((a << 16) + (b >> 1)) / b : 0x7FFFFFFFL;
294 return (s < 0 ? -q : q);
295}
296
297
298int64_t mathMulDiv(int64_t a, int64_t b, int64_t c)
299{
300 int32_t s = 1;
301
302 //move sign
303 if (a < 0) {
304 a = -a;
305 s = -s;
306 }
307 if (b < 0) {
308 b = -b;
309 s = -s;
310 }
311 if (c < 0) {
312 c = -c;
313 s = -s;
314 }
315 int64_t d = c > 0 ? (a * b + (c >> 1)) / c : 0x7FFFFFFFL;
316
317 return (s > 0 ? d : -d);
318}
319
320
321void mathRotate(SwPoint& pt, SwFixed angle)
322{
323 if (angle == 0 || (pt.x == 0 && pt.y == 0)) return;
324
325 auto v = pt;
326 auto shift = _normalize(v);
327
328 auto theta = angle;
329 _rotate(v, theta);
330
331 v.x = _downscale(v.x);
332 v.y = _downscale(v.y);
333
334 if (shift > 0) {
335 auto half = static_cast<int32_t>(1L << (shift - 1));
336 pt.x = (v.x + half + SATURATE(v.x)) >> shift;
337 pt.y = (v.y + half + SATURATE(v.y)) >> shift;
338 } else {
339 shift = -shift;
340 pt.x = static_cast<SwCoord>((unsigned long)v.x << shift);
341 pt.y = static_cast<SwCoord>((unsigned long)v.y << shift);
342 }
343}
344
345SwFixed mathTan(SwFixed angle)
346{
347 SwPoint v = {CORDIC_FACTOR >> 8, 0};
348 _rotate(v, angle);
349 return mathDivide(v.y, v.x);
350}
351
352
353SwFixed mathAtan(const SwPoint& pt)
354{
355 if (pt.x == 0 && pt.y == 0) return 0;
356
357 auto v = pt;
358 _normalize(v);
359 _polarize(v);
360
361 return v.y;
362}
363
364
365SwFixed mathSin(SwFixed angle)
366{
367 return mathCos(SW_ANGLE_PI2 - angle);
368}
369
370
371SwFixed mathCos(SwFixed angle)
372{
373 SwPoint v = {CORDIC_FACTOR >> 8, 0};
374 _rotate(v, angle);
375 return (v.x + 0x80L) >> 8;
376}
377
378
379SwFixed mathLength(const SwPoint& pt)
380{
381 auto v = pt;
382
383 //trivial case
384 if (v.x == 0) return abs(v.y);
385 if (v.y == 0) return abs(v.x);
386
387 //general case
388 auto shift = _normalize(v);
389 _polarize(v);
390 v.x = _downscale(v.x);
391
392 if (shift > 0) return (v.x + (static_cast<SwFixed>(1) << (shift -1))) >> shift;
393 return static_cast<SwFixed>((uint32_t)v.x << -shift);
394}
395
396
397void mathSplitCubic(SwPoint* base)
398{
399 SwCoord a, b, c, d;
400
401 base[6].x = base[3].x;
402 c = base[1].x;
403 d = base[2].x;
404 base[1].x = a = (base[0].x + c) / 2;
405 base[5].x = b = (base[3].x + d) / 2;
406 c = (c + d) / 2;
407 base[2].x = a = (a + c) / 2;
408 base[4].x = b = (b + c) / 2;
409 base[3].x = (a + b) / 2;
410
411 base[6].y = base[3].y;
412 c = base[1].y;
413 d = base[2].y;
414 base[1].y = a = (base[0].y + c) / 2;
415 base[5].y = b = (base[3].y + d) / 2;
416 c = (c + d) / 2;
417 base[2].y = a = (a + c) / 2;
418 base[4].y = b = (b + c) / 2;
419 base[3].y = (a + b) / 2;
420}
421
422
423SwFixed mathDiff(SwFixed angle1, SwFixed angle2)
424{
425 auto delta = angle2 - angle1;
426
427 delta %= SW_ANGLE_2PI;
428 if (delta < 0) delta += SW_ANGLE_2PI;
429 if (delta > SW_ANGLE_PI) delta -= SW_ANGLE_2PI;
430
431 return delta;
432}
433
434
435SwPoint mathTransform(const Point* to, const Matrix* transform)
436{
437 if (!transform) return {TO_SWCOORD(to->x), TO_SWCOORD(to->y)};
438
439 auto tx = to->x * transform->e11 + to->y * transform->e12 + transform->e13;
440 auto ty = to->x * transform->e21 + to->y * transform->e22 + transform->e23;
441
442 return {TO_SWCOORD(tx), TO_SWCOORD(ty)};
443}
444
445
446bool mathClipBBox(const SwBBox& clipper, SwBBox& clipee)
447{
448 clipee.max.x = (clipee.max.x < clipper.max.x) ? clipee.max.x : clipper.max.x;
449 clipee.max.y = (clipee.max.y < clipper.max.y) ? clipee.max.y : clipper.max.y;
450 clipee.min.x = (clipee.min.x > clipper.min.x) ? clipee.min.x : clipper.min.x;
451 clipee.min.y = (clipee.min.y > clipper.min.y) ? clipee.min.y : clipper.min.y;
452
453 //Check valid region
454 if (clipee.max.x - clipee.min.x < 1 && clipee.max.y - clipee.min.y < 1) return false;
455
456 //Check boundary
457 if (clipee.min.x >= clipper.max.x || clipee.min.y >= clipper.max.y ||
458 clipee.max.x <= clipper.min.x || clipee.max.y <= clipper.min.y) return false;
459
460 return true;
461}
462
463
464bool mathUpdateOutlineBBox(const SwOutline* outline, const SwBBox& clipRegion, SwBBox& renderRegion, bool fastTrack)
465{
466 if (!outline) return false;
467
468 auto pt = outline->pts.data;
469
470 if (outline->pts.empty() || outline->cntrs.empty()) {
471 renderRegion.reset();
472 return false;
473 }
474
475 auto xMin = pt->x;
476 auto xMax = pt->x;
477 auto yMin = pt->y;
478 auto yMax = pt->y;
479
480 for (++pt; pt < outline->pts.end(); ++pt) {
481 if (xMin > pt->x) xMin = pt->x;
482 if (xMax < pt->x) xMax = pt->x;
483 if (yMin > pt->y) yMin = pt->y;
484 if (yMax < pt->y) yMax = pt->y;
485 }
486 //Since no antialiasing is applied in the Fast Track case,
487 //the rasterization region has to be rearranged.
488 //https://github.com/Samsung/thorvg/issues/916
489 if (fastTrack) {
490 renderRegion.min.x = static_cast<SwCoord>(round(xMin / 64.0f));
491 renderRegion.max.x = static_cast<SwCoord>(round(xMax / 64.0f));
492 renderRegion.min.y = static_cast<SwCoord>(round(yMin / 64.0f));
493 renderRegion.max.y = static_cast<SwCoord>(round(yMax / 64.0f));
494 } else {
495 renderRegion.min.x = xMin >> 6;
496 renderRegion.max.x = (xMax + 63) >> 6;
497 renderRegion.min.y = yMin >> 6;
498 renderRegion.max.y = (yMax + 63) >> 6;
499 }
500 return mathClipBBox(clipRegion, renderRegion);
501}
502