STL to Voxel converter

user_libs/stltovoxel/__init__.py

@@ -0,0 +1,7 @@
+from .main import convert_file, convert_files,convert_meshes
+
+__all__ = [
+    'convert_file',
+    'convert_files',
+    'convert_meshes',
+]

user_libs/stltovoxel/license.md

@@ -0,0 +1,21 @@
+The MIT License (MIT)
+
+Copyright (c) 2015 Christian Pederkoff
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.

user_libs/stltovoxel/main.py

@@ -0,0 +1,29 @@
+from stl import mesh
+import numpy as np
+from . import slice
+
+def convert_meshes(meshes, discretization, parallel=True):
+    scale, shift, shape = slice.calculate_scale_shift(meshes, discretization)
+    vol = np.zeros(shape[::-1], dtype=np.int32)
+    for mesh_ind, org_mesh in enumerate(meshes):
+        slice.scale_and_shift_mesh(org_mesh, scale, shift)
+        cur_vol = slice.mesh_to_plane(org_mesh, shape, parallel)
+        vol[cur_vol] = mesh_ind + 1
+    return vol, scale, shift
+
+
+def convert_file(input_file_path, discretization, pad=1, parallel=False):
+    return convert_files([input_file_path], discretization, pad=pad, parallel=parallel)
+
+
+def convert_files(input_file_paths, discretization, colors=[(0, 0, 0)], pad=1, parallel=False):
+    meshes = []
+    for input_file_path in input_file_paths:
+        mesh_obj = mesh.Mesh.from_file(input_file_path)
+        org_mesh = np.hstack((mesh_obj.v0[:, np.newaxis], mesh_obj.v1[:, np.newaxis], mesh_obj.v2[:, np.newaxis]))
+        meshes.append(org_mesh)
+    vol, scale, shift = convert_meshes(meshes, discretization, parallel)
+    vol=np.transpose(vol)
+    return vol
+
+

user_libs/stltovoxel/perimeter.py

@@ -0,0 +1,63 @@
+from functools import reduce
+
+def lines_to_voxels(line_list, pixels):
+    current_line_indices = set()
+    x = 0
+    for (event_x, status, line_ind) in generate_line_events(line_list):
+        while event_x - x >= 0:
+            lines = reduce(lambda acc, cur: acc + [line_list[cur]], current_line_indices, [])
+            paint_y_axis(lines, pixels, x)
+            x += 1
+
+        if status == 'start':
+            assert line_ind not in current_line_indices
+            current_line_indices.add(line_ind)
+        elif status == 'end':
+            assert line_ind in current_line_indices
+            current_line_indices.remove(line_ind)
+
+
+def slope_intercept(p1, p2):
+    x1, y1 = p1[:2]
+    x2, y2 = p2[:2]
+    slope = (y2 - y1) / (x2 - x1)
+    intercept = y1 - slope * x1
+    return slope, intercept
+
+
+def generate_y(p1, p2, x):
+    slope, intercept = slope_intercept(p1, p2)
+    y = slope * x + intercept
+    return y
+
+
+def paint_y_axis(lines, pixels, x):
+    is_black = False
+    target_ys = list(map(lambda line: int(generate_y(line[0], line[1], x)), lines))
+    target_ys.sort()
+    if len(target_ys) % 2:
+        distances = []
+        for i in range(len(target_ys) - 1):
+            distances.append(target_ys[i+1] - target_ys[i])
+        # https://stackoverflow.com/a/17952763
+        min_idx = -min((x, -i) for i, x in enumerate(distances))[1]
+        del target_ys[min_idx]
+
+    yi = 0
+    for target_y in target_ys:
+        if is_black:
+            # Bulk assign all pixels between yi -> target_y
+            pixels[yi:target_y, x] = True
+        pixels[target_y][x] = True
+        is_black = not is_black
+        yi = target_y
+    assert is_black is False, 'an error has occured at x%s' % x
+
+
+def generate_line_events(line_list):
+    events = []
+    for i, line in enumerate(line_list):
+        first, second = sorted(line, key=lambda pt: pt[0])
+        events.append((first[0], 'start', i))
+        events.append((second[0], 'end', i))
+    return sorted(events, key=lambda tup: tup[0])

user_libs/stltovoxel/slice.py

