import struct, numpy as np

path = '/app/projects/esp32-box/original/Window 20x60mm.stl'

# Load binary STL
with open(path, 'rb') as f:
    f.read(80)
    n = struct.unpack('<I', f.read(4))[0]
    verts = []
    normals = []
    for i in range(n):
        nx, ny, nz = struct.unpack('<fff', f.read(12))
        normals.append([nx, ny, nz])
        tri = []
        for j in range(3):
            x, y, z = struct.unpack('<fff', f.read(12))
            tri.append([x, y, z])
        f.read(2)
        verts.append(tri)

verts = np.array(verts)
print('Total triangles:', n)

# Bounding box
all_pts = verts.reshape(-1, 3)
print('Bounding box:')
print('  X: %.2f to %.2f (%.2f)' % (all_pts[:,0].min(), all_pts[:,0].max(), all_pts[:,0].max()-all_pts[:,0].min()))
print('  Y: %.2f to %.2f (%.2f)' % (all_pts[:,1].min(), all_pts[:,1].max(), all_pts[:,1].max()-all_pts[:,1].min()))
print('  Z: %.2f to %.2f (%.2f)' % (all_pts[:,2].min(), all_pts[:,2].max(), all_pts[:,2].max()-all_pts[:,2].min()))

# Cross-sections at multiple X positions
print()
print('=== Cross-sections at various X positions ===')
x_min, x_max = all_pts[:,0].min(), all_pts[:,0].max()
x_positions = np.linspace(x_min+0.1, x_max-0.1, 15)

for x_pos in x_positions:
    # Find all triangles that cross this X plane
    yz_points = set()
    for tri in verts:
        xs = tri[:,0]
        if (xs >= x_pos).any() and (xs <= x_pos).any():
            # Find intersection edges
            for i in range(3):
                j = (i+1) % 3
                if (xs[i] - x_pos) * (xs[j] - x_pos) < 0:
                    t = (x_pos - xs[i]) / (xs[j] - xs[i])
                    y = tri[i,1] + t * (tri[j,1] - tri[i,1])
                    z = tri[i,2] + t * (tri[j,2] - tri[i,2])
                    yz_points.add((round(y,2), round(z,2)))
    if yz_points:
        yz = np.array(sorted(yz_points))
        print('X=%.2f: Y=[%.2f..%.2f] Z=[%.2f..%.2f] (%d pts)' % (
            x_pos, yz[:,0].min(), yz[:,0].max(), yz[:,1].min(), yz[:,1].max(), len(yz)))

# Cross-sections at multiple Y positions
print()
print('=== Cross-sections at various Y positions ===')
y_min, y_max = all_pts[:,1].min(), all_pts[:,1].max()
y_positions = np.linspace(y_min+0.1, y_max-0.1, 15)

for y_pos in y_positions:
    xz_points = set()
    for tri in verts:
        ys = tri[:,1]
        if (ys >= y_pos).any() and (ys <= y_pos).any():
            for i in range(3):
                j = (i+1) % 3
                if (ys[i] - y_pos) * (ys[j] - y_pos) < 0:
                    t = (y_pos - ys[i]) / (ys[j] - ys[i])
                    x = tri[i,0] + t * (tri[j,0] - tri[i,0])
                    z = tri[i,2] + t * (tri[j,2] - tri[i,2])
                    xz_points.add((round(x,2), round(z,2)))
    if xz_points:
        xz = np.array(sorted(xz_points))
        print('Y=%.2f: X=[%.2f..%.2f] Z=[%.2f..%.2f] (%d pts)' % (
            y_pos, xz[:,0].min(), xz[:,0].max(), xz[:,1].min(), xz[:,1].max(), len(xz)))

# Cross-sections at Z positions (top-down view)
print()
print('=== Cross-sections at various Z positions ===')
z_min, z_max = all_pts[:,2].min(), all_pts[:,2].max()
z_positions = np.linspace(z_min+0.05, z_max-0.05, 10)

for z_pos in z_positions:
    xy_points = set()
    for tri in verts:
        zs = tri[:,2]
        if (zs >= z_pos).any() and (zs <= z_pos).any():
            for i in range(3):
                j = (i+1) % 3
                if (zs[i] - z_pos) * (zs[j] - z_pos) < 0:
                    t = (z_pos - zs[i]) / (zs[j] - zs[i])
                    x = tri[i,0] + t * (tri[j,0] - tri[i,0])
                    y = tri[i,1] + t * (tri[j,1] - tri[i,1])
                    xy_points.add((round(x,2), round(y,2)))
    if xy_points:
        xy = np.array(sorted(xy_points))
        print('Z=%.2f: X=[%.2f..%.2f] Y=[%.2f..%.2f] (%d pts)' % (
            z_pos, xy[:,0].min(), xy[:,0].max(), xy[:,1].min(), xy[:,1].max(), len(xy)))

# Identify unique Z levels (peaks, flat surfaces)
print()
print('=== Z level analysis ===')
z_vals = np.round(all_pts[:,2], 2)
unique_z = np.unique(z_vals)
print('Unique Z values:', unique_z)

# Count vertices at each Z level
for z in unique_z:
    count = np.sum(z_vals == z)
    print('  Z=%.2f: %d vertices' % (z, count))