seat-cad-bim.mjs 6.6 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173
  1. // Seat every CAD BIM building on the original DEM.
  2. // Flat ground is seated directly; a gently sloped footprint gets a bounded stone plinth; steeper
  3. // footprints are refused rather than hidden behind a fabricated platform.
  4. import { readFileSync, writeFileSync } from 'node:fs'
  5. import { Cartesian3, Matrix4 } from 'cesium'
  6. const dir = 'tools/city-massing/city/cad-bim/'
  7. const plan = JSON.parse(readFileSync(dir + 'plan.json'))
  8. const sampled = JSON.parse(readFileSync(dir + 'terrain-samples.json'))
  9. if (plan.samples.length !== sampled.heights.length) throw new Error('DEM sample count does not match the plan')
  10. const FLAT_RELIEF = 2.5
  11. const MAX_PLINTH = 3.2
  12. const PLINTH_HEIGHT_RATIO = 0.4
  13. const FLOOR_HEIGHT_M = 3.2
  14. const MAX_HEIGHT_M = 30
  15. const inverse = Matrix4.inverseTransformation(Matrix4.fromArray(plan.transform), new Matrix4())
  16. const local = (lon, lat, h) => Matrix4.multiplyByPoint(inverse, Cartesian3.fromDegrees(lon, lat, h), new Cartesian3())
  17. // Survey outlines occasionally fold back on themselves (a "bow tie"), repeat a vertex, or enclose a
  18. // negligible area. Those cannot be extruded into a solid, so they are rejected here — with a reason —
  19. // instead of reaching the exporter and aborting a full-city run.
  20. const MIN_MASSING_AREA = 4
  21. const DUPLICATE_VERTEX_M = 0.05
  22. function cleanRing(points) {
  23. const out = []
  24. for (const point of points) {
  25. const previous = out[out.length - 1]
  26. if (previous && Math.hypot(point[0] - previous[0], point[1] - previous[1]) < DUPLICATE_VERTEX_M) continue
  27. out.push(point)
  28. }
  29. if (out.length > 2) {
  30. const first = out[0], last = out[out.length - 1]
  31. if (Math.hypot(first[0] - last[0], first[1] - last[1]) < DUPLICATE_VERTEX_M) out.pop()
  32. }
  33. return out
  34. }
  35. function signedArea(points) {
  36. let total = 0
  37. for (let i = 0; i < points.length; i++) {
  38. const [x1, y1] = points[i]
  39. const [x2, y2] = points[(i + 1) % points.length]
  40. total += x1 * y2 - x2 * y1
  41. }
  42. return total / 2
  43. }
  44. function segmentsCross(a, b) {
  45. const cross = (o, p, q) => (p[0] - o[0]) * (q[1] - o[1]) - (p[1] - o[1]) * (q[0] - o[0])
  46. const d1 = cross(a[0], a[1], b[0]), d2 = cross(a[0], a[1], b[1])
  47. const d3 = cross(b[0], b[1], a[0]), d4 = cross(b[0], b[1], a[1])
  48. return ((d1 > 0) !== (d2 > 0)) && ((d3 > 0) !== (d4 > 0))
  49. }
  50. /** True when no two non-adjacent edges of the ring intersect. */
  51. function isSimpleRing(points) {
  52. const n = points.length
  53. for (let i = 0; i < n; i++) {
  54. const a = [points[i], points[(i + 1) % n]]
  55. for (let j = i + 1; j < n; j++) {
  56. if (j === i || (j + 1) % n === i || j === (i + 1) % n) continue
  57. const b = [points[j], points[(j + 1) % n]]
  58. if (segmentsCross(a, b)) return false
  59. }
  60. }
  61. return true
  62. }
  63. const buildings = []
  64. const excluded = []
  65. const construction = { flat: 0, plinth: 0 }
  66. let missing = 0
  67. let invalidGeometry = 0
  68. for (const building of plan.buildings) {
  69. const heights = building.sampleIndices.map((i) => sampled.heights[i])
  70. if (heights.some((h) => !Number.isFinite(h) || h <= 1)) {
  71. missing++
  72. excluded.push({ id: building.id, reason: 'missing or invalid DEM sample' })
  73. continue
  74. }
  75. let ring = building.footprintENU
  76. if (building.kind === 'massing') {
  77. ring = cleanRing(ring)
  78. if (ring.length < 3) {
  79. invalidGeometry++
  80. excluded.push({ id: building.id, tier: building.tier, reason: 'outline collapses to fewer than 3 distinct points' })
  81. continue
  82. }
  83. if (!isSimpleRing(ring)) {
  84. invalidGeometry++
  85. excluded.push({ id: building.id, tier: building.tier, reason: 'self-intersecting survey outline cannot be extruded' })
  86. continue
  87. }
  88. if (Math.abs(signedArea(ring)) < MIN_MASSING_AREA) {
  89. invalidGeometry++
  90. excluded.push({ id: building.id, tier: building.tier, reason: 'outline encloses negligible area' })
  91. continue
  92. }
  93. }
  94. const min = Math.min(...heights)
  95. const max = Math.max(...heights)
  96. const relief = max - min
  97. const heightM = Math.min(MAX_HEIGHT_M, Math.max(FLOOR_HEIGHT_M, building.floors * FLOOR_HEIGHT_M))
  98. if (relief > FLAT_RELIEF) {
  99. if (relief > MAX_PLINTH || relief > heightM * PLINTH_HEIGHT_RATIO) {
  100. excluded.push({
  101. id: building.id, tier: building.tier,
  102. reason: 'terrain relief exceeds the allowed plinth; needs stepped construction',
  103. relief: Math.round(relief * 100) / 100,
  104. })
  105. continue
  106. }
  107. construction.plinth++
  108. } else {
  109. construction.flat++
  110. }
  111. // Floor sits just below the highest measured ground so the building never floats.
  112. const floorZ = max - 0.25
  113. const centroid = building.sampleIndices[building.sampleIndices.length - 1]
  114. const anchor = plan.samples[centroid]
  115. buildings.push({
  116. id: building.id,
  117. kind: building.kind,
  118. tier: building.tier,
  119. maskId: building.maskId,
  120. footprintENU: ring,
  121. width: building.width,
  122. depth: building.depth,
  123. heading: building.heading,
  124. archetype: building.archetype,
  125. preferredAssets: building.preferredAssets,
  126. floors: building.floors,
  127. heightM: Math.round(heightM * 100) / 100,
  128. heightSource: building.heightSource,
  129. positionENU: [
  130. Math.round(local(anchor.lon, anchor.lat, floorZ).x * 100) / 100,
  131. Math.round(local(anchor.lon, anchor.lat, floorZ).y * 100) / 100,
  132. Math.round(floorZ * 100) / 100,
  133. ],
  134. foundationBottom: Math.round((min - 3) * 100) / 100,
  135. terrainMin: Math.round(min * 100) / 100,
  136. terrainMax: Math.round(max * 100) / 100,
  137. terrainRelief: Math.round(relief * 100) / 100,
  138. construction: relief > FLAT_RELIEF ? 'plinth' : 'flat',
  139. })
  140. }
  141. const out = {
  142. ...plan,
  143. buildings,
  144. excluded: [...(plan.excluded || []), ...excluded],
  145. construction,
  146. notice: plan.disclaimer + ' Seated on the original quantized-mesh DEM; steep footprints stay excluded.',
  147. }
  148. writeFileSync(dir + 'seated-plan.json', JSON.stringify(out, null, 2))
  149. console.log(JSON.stringify({
  150. planned: plan.buildings.length,
  151. seated: buildings.length,
  152. byKind: {
  153. asset: buildings.filter((b) => b.kind === 'asset').length,
  154. massing: buildings.filter((b) => b.kind === 'massing').length,
  155. },
  156. byConstruction: construction,
  157. missingSamples: missing,
  158. invalidGeometry,
  159. excludedRelief: excluded.filter((e) => e.relief !== undefined).length,
  160. urbanSeated: buildings.filter((b) => b.tier === 'urban').length,
  161. ruralSeated: buildings.filter((b) => b.tier === 'rural').length,
  162. }))