modeling/src/maths/plane/fromNoisyPoints.js

  1. const vec3 = require('../vec3')
  2. const fromNormalAndPoint = require('./fromNormalAndPoint')
  3. /**
  4. * Create a best-fit plane from the given noisy vertices.
  5. *
  6. * NOTE: There are two possible orientations for every plane.
  7. * This function always produces positive orientations.
  8. *
  9. * See http://www.ilikebigbits.com for the original discussion
  10. *
  11. * @param {Plane} out - receiving plane
  12. * @param {Array} vertices - list of vertices in any order or position
  13. * @returns {Plane} out
  14. * @alias module:modeling/maths/plane.fromNoisyPoints
  15. */
  16. const fromNoisyPoints = (out, ...vertices) => {
  17. out[0] = 0.0
  18. out[1] = 0.0
  19. out[2] = 0.0
  20. out[3] = 0.0
  21. // calculate the centroid of the vertices
  22. // NOTE: out is the centriod
  23. const n = vertices.length
  24. vertices.forEach((v) => {
  25. vec3.add(out, out, v)
  26. })
  27. vec3.scale(out, out, 1.0 / n)
  28. // Calculate full 3x3 covariance matrix, excluding symmetries
  29. let xx = 0.0
  30. let xy = 0.0
  31. let xz = 0.0
  32. let yy = 0.0
  33. let yz = 0.0
  34. let zz = 0.0
  35. const vn = vec3.create()
  36. vertices.forEach((v) => {
  37. // NOTE: out is the centriod
  38. vec3.subtract(vn, v, out)
  39. xx += vn[0] * vn[0]
  40. xy += vn[0] * vn[1]
  41. xz += vn[0] * vn[2]
  42. yy += vn[1] * vn[1]
  43. yz += vn[1] * vn[2]
  44. zz += vn[2] * vn[2]
  45. })
  46. xx /= n
  47. xy /= n
  48. xz /= n
  49. yy /= n
  50. yz /= n
  51. zz /= n
  52. // Calculate the smallest Eigenvector of the covariance matrix
  53. // which becomes the plane normal
  54. vn[0] = 0.0
  55. vn[1] = 0.0
  56. vn[2] = 0.0
  57. // weighted directional vector
  58. const wdv = vec3.create()
  59. // X axis
  60. let det = yy * zz - yz * yz
  61. wdv[0] = det
  62. wdv[1] = xz * yz - xy * zz
  63. wdv[2] = xy * yz - xz * yy
  64. let weight = det * det
  65. vec3.add(vn, vn, vec3.scale(wdv, wdv, weight))
  66. // Y axis
  67. det = xx * zz - xz * xz
  68. wdv[0] = xz * yz - xy * zz
  69. wdv[1] = det
  70. wdv[2] = xy * xz - yz * xx
  71. weight = det * det
  72. if (vec3.dot(vn, wdv) < 0.0) {
  73. weight = -weight
  74. }
  75. vec3.add(vn, vn, vec3.scale(wdv, wdv, weight))
  76. // Z axis
  77. det = xx * yy - xy * xy
  78. wdv[0] = xy * yz - xz * yy
  79. wdv[1] = xy * xz - yz * xx
  80. wdv[2] = det
  81. weight = det * det
  82. if (vec3.dot(vn, wdv) < 0.0) {
  83. weight = -weight
  84. }
  85. vec3.add(vn, vn, vec3.scale(wdv, wdv, weight))
  86. // create the plane from normal and centriod
  87. // NOTE: out is the centriod
  88. return fromNormalAndPoint(out, vn, out)
  89. }
  90. module.exports = fromNoisyPoints