TrianglePredicate.js 5.6 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106
  1. import WKTWriter from '../../io/WKTWriter'
  2. import CoordinateArraySequence from '../../geom/impl/CoordinateArraySequence'
  3. import DD from '../../math/DD'
  4. import System from '../../../../../java/lang/System'
  5. import Triangle from '../../geom/Triangle'
  6. export default class TrianglePredicate {
  7. static triArea(a, b, c) {
  8. return (b.x - a.x) * (c.y - a.y) - (b.y - a.y) * (c.x - a.x)
  9. }
  10. static isInCircleDDNormalized(a, b, c, p) {
  11. const adx = DD.valueOf(a.x).selfSubtract(p.x)
  12. const ady = DD.valueOf(a.y).selfSubtract(p.y)
  13. const bdx = DD.valueOf(b.x).selfSubtract(p.x)
  14. const bdy = DD.valueOf(b.y).selfSubtract(p.y)
  15. const cdx = DD.valueOf(c.x).selfSubtract(p.x)
  16. const cdy = DD.valueOf(c.y).selfSubtract(p.y)
  17. const abdet = adx.multiply(bdy).selfSubtract(bdx.multiply(ady))
  18. const bcdet = bdx.multiply(cdy).selfSubtract(cdx.multiply(bdy))
  19. const cadet = cdx.multiply(ady).selfSubtract(adx.multiply(cdy))
  20. const alift = adx.multiply(adx).selfAdd(ady.multiply(ady))
  21. const blift = bdx.multiply(bdx).selfAdd(bdy.multiply(bdy))
  22. const clift = cdx.multiply(cdx).selfAdd(cdy.multiply(cdy))
  23. const sum = alift.selfMultiply(bcdet).selfAdd(blift.selfMultiply(cadet)).selfAdd(clift.selfMultiply(abdet))
  24. const isInCircle = sum.doubleValue() > 0
  25. return isInCircle
  26. }
  27. static checkRobustInCircle(a, b, c, p) {
  28. const nonRobustInCircle = TrianglePredicate.isInCircleNonRobust(a, b, c, p)
  29. const isInCircleDD = TrianglePredicate.isInCircleDDSlow(a, b, c, p)
  30. const isInCircleCC = TrianglePredicate.isInCircleCC(a, b, c, p)
  31. const circumCentre = Triangle.circumcentre(a, b, c)
  32. System.out.println('p radius diff a = ' + Math.abs(p.distance(circumCentre) - a.distance(circumCentre)) / a.distance(circumCentre))
  33. if (nonRobustInCircle !== isInCircleDD || nonRobustInCircle !== isInCircleCC) {
  34. System.out.println('inCircle robustness failure (double result = ' + nonRobustInCircle + ', DD result = ' + isInCircleDD + ', CC result = ' + isInCircleCC + ')')
  35. System.out.println(WKTWriter.toLineString(new CoordinateArraySequence([a, b, c, p])))
  36. System.out.println('Circumcentre = ' + WKTWriter.toPoint(circumCentre) + ' radius = ' + a.distance(circumCentre))
  37. System.out.println('p radius diff a = ' + Math.abs(p.distance(circumCentre) / a.distance(circumCentre) - 1))
  38. System.out.println('p radius diff b = ' + Math.abs(p.distance(circumCentre) / b.distance(circumCentre) - 1))
  39. System.out.println('p radius diff c = ' + Math.abs(p.distance(circumCentre) / c.distance(circumCentre) - 1))
  40. System.out.println()
  41. }
  42. }
  43. static isInCircleDDFast(a, b, c, p) {
  44. const aTerm = DD.sqr(a.x).selfAdd(DD.sqr(a.y)).selfMultiply(TrianglePredicate.triAreaDDFast(b, c, p))
  45. const bTerm = DD.sqr(b.x).selfAdd(DD.sqr(b.y)).selfMultiply(TrianglePredicate.triAreaDDFast(a, c, p))
  46. const cTerm = DD.sqr(c.x).selfAdd(DD.sqr(c.y)).selfMultiply(TrianglePredicate.triAreaDDFast(a, b, p))
  47. const pTerm = DD.sqr(p.x).selfAdd(DD.sqr(p.y)).selfMultiply(TrianglePredicate.triAreaDDFast(a, b, c))
  48. const sum = aTerm.selfSubtract(bTerm).selfAdd(cTerm).selfSubtract(pTerm)
  49. const isInCircle = sum.doubleValue() > 0
  50. return isInCircle
  51. }
  52. static isInCircleCC(a, b, c, p) {
  53. const cc = Triangle.circumcentre(a, b, c)
  54. const ccRadius = a.distance(cc)
  55. const pRadiusDiff = p.distance(cc) - ccRadius
  56. return pRadiusDiff <= 0
  57. }
  58. static isInCircleNormalized(a, b, c, p) {
  59. const adx = a.x - p.x
  60. const ady = a.y - p.y
  61. const bdx = b.x - p.x
  62. const bdy = b.y - p.y
  63. const cdx = c.x - p.x
  64. const cdy = c.y - p.y
  65. const abdet = adx * bdy - bdx * ady
  66. const bcdet = bdx * cdy - cdx * bdy
  67. const cadet = cdx * ady - adx * cdy
  68. const alift = adx * adx + ady * ady
  69. const blift = bdx * bdx + bdy * bdy
  70. const clift = cdx * cdx + cdy * cdy
  71. const disc = alift * bcdet + blift * cadet + clift * abdet
  72. return disc > 0
  73. }
  74. static isInCircleDDSlow(a, b, c, p) {
  75. const px = DD.valueOf(p.x)
  76. const py = DD.valueOf(p.y)
  77. const ax = DD.valueOf(a.x)
  78. const ay = DD.valueOf(a.y)
  79. const bx = DD.valueOf(b.x)
  80. const by = DD.valueOf(b.y)
  81. const cx = DD.valueOf(c.x)
  82. const cy = DD.valueOf(c.y)
  83. const aTerm = ax.multiply(ax).add(ay.multiply(ay)).multiply(TrianglePredicate.triAreaDDSlow(bx, by, cx, cy, px, py))
  84. const bTerm = bx.multiply(bx).add(by.multiply(by)).multiply(TrianglePredicate.triAreaDDSlow(ax, ay, cx, cy, px, py))
  85. const cTerm = cx.multiply(cx).add(cy.multiply(cy)).multiply(TrianglePredicate.triAreaDDSlow(ax, ay, bx, by, px, py))
  86. const pTerm = px.multiply(px).add(py.multiply(py)).multiply(TrianglePredicate.triAreaDDSlow(ax, ay, bx, by, cx, cy))
  87. const sum = aTerm.subtract(bTerm).add(cTerm).subtract(pTerm)
  88. const isInCircle = sum.doubleValue() > 0
  89. return isInCircle
  90. }
  91. static isInCircleNonRobust(a, b, c, p) {
  92. const isInCircle = (a.x * a.x + a.y * a.y) * TrianglePredicate.triArea(b, c, p) - (b.x * b.x + b.y * b.y) * TrianglePredicate.triArea(a, c, p) + (c.x * c.x + c.y * c.y) * TrianglePredicate.triArea(a, b, p) - (p.x * p.x + p.y * p.y) * TrianglePredicate.triArea(a, b, c) > 0
  93. return isInCircle
  94. }
  95. static isInCircleRobust(a, b, c, p) {
  96. return TrianglePredicate.isInCircleNormalized(a, b, c, p)
  97. }
  98. static triAreaDDSlow(ax, ay, bx, by, cx, cy) {
  99. return bx.subtract(ax).multiply(cy.subtract(ay)).subtract(by.subtract(ay).multiply(cx.subtract(ax)))
  100. }
  101. static triAreaDDFast(a, b, c) {
  102. const t1 = DD.valueOf(b.x).selfSubtract(a.x).selfMultiply(DD.valueOf(c.y).selfSubtract(a.y))
  103. const t2 = DD.valueOf(b.y).selfSubtract(a.y).selfMultiply(DD.valueOf(c.x).selfSubtract(a.x))
  104. return t1.selfSubtract(t2)
  105. }
  106. }