cdshelat/MorphAI

1import io2import struct3import numpy as np4from scipy.interpolate import RegularGridInterpolator5from skimage import measure6from materials import INFILL_PATTERNS7 8 9def generate_infill(xPhys, box_w, box_h, box_d, pattern_name,10                    period_mm=12.0, void_threshold=0.15, solid_threshold=0.75,11                    fine_nx=60, fine_ny=36, fine_nz=20):12    nely, nelx = xPhys.shape13    y_c = ((np.arange(nely) + 0.5) * box_h / nely)[::-1]14    x_c = (np.arange(nelx) + 0.5) * box_w / nelx15    interp = RegularGridInterpolator(16        (y_c[::-1], x_c), xPhys[::-1, :],17        method='linear', bounds_error=False, fill_value=0.018    )19    xs = np.linspace(0, box_w, fine_nx)20    ys = np.linspace(0, box_h, fine_ny)21    zs = np.linspace(0, box_d, fine_nz)22    X3, Y3, Z3 = np.meshgrid(xs, ys, zs, indexing='ij')23    dens = interp(np.stack([Y3.ravel(), X3.ravel()], axis=-1)).reshape(fine_nx, fine_ny, fine_nz)24 25    fn = INFILL_PATTERNS[pattern_name]["fn"]26    tpms = fn(X3, Y3, Z3, period_mm)27    tpms_norm = tpms / (np.abs(tpms).max() + 1e-9)28 29    sf = np.zeros_like(dens)30    vm = dens < void_threshold31    sm = dens > solid_threshold32    im = (~vm) & (~sm)33    sf[im] = tpms_norm[im] - (1.0 - 2.0 * dens[im])34    sf[sm] = 2.035    sf[vm] = -2.036 37    if 0. <= sf.min() or 0. >= sf.max():38        return None, 0, 0.0, sf39 40    try:41        v, f, _, _ = measure.marching_cubes(42            sf, level=0.0,43            spacing=(box_w / fine_nx, box_h / fine_ny, box_d / fine_nz)44        )45    except Exception:46        return None, 0, 0.0, sf47 48    stl_bytes = _triangles_to_stl(v, f)49    return stl_bytes, len(f), float((sf > 0).sum() / sf.size), sf50 51 52def to_stl_bytes(vol_or_xphys, box_w, box_h, box_d, iso, is_xphys=True):53    if is_xphys:54        nely, nelx = vol_or_xphys.shape55        n_depth = 856        # Build 3-D volume: axis-0=nely(Y), axis-1=nelx(X), axis-2=depth57        vol = np.zeros((nely, nelx, n_depth))58        for d in range(n_depth):59            vol[:, :, d] = vol_or_xphys60 61        # Pad with a zero-shell on all sides so that marching-cubes produces a62        # *closed* watertight surface even where the part touches the bounding box.63        # Without padding, those faces are open edges → Gmsh cannot tet-mesh them.64        vol_pad = np.pad(vol, 1, mode='constant', constant_values=0.0)65 66        # Spacing along (axis-0=Y, axis-1=X, axis-2=Z)67        dy, dx, dz = box_h / nely, box_w / nelx, box_d / n_depth68 69        if iso <= vol_pad.min() or iso >= vol_pad.max():70            return None, 071 72        try:73            verts, faces, _, _ = measure.marching_cubes(74                vol_pad, level=iso, spacing=(dy, dx, dz)75            )76        except Exception:77            return None, 078 79        # Subtract the 1-voxel padding offset so the part starts at (0,0,0)80        verts[:, 0] -= dy81        verts[:, 1] -= dx82        verts[:, 2] -= dz83        np.clip(verts[:, 0], 0, box_h, out=verts[:, 0])84        np.clip(verts[:, 1], 0, box_w, out=verts[:, 1])85        np.clip(verts[:, 2], 0, box_d, out=verts[:, 2])86 87        # marching_cubes returns (Y, X, Z). Swap to (X, Y, Z).88        # Swapping axes 0↔1 is a reflection (det=-1) → flip winding to preserve89        # outward normals so Gmsh classifySurfaces works correctly.90        verts = verts[:, [1, 0, 2]]91        faces = faces[:, [0, 2, 1]]92    else:93        vol = vol_or_xphys94        spacing = (box_w / vol.shape[0], box_h / vol.shape[1], box_d / vol.shape[2])95 96        if iso <= vol.min() or iso >= vol.max():97            return None, 098 99        try:100            verts, faces, _, _ = measure.marching_cubes(vol, level=iso, spacing=spacing)101        except Exception:102            return None, 0103 104    return _triangles_to_stl(verts, faces), len(faces)105 106 107def _triangles_to_stl(verts, faces):108    buf = io.BytesIO()109    buf.write(b'\x00' * 80)110    buf.write(struct.pack('<I', len(faces)))111    for tri in faces:112        v0, v1, v2 = verts[tri[0]], verts[tri[1]], verts[tri[2]]113        n = np.cross(v1 - v0, v2 - v0)114        nl = np.linalg.norm(n)115        n = (n / nl if nl > 0 else np.zeros(3)).astype(float)116        buf.write(struct.pack('<fff', *n))117        for vi in tri:118            buf.write(struct.pack('<fff', *verts[vi].astype(float)))119        buf.write(struct.pack('<H', 0))120    return buf.getvalue()121 122 123def voxelize_mesh(stl_bytes: bytes, nelx: int, nely: int, nelz: int) -> np.ndarray:124    """Voxelize an uploaded STL mesh onto a (nely, nelx, nelz) boolean grid.125 126    Returns a bool array where True = inside the design space.127    Falls back to all-True (full box) if trimesh is unavailable or mesh is not watertight.128    """129    try:130        import trimesh131        import trimesh.util132 133        mesh = trimesh.load(134            trimesh.util.wrap_as_stream(stl_bytes),135            file_type='stl',136            force='mesh'137        )138        if not isinstance(mesh, trimesh.Trimesh):139            raise ValueError("Uploaded file did not produce a single mesh.")140 141        pitch = max(mesh.extents) / max(nelx, nely, nelz)142        vox = mesh.voxelized(pitch=pitch).fill()143        raw = vox.matrix  # (nx, ny, nz) bool, may differ from target resolution144 145        # Resample to target resolution using nearest-neighbour146        from scipy.ndimage import zoom147        scale = (nely / raw.shape[0], nelx / raw.shape[1], nelz / raw.shape[2])148        resampled = zoom(raw.astype(float), scale, order=0) > 0.5149        # Ensure exact shape150        mask = np.zeros((nely, nelx, nelz), dtype=bool)151        sy = min(resampled.shape[0], nely)152        sx = min(resampled.shape[1], nelx)153        sz = min(resampled.shape[2], nelz)154        mask[:sy, :sx, :sz] = resampled[:sy, :sx, :sz]155        return mask156 157    except Exception:158        # Graceful fallback: treat entire bounding box as design space159        return np.ones((nely, nelx, nelz), dtype=bool)160