City-Lifeline-Cockpit/.scratch/cad-bim-probe.py

"""Diagnose the massing prism volume check for one building. Run with Blender:
blender --background --factory-startup --python-exit-code 1 --python <this> -- <plan.json> <building-id>
"""
import json
import sys
from pathlib import Path

import bpy
import bmesh
from mathutils import Vector

argv = sys.argv[sys.argv.index('--') + 1:]
plan = json.loads(Path(argv[0]).read_text(encoding='utf-8'))
target = argv[1] if len(argv) > 1 else 'cad-210'
building = next(b for b in plan['buildings'] if b['id'] == target)


def ring_area(points):
    total = 0.0
    for index in range(len(points)):
        x1, y1 = points[index]
        x2, y2 = points[(index + 1) % len(points)]
        total += x1 * y2 - x2 * y1
    return total / 2.0


def as_ccw(points):
    return list(reversed(points)) if ring_area(points) < 0 else list(points)


def volume_from(obj, anchor):
    total = 0.0
    for polygon in obj.data.polygons:
        pts = [obj.data.vertices[i].co - anchor for i in polygon.vertices]
        for k in range(1, len(pts) - 1):
            total += pts[0].cross(pts[k]).dot(pts[k + 1])
    return total


outline = [(float(x), float(y)) for x, y in building['footprintENU']]
print('raw ring area   ', round(ring_area(outline), 3))
outline = as_ccw(outline)
print('ccw ring area   ', round(ring_area(outline), 3))
print('positionENU     ', building['positionENU'])
print('foundationBottom', building['foundationBottom'], 'heightM', building['heightM'])

bottom = building['foundationBottom']
top = building['positionENU'][2] + building['heightM']
ax, ay, az = building['positionENU']
n = len(outline)
vertices = [(x - ax, y - ay, bottom - az) for x, y in outline] + [(x - ax, y - ay, top - az) for x, y in outline]

bpy.ops.wm.read_factory_settings(use_empty=True)
faces = [(i, (i + 1) % n, n + (i + 1) % n, n + i) for i in range(n)]
faces += [tuple(range(n)), tuple(reversed(range(n, 2 * n)))]
faces += [(n + i, n + (i + 1) % n, (i + 1) % n, i) for i in range(n)]
print('faces           ', faces)

mesh = bpy.data.meshes.new('probe')
mesh.from_pydata(vertices, [], faces)
mesh.update()
obj = bpy.data.objects.new('probe', mesh)
bpy.context.scene.collection.objects.link(obj)

anchor = Vector((0.0, 0.0, 0.0))
for vertex in obj.data.vertices:
    anchor += vertex.co
anchor /= len(obj.data.vertices)
print('local anchor    ', tuple(round(v, 3) for v in anchor))
print('local verts     ', [tuple(round(c, 2) for c in v.co) for v in obj.data.vertices])
print('volume pre-recalc ', round(volume_from(obj, anchor), 3))

edit = bmesh.new()
edit.from_mesh(mesh)
bmesh.ops.recalc_face_normals(edit, faces=edit.faces)
edit.to_mesh(mesh)
edit.free()
mesh.update()
print('volume post-recalc', round(volume_from(obj, anchor), 3))

print('polygons        ', len(mesh.polygons), 'expected', n + 2)
bm = bmesh.new()
bm.from_mesh(mesh)
print('non-manifold edges', [e.index for e in bm.edges if len(e.link_faces) != 2])
bm.free()