edge-segments.cpp 16 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505
  1. #include "edge-segments.h"
  2. #include "arithmetics.hpp"
  3. #include "equation-solver.h"
  4. namespace msdfgen {
  5. void EdgeSegment::distanceToPseudoDistance(SignedDistance &distance, Point2 origin, double param) const {
  6. if (param < 0) {
  7. Vector2 dir = direction(0).normalize();
  8. Vector2 aq = origin-point(0);
  9. double ts = dotProduct(aq, dir);
  10. if (ts < 0) {
  11. double pseudoDistance = crossProduct(aq, dir);
  12. if (fabs(pseudoDistance) <= fabs(distance.distance)) {
  13. distance.distance = pseudoDistance;
  14. distance.dot = 0;
  15. }
  16. }
  17. } else if (param > 1) {
  18. Vector2 dir = direction(1).normalize();
  19. Vector2 bq = origin-point(1);
  20. double ts = dotProduct(bq, dir);
  21. if (ts > 0) {
  22. double pseudoDistance = crossProduct(bq, dir);
  23. if (fabs(pseudoDistance) <= fabs(distance.distance)) {
  24. distance.distance = pseudoDistance;
  25. distance.dot = 0;
  26. }
  27. }
  28. }
  29. }
  30. LinearSegment::LinearSegment(Point2 p0, Point2 p1, EdgeColor edgeColor) : EdgeSegment(edgeColor) {
  31. p[0] = p0;
  32. p[1] = p1;
  33. }
  34. QuadraticSegment::QuadraticSegment(Point2 p0, Point2 p1, Point2 p2, EdgeColor edgeColor) : EdgeSegment(edgeColor) {
  35. if (p1 == p0 || p1 == p2)
  36. p1 = 0.5*(p0+p2);
  37. p[0] = p0;
  38. p[1] = p1;
  39. p[2] = p2;
  40. }
  41. CubicSegment::CubicSegment(Point2 p0, Point2 p1, Point2 p2, Point2 p3, EdgeColor edgeColor) : EdgeSegment(edgeColor) {
  42. if ((p1 == p0 || p1 == p3) && (p2 == p0 || p2 == p3)) {
  43. p1 = mix(p0, p3, 1/3.);
  44. p2 = mix(p0, p3, 2/3.);
  45. }
  46. p[0] = p0;
  47. p[1] = p1;
  48. p[2] = p2;
  49. p[3] = p3;
  50. }
  51. LinearSegment * LinearSegment::clone() const {
  52. return new LinearSegment(p[0], p[1], color);
  53. }
  54. QuadraticSegment * QuadraticSegment::clone() const {
  55. return new QuadraticSegment(p[0], p[1], p[2], color);
  56. }
  57. CubicSegment * CubicSegment::clone() const {
  58. return new CubicSegment(p[0], p[1], p[2], p[3], color);
  59. }
  60. Point2 LinearSegment::point(double param) const {
  61. return mix(p[0], p[1], param);
  62. }
  63. Point2 QuadraticSegment::point(double param) const {
  64. return mix(mix(p[0], p[1], param), mix(p[1], p[2], param), param);
  65. }
  66. Point2 CubicSegment::point(double param) const {
  67. Vector2 p12 = mix(p[1], p[2], param);
  68. return mix(mix(mix(p[0], p[1], param), p12, param), mix(p12, mix(p[2], p[3], param), param), param);
  69. }
  70. Vector2 LinearSegment::direction(double param) const {
  71. return p[1]-p[0];
  72. }
  73. Vector2 QuadraticSegment::direction(double param) const {
  74. Vector2 tangent = mix(p[1]-p[0], p[2]-p[1], param);
  75. if (!tangent)
  76. return p[2]-p[0];
  77. return tangent;
  78. }
  79. Vector2 CubicSegment::direction(double param) const {
  80. Vector2 tangent = mix(mix(p[1]-p[0], p[2]-p[1], param), mix(p[2]-p[1], p[3]-p[2], param), param);
  81. if (!tangent) {
  82. if (param == 0) return p[2]-p[0];
  83. if (param == 1) return p[3]-p[1];
  84. }
  85. return tangent;
  86. }
  87. Vector2 LinearSegment::directionChange(double param) const {
  88. return Vector2();
  89. }
  90. Vector2 QuadraticSegment::directionChange(double param) const {
  91. return (p[2]-p[1])-(p[1]-p[0]);
  92. }
  93. Vector2 CubicSegment::directionChange(double param) const {
  94. return mix((p[2]-p[1])-(p[1]-p[0]), (p[3]-p[2])-(p[2]-p[1]), param);
  95. }
  96. double LinearSegment::length() const {
  97. return (p[1]-p[0]).length();
  98. }
  99. double QuadraticSegment::length() const {
  100. Vector2 ab = p[1]-p[0];
  101. Vector2 br = p[2]-p[1]-ab;
  102. double abab = dotProduct(ab, ab);
  103. double abbr = dotProduct(ab, br);
  104. double brbr = dotProduct(br, br);
  105. double abLen = sqrt(abab);
  106. double brLen = sqrt(brbr);
  107. double crs = crossProduct(ab, br);
  108. double h = sqrt(abab+abbr+abbr+brbr);
  109. return (
  110. brLen*((abbr+brbr)*h-abbr*abLen)+
  111. crs*crs*log((brLen*h+abbr+brbr)/(brLen*abLen+abbr))
  112. )/(brbr*brLen);
  113. }
  114. SignedDistance LinearSegment::signedDistance(Point2 origin, double &param) const {
  115. Vector2 aq = origin-p[0];
  116. Vector2 ab = p[1]-p[0];
  117. param = dotProduct(aq, ab)/dotProduct(ab, ab);
  118. Vector2 eq = p[param > .5]-origin;
  119. double endpointDistance = eq.length();
  120. if (param > 0 && param < 1) {
  121. double orthoDistance = dotProduct(ab.getOrthonormal(false), aq);
  122. if (fabs(orthoDistance) < endpointDistance)
  123. return SignedDistance(orthoDistance, 0);
  124. }
  125. return SignedDistance(nonZeroSign(crossProduct(aq, ab))*endpointDistance, fabs(dotProduct(ab.normalize(), eq.normalize())));
  126. }
  127. SignedDistance QuadraticSegment::signedDistance(Point2 origin, double &param) const {
  128. Vector2 qa = p[0]-origin;
  129. Vector2 ab = p[1]-p[0];
  130. Vector2 br = p[2]-p[1]-ab;
  131. double a = dotProduct(br, br);
  132. double b = 3*dotProduct(ab, br);
  133. double c = 2*dotProduct(ab, ab)+dotProduct(qa, br);
  134. double d = dotProduct(qa, ab);
  135. double t[3];
  136. int solutions = solveCubic(t, a, b, c, d);
  137. Vector2 epDir = direction(0);
  138. double minDistance = nonZeroSign(crossProduct(epDir, qa))*qa.length(); // distance from A
  139. param = -dotProduct(qa, epDir)/dotProduct(epDir, epDir);
  140. {
  141. epDir = direction(1);
  142. double distance = (p[2]-origin).length(); // distance from B
  143. if (distance < fabs(minDistance)) {
  144. minDistance = nonZeroSign(crossProduct(epDir, p[2]-origin))*distance;
  145. param = dotProduct(origin-p[1], epDir)/dotProduct(epDir, epDir);
  146. }
  147. }
  148. for (int i = 0; i < solutions; ++i) {
  149. if (t[i] > 0 && t[i] < 1) {
  150. Point2 qe = qa+2*t[i]*ab+t[i]*t[i]*br;
  151. double distance = qe.length();
  152. if (distance <= fabs(minDistance)) {
  153. minDistance = nonZeroSign(crossProduct(ab+t[i]*br, qe))*distance;
  154. param = t[i];
  155. }
  156. }
  157. }
  158. if (param >= 0 && param <= 1)
  159. return SignedDistance(minDistance, 0);
  160. if (param < .5)
  161. return SignedDistance(minDistance, fabs(dotProduct(direction(0).normalize(), qa.normalize())));
  162. else
  163. return SignedDistance(minDistance, fabs(dotProduct(direction(1).normalize(), (p[2]-origin).normalize())));
  164. }
  165. SignedDistance CubicSegment::signedDistance(Point2 origin, double &param) const {
  166. Vector2 qa = p[0]-origin;
  167. Vector2 ab = p[1]-p[0];
  168. Vector2 br = p[2]-p[1]-ab;
  169. Vector2 as = (p[3]-p[2])-(p[2]-p[1])-br;
  170. Vector2 epDir = direction(0);
  171. double minDistance = nonZeroSign(crossProduct(epDir, qa))*qa.length(); // distance from A
  172. param = -dotProduct(qa, epDir)/dotProduct(epDir, epDir);
  173. {
  174. epDir = direction(1);
  175. double distance = (p[3]-origin).length(); // distance from B
  176. if (distance < fabs(minDistance)) {
  177. minDistance = nonZeroSign(crossProduct(epDir, p[3]-origin))*distance;
  178. param = dotProduct(epDir-(p[3]-origin), epDir)/dotProduct(epDir, epDir);
  179. }
  180. }
  181. // Iterative minimum distance search
  182. for (int i = 0; i <= MSDFGEN_CUBIC_SEARCH_STARTS; ++i) {
  183. double t = (double) i/MSDFGEN_CUBIC_SEARCH_STARTS;
  184. Vector2 qe = qa+3*t*ab+3*t*t*br+t*t*t*as;
  185. for (int step = 0; step < MSDFGEN_CUBIC_SEARCH_STEPS; ++step) {
  186. // Improve t
  187. Vector2 d1 = 3*ab+6*t*br+3*t*t*as;
  188. Vector2 d2 = 6*br+6*t*as;
  189. t -= dotProduct(qe, d1)/(dotProduct(d1, d1)+dotProduct(qe, d2));
  190. if (t <= 0 || t >= 1)
  191. break;
  192. qe = qa+3*t*ab+3*t*t*br+t*t*t*as;
  193. double distance = qe.length();
  194. if (distance < fabs(minDistance)) {
  195. minDistance = nonZeroSign(crossProduct(d1, qe))*distance;
  196. param = t;
  197. }
  198. }
  199. }
  200. if (param >= 0 && param <= 1)
  201. return SignedDistance(minDistance, 0);
  202. if (param < .5)
  203. return SignedDistance(minDistance, fabs(dotProduct(direction(0).normalize(), qa.normalize())));
  204. else
  205. return SignedDistance(minDistance, fabs(dotProduct(direction(1).normalize(), (p[3]-origin).normalize())));
  206. }
  207. int LinearSegment::scanlineIntersections(double x[3], int dy[3], double y) const {
  208. if ((y >= p[0].y && y < p[1].y) || (y >= p[1].y && y < p[0].y)) {
  209. double param = (y-p[0].y)/(p[1].y-p[0].y);
  210. x[0] = mix(p[0].x, p[1].x, param);
  211. dy[0] = sign(p[1].y-p[0].y);
  212. return 1;
  213. }
  214. return 0;
  215. }
  216. int QuadraticSegment::scanlineIntersections(double x[3], int dy[3], double y) const {
  217. int total = 0;
  218. int nextDY = y > p[0].y ? 1 : -1;
  219. x[total] = p[0].x;
  220. if (p[0].y == y) {
  221. if (p[0].y < p[1].y || (p[0].y == p[1].y && p[0].y < p[2].y))
  222. dy[total++] = 1;
  223. else
  224. nextDY = 1;
  225. }
  226. {
  227. Vector2 ab = p[1]-p[0];
  228. Vector2 br = p[2]-p[1]-ab;
  229. double t[2];
  230. int solutions = solveQuadratic(t, br.y, 2*ab.y, p[0].y-y);
  231. // Sort solutions
  232. double tmp;
  233. if (solutions >= 2 && t[0] > t[1])
  234. tmp = t[0], t[0] = t[1], t[1] = tmp;
  235. for (int i = 0; i < solutions && total < 2; ++i) {
  236. if (t[i] >= 0 && t[i] <= 1) {
  237. x[total] = p[0].x+2*t[i]*ab.x+t[i]*t[i]*br.x;
  238. if (nextDY*(ab.y+t[i]*br.y) >= 0) {
  239. dy[total++] = nextDY;
  240. nextDY = -nextDY;
  241. }
  242. }
  243. }
  244. }
  245. if (p[2].y == y) {
  246. if (nextDY > 0 && total > 0) {
  247. --total;
  248. nextDY = -1;
  249. }
  250. if ((p[2].y < p[1].y || (p[2].y == p[1].y && p[2].y < p[0].y)) && total < 2) {
  251. x[total] = p[2].x;
  252. if (nextDY < 0) {
  253. dy[total++] = -1;
  254. nextDY = 1;
  255. }
  256. }
  257. }
  258. if (nextDY != (y >= p[2].y ? 1 : -1)) {
  259. if (total > 0)
  260. --total;
  261. else {
  262. if (fabs(p[2].y-y) < fabs(p[0].y-y))
  263. x[total] = p[2].x;
  264. dy[total++] = nextDY;
  265. }
  266. }
  267. return total;
  268. }
  269. int CubicSegment::scanlineIntersections(double x[3], int dy[3], double y) const {
  270. int total = 0;
  271. int nextDY = y > p[0].y ? 1 : -1;
  272. x[total] = p[0].x;
  273. if (p[0].y == y) {
  274. if (p[0].y < p[1].y || (p[0].y == p[1].y && (p[0].y < p[2].y || (p[0].y == p[2].y && p[0].y < p[3].y))))
  275. dy[total++] = 1;
  276. else
  277. nextDY = 1;
  278. }
  279. {
  280. Vector2 ab = p[1]-p[0];
  281. Vector2 br = p[2]-p[1]-ab;
  282. Vector2 as = (p[3]-p[2])-(p[2]-p[1])-br;
  283. double t[3];
  284. int solutions = solveCubic(t, as.y, 3*br.y, 3*ab.y, p[0].y-y);
  285. // Sort solutions
  286. double tmp;
  287. if (solutions >= 2) {
  288. if (t[0] > t[1])
  289. tmp = t[0], t[0] = t[1], t[1] = tmp;
  290. if (solutions >= 3 && t[1] > t[2]) {
  291. tmp = t[1], t[1] = t[2], t[2] = tmp;
  292. if (t[0] > t[1])
  293. tmp = t[0], t[0] = t[1], t[1] = tmp;
  294. }
  295. }
  296. for (int i = 0; i < solutions && total < 3; ++i) {
  297. if (t[i] >= 0 && t[i] <= 1) {
  298. x[total] = p[0].x+3*t[i]*ab.x+3*t[i]*t[i]*br.x+t[i]*t[i]*t[i]*as.x;
  299. if (nextDY*(ab.y+2*t[i]*br.y+t[i]*t[i]*as.y) >= 0) {
  300. dy[total++] = nextDY;
  301. nextDY = -nextDY;
  302. }
  303. }
  304. }
  305. }
  306. if (p[3].y == y) {
  307. if (nextDY > 0 && total > 0) {
  308. --total;
  309. nextDY = -1;
  310. }
  311. if ((p[3].y < p[2].y || (p[3].y == p[2].y && (p[3].y < p[1].y || (p[3].y == p[1].y && p[3].y < p[0].y)))) && total < 3) {
  312. x[total] = p[3].x;
  313. if (nextDY < 0) {
  314. dy[total++] = -1;
  315. nextDY = 1;
  316. }
  317. }
  318. }
  319. if (nextDY != (y >= p[3].y ? 1 : -1)) {
  320. if (total > 0)
  321. --total;
  322. else {
  323. if (fabs(p[3].y-y) < fabs(p[0].y-y))
  324. x[total] = p[3].x;
  325. dy[total++] = nextDY;
  326. }
  327. }
  328. return total;
  329. }
  330. static void pointBounds(Point2 p, double &l, double &b, double &r, double &t) {
  331. if (p.x < l) l = p.x;
  332. if (p.y < b) b = p.y;
  333. if (p.x > r) r = p.x;
  334. if (p.y > t) t = p.y;
  335. }
  336. void LinearSegment::bound(double &l, double &b, double &r, double &t) const {
  337. pointBounds(p[0], l, b, r, t);
  338. pointBounds(p[1], l, b, r, t);
  339. }
  340. void QuadraticSegment::bound(double &l, double &b, double &r, double &t) const {
  341. pointBounds(p[0], l, b, r, t);
  342. pointBounds(p[2], l, b, r, t);
  343. Vector2 bot = (p[1]-p[0])-(p[2]-p[1]);
  344. if (bot.x) {
  345. double param = (p[1].x-p[0].x)/bot.x;
  346. if (param > 0 && param < 1)
  347. pointBounds(point(param), l, b, r, t);
  348. }
  349. if (bot.y) {
  350. double param = (p[1].y-p[0].y)/bot.y;
  351. if (param > 0 && param < 1)
  352. pointBounds(point(param), l, b, r, t);
  353. }
  354. }
  355. void CubicSegment::bound(double &l, double &b, double &r, double &t) const {
  356. pointBounds(p[0], l, b, r, t);
  357. pointBounds(p[3], l, b, r, t);
  358. Vector2 a0 = p[1]-p[0];
  359. Vector2 a1 = 2*(p[2]-p[1]-a0);
  360. Vector2 a2 = p[3]-3*p[2]+3*p[1]-p[0];
  361. double params[2];
  362. int solutions;
  363. solutions = solveQuadratic(params, a2.x, a1.x, a0.x);
  364. for (int i = 0; i < solutions; ++i)
  365. if (params[i] > 0 && params[i] < 1)
  366. pointBounds(point(params[i]), l, b, r, t);
  367. solutions = solveQuadratic(params, a2.y, a1.y, a0.y);
  368. for (int i = 0; i < solutions; ++i)
  369. if (params[i] > 0 && params[i] < 1)
  370. pointBounds(point(params[i]), l, b, r, t);
  371. }
  372. void LinearSegment::reverse() {
  373. Point2 tmp = p[0];
  374. p[0] = p[1];
  375. p[1] = tmp;
  376. }
  377. void QuadraticSegment::reverse() {
  378. Point2 tmp = p[0];
  379. p[0] = p[2];
  380. p[2] = tmp;
  381. }
  382. void CubicSegment::reverse() {
  383. Point2 tmp = p[0];
  384. p[0] = p[3];
  385. p[3] = tmp;
  386. tmp = p[1];
  387. p[1] = p[2];
  388. p[2] = tmp;
  389. }
  390. void LinearSegment::moveStartPoint(Point2 to) {
  391. p[0] = to;
  392. }
  393. void QuadraticSegment::moveStartPoint(Point2 to) {
  394. Vector2 origSDir = p[0]-p[1];
  395. Point2 origP1 = p[1];
  396. p[1] += crossProduct(p[0]-p[1], to-p[0])/crossProduct(p[0]-p[1], p[2]-p[1])*(p[2]-p[1]);
  397. p[0] = to;
  398. if (dotProduct(origSDir, p[0]-p[1]) < 0)
  399. p[1] = origP1;
  400. }
  401. void CubicSegment::moveStartPoint(Point2 to) {
  402. p[1] += to-p[0];
  403. p[0] = to;
  404. }
  405. void LinearSegment::moveEndPoint(Point2 to) {
  406. p[1] = to;
  407. }
  408. void QuadraticSegment::moveEndPoint(Point2 to) {
  409. Vector2 origEDir = p[2]-p[1];
  410. Point2 origP1 = p[1];
  411. p[1] += crossProduct(p[2]-p[1], to-p[2])/crossProduct(p[2]-p[1], p[0]-p[1])*(p[0]-p[1]);
  412. p[2] = to;
  413. if (dotProduct(origEDir, p[2]-p[1]) < 0)
  414. p[1] = origP1;
  415. }
  416. void CubicSegment::moveEndPoint(Point2 to) {
  417. p[2] += to-p[3];
  418. p[3] = to;
  419. }
  420. void LinearSegment::splitInThirds(EdgeSegment *&part1, EdgeSegment *&part2, EdgeSegment *&part3) const {
  421. part1 = new LinearSegment(p[0], point(1/3.), color);
  422. part2 = new LinearSegment(point(1/3.), point(2/3.), color);
  423. part3 = new LinearSegment(point(2/3.), p[1], color);
  424. }
  425. void QuadraticSegment::splitInThirds(EdgeSegment *&part1, EdgeSegment *&part2, EdgeSegment *&part3) const {
  426. part1 = new QuadraticSegment(p[0], mix(p[0], p[1], 1/3.), point(1/3.), color);
  427. part2 = new QuadraticSegment(point(1/3.), mix(mix(p[0], p[1], 5/9.), mix(p[1], p[2], 4/9.), .5), point(2/3.), color);
  428. part3 = new QuadraticSegment(point(2/3.), mix(p[1], p[2], 2/3.), p[2], color);
  429. }
  430. void CubicSegment::splitInThirds(EdgeSegment *&part1, EdgeSegment *&part2, EdgeSegment *&part3) const {
  431. part1 = new CubicSegment(p[0], p[0] == p[1] ? p[0] : mix(p[0], p[1], 1/3.), mix(mix(p[0], p[1], 1/3.), mix(p[1], p[2], 1/3.), 1/3.), point(1/3.), color);
  432. part2 = new CubicSegment(point(1/3.),
  433. mix(mix(mix(p[0], p[1], 1/3.), mix(p[1], p[2], 1/3.), 1/3.), mix(mix(p[1], p[2], 1/3.), mix(p[2], p[3], 1/3.), 1/3.), 2/3.),
  434. mix(mix(mix(p[0], p[1], 2/3.), mix(p[1], p[2], 2/3.), 2/3.), mix(mix(p[1], p[2], 2/3.), mix(p[2], p[3], 2/3.), 2/3.), 1/3.),
  435. point(2/3.), color);
  436. part3 = new CubicSegment(point(2/3.), mix(mix(p[1], p[2], 2/3.), mix(p[2], p[3], 2/3.), 2/3.), p[2] == p[3] ? p[3] : mix(p[2], p[3], 2/3.), p[3], color);
  437. }
  438. EdgeSegment * QuadraticSegment::convertToCubic() const {
  439. return new CubicSegment(p[0], mix(p[0], p[1], 2/3.), mix(p[1], p[2], 1/3.), p[2], color);
  440. }
  441. void CubicSegment::deconverge(int param, double amount) {
  442. Vector2 dir = direction(param);
  443. Vector2 normal = dir.getOrthonormal();
  444. double h = dotProduct(directionChange(param)-dir, normal);
  445. switch (param) {
  446. case 0:
  447. p[1] += amount*(dir+sign(h)*sqrt(fabs(h))*normal);
  448. break;
  449. case 1:
  450. p[2] -= amount*(dir-sign(h)*sqrt(fabs(h))*normal);
  451. break;
  452. }
  453. }
  454. }