import struct
import math
import os

def merge_stls(input_files, output_path):
    try:
        from stl import mesh
        import numpy as np
        
        meshes = []
        for path, matrix in input_files:
            try:
                # 重新换回轻量级的 numpy-stl 以防内存溢出 (OOM)
                m = mesh.Mesh.from_file(path)
                
                mat = np.array(matrix, dtype=np.float64).reshape((4, 4)).T
                
                vectors = m.vectors.reshape(-1, 3)
                hom_vectors = np.hstack((vectors, np.ones((len(vectors), 1), dtype=np.float32)))
                transformed = (mat @ hom_vectors.T).T
                m.vectors = transformed[:, :3].reshape(-1, 3, 3)
                
                # 检测缩放矩阵是否引发镜像翻转 (行列式为负数)
                det = np.linalg.det(mat[:3, :3])
                if det < 0:
                    # 发生镜像反转，不仅法线会反向，三角形三个顶点的顺逆时针（Winding Order）也会错乱
                    # 强行交换每个三角形的顶点2和顶点3以纠正渲染正反面
                    m.vectors[:, [1, 2]] = m.vectors[:, [2, 1]]
                
                m.update_normals()
                meshes.append(m)
            except Exception as e:
                print(f"Error processing path {path} with stl mesh: {e}")
                
        if not meshes:
            return
            
        if len(meshes) == 1:
            meshes[0].save(output_path)
            return
            
        merged_data = np.concatenate([m.data for m in meshes])
        merged_mesh = mesh.Mesh(merged_data)
        merged_mesh.save(output_path)
        return
        
    except Exception as e:
        print(f"Mesh fast-merge failed: {e}. Falling back to struct parsing.")

    # Extreme fallback just in case no stl libraries work
    total_faces = 0
    meshes_data = []

    for path, matrix in input_files:
        with open(path, 'rb') as f:
            f.read(80)
            faces = struct.unpack('<I', f.read(4))[0]
            data = f.read(faces * 50)
            
            def apply_m(_x, _y, _z):
                w = _x * matrix[3] + _y * matrix[7] + _z * matrix[11] + matrix[15]
                nx = (_x * matrix[0] + _y * matrix[4] + _z * matrix[8] + matrix[12]) / w
                ny = (_x * matrix[1] + _y * matrix[5] + _z * matrix[9] + matrix[13]) / w
                nz = (_x * matrix[2] + _y * matrix[6] + _z * matrix[10] + matrix[14]) / w
                return nx, ny, nz
            
            new_data = bytearray(faces * 50)
            src_offset = 0
            dst_offset = 0
            for _ in range(faces):
                n_x, n_y, n_z, v1x, v1y, v1z, v2x, v2y, v2z, v3x, v3y, v3z, attr = struct.unpack_from('<12fH', data, src_offset)
                
                nv1x, nv1y, nv1z = apply_m(v1x, v1y, v1z)
                nv2x, nv2y, nv2z = apply_m(v2x, v2y, v2z)
                nv3x, nv3y, nv3z = apply_m(v3x, v3y, v3z)
                
                # Recalculate normal properly using cross product to fix flipped/sheared surfaces
                Ux = nv2x - nv1x
                Uy = nv2y - nv1y
                Uz = nv2z - nv1z
                Vx = nv3x - nv1x
                Vy = nv3y - nv1y
                Vz = nv3z - nv1z
                
                nnx = Uy * Vz - Uz * Vy
                nny = Uz * Vx - Ux * Vz
                nnz = Ux * Vy - Uy * Vx
                
                l = math.sqrt(nnx**2 + nny**2 + nnz**2)
                if l > 1e-8:
                    nnx, nny, nnz = nnx/l, nny/l, nnz/l
                else:
                    nnx, nny, nnz = 0.0, 0.0, 0.0
                
                struct.pack_into('<12fH', new_data, dst_offset, nnx, nny, nnz, nv1x, nv1y, nv1z, nv2x, nv2y, nv2z, nv3x, nv3y, nv3z, attr)
                src_offset += 50
                dst_offset += 50
                
            meshes_data.append(new_data)
            total_faces += faces
            
    with open(output_path, 'wb') as f:
        f.write(b'\0' * 80)
        f.write(struct.pack('<I', total_faces))
        for d in meshes_data:
            f.write(d)