centroid.js 3.7 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140
  1. import {asin, atan2, cos, degrees, epsilon, epsilon2, radians, sin, sqrt} from "./math";
  2. import noop from "./noop";
  3. import stream from "./stream";
  4. var W0, W1,
  5. X0, Y0, Z0,
  6. X1, Y1, Z1,
  7. X2, Y2, Z2,
  8. lambda00, phi00, // first point
  9. x0, y0, z0; // previous point
  10. var centroidStream = {
  11. sphere: noop,
  12. point: centroidPoint,
  13. lineStart: centroidLineStart,
  14. lineEnd: centroidLineEnd,
  15. polygonStart: function() {
  16. centroidStream.lineStart = centroidRingStart;
  17. centroidStream.lineEnd = centroidRingEnd;
  18. },
  19. polygonEnd: function() {
  20. centroidStream.lineStart = centroidLineStart;
  21. centroidStream.lineEnd = centroidLineEnd;
  22. }
  23. };
  24. // Arithmetic mean of Cartesian vectors.
  25. function centroidPoint(lambda, phi) {
  26. lambda *= radians, phi *= radians;
  27. var cosPhi = cos(phi);
  28. centroidPointCartesian(cosPhi * cos(lambda), cosPhi * sin(lambda), sin(phi));
  29. }
  30. function centroidPointCartesian(x, y, z) {
  31. ++W0;
  32. X0 += (x - X0) / W0;
  33. Y0 += (y - Y0) / W0;
  34. Z0 += (z - Z0) / W0;
  35. }
  36. function centroidLineStart() {
  37. centroidStream.point = centroidLinePointFirst;
  38. }
  39. function centroidLinePointFirst(lambda, phi) {
  40. lambda *= radians, phi *= radians;
  41. var cosPhi = cos(phi);
  42. x0 = cosPhi * cos(lambda);
  43. y0 = cosPhi * sin(lambda);
  44. z0 = sin(phi);
  45. centroidStream.point = centroidLinePoint;
  46. centroidPointCartesian(x0, y0, z0);
  47. }
  48. function centroidLinePoint(lambda, phi) {
  49. lambda *= radians, phi *= radians;
  50. var cosPhi = cos(phi),
  51. x = cosPhi * cos(lambda),
  52. y = cosPhi * sin(lambda),
  53. z = sin(phi),
  54. w = atan2(sqrt((w = y0 * z - z0 * y) * w + (w = z0 * x - x0 * z) * w + (w = x0 * y - y0 * x) * w), x0 * x + y0 * y + z0 * z);
  55. W1 += w;
  56. X1 += w * (x0 + (x0 = x));
  57. Y1 += w * (y0 + (y0 = y));
  58. Z1 += w * (z0 + (z0 = z));
  59. centroidPointCartesian(x0, y0, z0);
  60. }
  61. function centroidLineEnd() {
  62. centroidStream.point = centroidPoint;
  63. }
  64. // See J. E. Brock, The Inertia Tensor for a Spherical Triangle,
  65. // J. Applied Mechanics 42, 239 (1975).
  66. function centroidRingStart() {
  67. centroidStream.point = centroidRingPointFirst;
  68. }
  69. function centroidRingEnd() {
  70. centroidRingPoint(lambda00, phi00);
  71. centroidStream.point = centroidPoint;
  72. }
  73. function centroidRingPointFirst(lambda, phi) {
  74. lambda00 = lambda, phi00 = phi;
  75. lambda *= radians, phi *= radians;
  76. centroidStream.point = centroidRingPoint;
  77. var cosPhi = cos(phi);
  78. x0 = cosPhi * cos(lambda);
  79. y0 = cosPhi * sin(lambda);
  80. z0 = sin(phi);
  81. centroidPointCartesian(x0, y0, z0);
  82. }
  83. function centroidRingPoint(lambda, phi) {
  84. lambda *= radians, phi *= radians;
  85. var cosPhi = cos(phi),
  86. x = cosPhi * cos(lambda),
  87. y = cosPhi * sin(lambda),
  88. z = sin(phi),
  89. cx = y0 * z - z0 * y,
  90. cy = z0 * x - x0 * z,
  91. cz = x0 * y - y0 * x,
  92. m = sqrt(cx * cx + cy * cy + cz * cz),
  93. w = asin(m), // line weight = angle
  94. v = m && -w / m; // area weight multiplier
  95. X2 += v * cx;
  96. Y2 += v * cy;
  97. Z2 += v * cz;
  98. W1 += w;
  99. X1 += w * (x0 + (x0 = x));
  100. Y1 += w * (y0 + (y0 = y));
  101. Z1 += w * (z0 + (z0 = z));
  102. centroidPointCartesian(x0, y0, z0);
  103. }
  104. export default function(object) {
  105. W0 = W1 =
  106. X0 = Y0 = Z0 =
  107. X1 = Y1 = Z1 =
  108. X2 = Y2 = Z2 = 0;
  109. stream(object, centroidStream);
  110. var x = X2,
  111. y = Y2,
  112. z = Z2,
  113. m = x * x + y * y + z * z;
  114. // If the area-weighted ccentroid is undefined, fall back to length-weighted ccentroid.
  115. if (m < epsilon2) {
  116. x = X1, y = Y1, z = Z1;
  117. // If the feature has zero length, fall back to arithmetic mean of point vectors.
  118. if (W1 < epsilon) x = X0, y = Y0, z = Z0;
  119. m = x * x + y * y + z * z;
  120. // If the feature still has an undefined ccentroid, then return.
  121. if (m < epsilon2) return [NaN, NaN];
  122. }
  123. return [atan2(y, x) * degrees, asin(z / sqrt(m)) * degrees];
  124. }