src/lib/sw_engine/tvgSwMath.cpp

   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
  45 constexpr 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
  48 constexpr static auto ATAN_MAX = 23;
  49 constexpr 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
  54 static inline SwCoord SATURATE(const SwCoord x)
  55 {
  56     return (x >> (sizeof(SwCoord) * 8 - 1));
  57 }
  58
  59
  60 static inline SwFixed PAD_ROUND(const SwFixed x, int32_t n)
  61 {
  62     return (((x) + ((n)/2)) & ~((n)-1));
  63 }
  64
  65
  66 static 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
  77 static 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
 102 static 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
 158 static 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
 205 SwFixed mathMean(SwFixed angle1, SwFixed angle2)
 206 {
 207     return angle1 + mathDiff(angle1, angle2) / 2;
 208 }
 209
 210
 211 bool 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
 262 int64_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
 280 int64_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
 298 int64_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
 321 void 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
 345 SwFixed 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
 353 SwFixed 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
 365 SwFixed mathSin(SwFixed angle)
 366 {
 367     return mathCos(SW_ANGLE_PI2 - angle);
 368 }
 369
 370
 371 SwFixed mathCos(SwFixed angle)
 372 {
 373     SwPoint v = {CORDIC_FACTOR >> 8, 0};
 374     _rotate(v, angle);
 375     return (v.x + 0x80L) >> 8;
 376 }
 377
 378
 379 SwFixed 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
 397 void 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
 423 SwFixed 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
 435 SwPoint 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
 446 bool 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
 464 bool mathUpdateOutlineBBox(const SwOutline* outline, const SwBBox& clipRegion, SwBBox& renderRegion, bool fastTrack)
 465 {
 466     if (!outline) return false;
 467
 468     auto pt = outline->pts;
 469
 470     if (outline->ptsCnt == 0 || outline->cntrsCnt <= 0) {
 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     ++pt;
 481
 482     for (uint32_t i = 1; i < outline->ptsCnt; ++i, ++pt) {
 483         if (xMin > pt->x) xMin = pt->x;
 484         if (xMax < pt->x) xMax = pt->x;
 485         if (yMin > pt->y) yMin = pt->y;
 486         if (yMax < pt->y) yMax = pt->y;
 487     }
 488     //Since no antialiasing is applied in the Fast Track case,
 489     //the rasterization region has to be rearranged.
 490     //https://github.com/Samsung/thorvg/issues/916
 491     if (fastTrack) {
 492         renderRegion.min.x = static_cast<SwCoord>(round(xMin / 64.0f));
 493         renderRegion.max.x = static_cast<SwCoord>(round(xMax / 64.0f));
 494         renderRegion.min.y = static_cast<SwCoord>(round(yMin / 64.0f));
 495         renderRegion.max.y = static_cast<SwCoord>(round(yMax / 64.0f));
 496     } else {
 497         renderRegion.min.x = xMin >> 6;
 498         renderRegion.max.x = (xMax + 63) >> 6;
 499         renderRegion.min.y = yMin >> 6;
 500         renderRegion.max.y = (yMax + 63) >> 6;
 501     }
 502     return mathClipBBox(clipRegion, renderRegion);
 503 }