@@ -0,0 +1,141 @@
+import numpy as np
+import multiprocessing as mp
+from tqdm import tqdm
+from . import perimeter
+
+def mesh_to_plane(mesh, bounding_box, parallel):
+    if parallel:
+        pool = mp.Pool(mp.cpu_count())
+        result_ids = []
+
+    vol = np.zeros(bounding_box[::-1], dtype=bool)
+
+    current_mesh_indices = set()
+    z = 0
+    with tqdm(total=bounding_box[2], desc="Processing Layers", bar_format="{l_bar}{bar:50}{r_bar}{bar:-50b}") as pbar:
+        for event_z, status, tri_ind in generate_tri_events(mesh):
+            while event_z - z >= 0:
+                mesh_subset = [mesh[ind] for ind in current_mesh_indices]
+
+                if parallel:
+                    pass
+                else:
+                    pbar.update(1)
+                    _, pixels = paint_z_plane(mesh_subset, z, bounding_box[1::-1])
+                    vol[z]=pixels
+                z += 1
+
+            if status == 'start':
+                assert tri_ind not in current_mesh_indices
+                current_mesh_indices.add(tri_ind)
+            elif status == 'end':
+                assert tri_ind in current_mesh_indices
+                current_mesh_indices.remove(tri_ind)
+
+    if parallel:
+        results = [r.get() for r in result_ids]
+
+        for z, pixels in results:
+            vol[z] = pixels
+
+        pool.close()
+        pool.join()
+
+    return vol
+
+
+def paint_z_plane(mesh, height, plane_shape):
+    pixels = np.zeros(plane_shape, dtype=bool)
+
+    lines = []
+    for triangle in mesh:
+        triangle_to_intersecting_lines(triangle, height, pixels, lines)
+    perimeter.lines_to_voxels(lines, pixels)
+
+    return height, pixels
+
+
+def linear_interpolation(p1, p2, distance):
+    '''
+    :param p1: Point 1
+    :param p2: Point 2
+    :param distance: Between 0 and 1, Lower numbers return points closer to p1.
+    :return: A point on the line between p1 and p2
+    '''
+    return p1 * (1-distance) + p2 * distance
+
+
+def triangle_to_intersecting_lines(triangle, height, pixels, lines):
+    assert (len(triangle) == 3)
+    above = list(filter(lambda pt: pt[2] > height, triangle))
+    below = list(filter(lambda pt: pt[2] < height, triangle))
+    same = list(filter(lambda pt: pt[2] == height, triangle))
+    if len(same) == 3:
+        for i in range(0, len(same) - 1):
+            for j in range(i + 1, len(same)):
+                lines.append((same[i], same[j]))
+    elif len(same) == 2:
+        lines.append((same[0], same[1]))
+    elif len(same) == 1:
+        if above and below:
+            side1 = where_line_crosses_z(above[0], below[0], height)
+            lines.append((side1, same[0]))
+        else:
+            x = int(same[0][0])
+            y = int(same[0][1])
+            pixels[y][x] = True
+    else:
+        cross_lines = []
+        for a in above:
+            for b in below:
+                cross_lines.append((b, a))
+        side1 = where_line_crosses_z(cross_lines[0][0], cross_lines[0][1], height)
+        side2 = where_line_crosses_z(cross_lines[1][0], cross_lines[1][1], height)
+        lines.append((side1, side2))
+
+
+def where_line_crosses_z(p1, p2, z):
+    if (p1[2] > p2[2]):
+        p1, p2 = p2, p1
+    # now p1 is below p2 in z
+    if p2[2] == p1[2]:
+        distance = 0
+    else:
+        distance = (z - p1[2]) / (p2[2] - p1[2])
+    return linear_interpolation(p1, p2, distance)
+
+
+def calculate_scale_shift(meshes, discretization):
+    mesh_min = meshes[0].min(axis=(0, 1))
+    mesh_max = meshes[0].max(axis=(0, 1))
+   
+    for mesh in meshes[1:]:
+        mesh_min = np.minimum(mesh_min, mesh.min(axis=(0, 1)))
+        mesh_max = np.maximum(mesh_max, mesh.max(axis=(0, 1)))
+    amplitude = mesh_max - mesh_min 
+    #Standard Unit of STL is mm
+    vx=discretization[0]*1000 
+    vy=discretization[1]*1000
+    vz=discretization[2]*1000
+    
+    bx=int(amplitude[0]/vx)
+    by=int(amplitude[1]/vy)
+    bz=int(amplitude[2]/vz)
+    print(amplitude)
+    bounding_box = [bx+1, by+1, bz+1]
+    return max(1/vx,1/vy,1/vz), mesh_min, bounding_box
+
+
+def scale_and_shift_mesh(mesh, scale, shift):
+    for i, dim_shift in enumerate(shift):
+        mesh[..., i] = (mesh[..., i] - dim_shift) * scale
+
+
+def generate_tri_events(mesh):
+    # Create data structure for plane sweep
+    events = []
+    for i, tri in enumerate(mesh):
+        bottom, middle, top = sorted(tri, key=lambda pt: pt[2])
+        events.append((bottom[2], 'start', i))
+        events.append((top[2], 'end', i))
+    return sorted(events, key=lambda tup: tup[0])