Overhaul GSoC stltovoxel work into a pure user_lib tool, i.e. not a gprMax command

conda_env.yml

@@ -17,9 +17,10 @@ dependencies:
 
 - pip:
 #  - mpi4py
+  - numpy-stl
 #  - pycuda
 #  - pyopencl
   - terminaltables
   - tqdm
   - git+https://github.com/craig-warren/PyEVTK.git
-  - numpy-stl
+  

gprMax/__init__.py

@@ -19,7 +19,6 @@ from .cmds_geometry.cylindrical_sector import CylindricalSector
 from .cmds_geometry.edge import Edge
 from .cmds_geometry.fractal_box import FractalBox
 from .cmds_geometry.geometry_objects_read import GeometryObjectsRead
-from .cmds_geometry.geometry_objects_read import GeometryObjectsReadSTL
 from .cmds_geometry.plate import Plate
 from .cmds_geometry.sphere import Sphere
 from .cmds_geometry.triangle import Triangle

gprMax/cmds_geometry/geometry_objects_read.py

@@ -29,58 +29,6 @@ from .cmds_geometry import UserObjectGeometry
 
 logger = logging.getLogger(__name__)
 
-import numpy as np
-import matplotlib.pyplot as plt
-from user_libs import stltovoxel
-
-class GeometryObjectsReadSTL(UserObjectGeometry):
-
-    def __init__(self, **kwargs):
-        super().__init__(**kwargs)
-        self.hash = '#geometry_objects_stl'
-    
-    def create(self, G, uip):
-
-        try:
-            stl_file = self.kwargs['stl_file']     # STL file of the model to be voxelized
-            mat_index = self.kwargs['mat_index']   #material index corresponding to the order of materials in the txt file
-            discretization = self.kwargs['discretization']  #spacial discretization of the model
-            
-        except KeyError:
-            logger.exception(self.__str__() + 'requires exactly three parameters')
-            raise
-
-        # See if stl file exists at specified path and if not try input
-        # file directory
-        
-        stl_file=Path(stl_file)
-        
-        if not stl_file.exists():
-            stl_file = Path(config.sim_config.input_file_path.parent, stl_file)
-        
-        model_array = stltovoxel.convert_file(stl_file,discretization)
-        print()
-        print("Resolution of the model along (X,Y,Z) axis is",model_array.shape)
-        parts = config.sim_config.input_file_path.with_suffix('').parts
-        filename_hdf5 = Path(*parts)
-        filename_hdf5 = filename_hdf5.with_suffix('.h5')
-
-        with h5py.File(filename_hdf5,'w') as hdf:
-            model_array[model_array==0] = -1
-            model_array[model_array==1] = mat_index
-            hdf.create_dataset('data',data=model_array)
-            hdf.attrs['dx_dy_dz']=(discretization[0],discretization[1],discretization[2])
-
-    def visualize(self,model_array):
-        #visualizer requires the model array to be built of only 0's and 1's
-        #if you want to use this function call it before the .h5 part
-        
-        from mpl_toolkits.mplot3d import Axes3D
-        fig = plt.figure()
-        ax = fig.add_subplot(111, projection='3d')
-        ax.voxels(model_array)
-        plt.show()
-
 
 class GeometryObjectsRead(UserObjectGeometry):
 
@@ -91,7 +39,7 @@ class GeometryObjectsRead(UserObjectGeometry):
     def rotate(self, axis, angle, origin=None):
         pass
 
-    def create(self, G, uip):
+    def create(self, grid, uip):
         """Create the object and add it to the grid."""
         try:
             p1 = self.kwargs['p1']
@@ -114,7 +62,7 @@ class GeometryObjectsRead(UserObjectGeometry):
             matfile = Path(config.sim_config.input_file_path.parent, matfile)
 
         matstr = matfile.with_suffix('').name
-        numexistmaterials = len(G.materials)
+        numexistmaterials = len(grid.materials)
 
         # Read materials from file
         with open(matfile, 'r') as f:
@@ -130,10 +78,10 @@ class GeometryObjectsRead(UserObjectGeometry):
             scene.add(material_obj)
 
         # Creates the internal simulation objects
-        scene.process_cmds(material_objs, G, sort=False)
+        scene.process_cmds(material_objs, grid, sort=False)
 
         # Update material type
-        for material in G.materials:
+        for material in grid.materials:
             if material.numID >= numexistmaterials:
                 if material.type:
                     material.type += ',\nimported'
@@ -148,7 +96,7 @@ class GeometryObjectsRead(UserObjectGeometry):
         # Open geometry object file and read/check spatial resolution attribute
         f = h5py.File(geofile, 'r')
         dx_dy_dz = f.attrs['dx_dy_dz']
-        if round_value(dx_dy_dz[0] / G.dx) != 1 or round_value(dx_dy_dz[1] / G.dy) != 1 or round_value(dx_dy_dz[2] / G.dz) != 1:
+        if round_value(dx_dy_dz[0] / grid.dx) != 1 or round_value(dx_dy_dz[1] / grid.dy) != 1 or round_value(dx_dy_dz[2] / grid.dz) != 1:
             logger.exception(self.__str__() + ' requires the spatial resolution of the geometry objects file to match the spatial resolution of the model')
             raise ValueError
 
@@ -165,12 +113,12 @@ class GeometryObjectsRead(UserObjectGeometry):
             rigidE = f['/rigidE'][:]
             rigidH = f['/rigidH'][:]
             ID = f['/ID'][:]
-            G.solid[xs:xs + data.shape[0], ys:ys + data.shape[1], zs:zs + data.shape[2]] = data + numexistmaterials
-            G.rigidE[:, xs:xs + rigidE.shape[1], ys:ys + rigidE.shape[2], zs:zs + rigidE.shape[3]] = rigidE
-            G.rigidH[:, xs:xs + rigidH.shape[1], ys:ys + rigidH.shape[2], zs:zs + rigidH.shape[3]] = rigidH
-            G.ID[:, xs:xs + ID.shape[1], ys:ys + ID.shape[2], zs:zs + ID.shape[3]] = ID + numexistmaterials
-            logger.info(self.grid_name(G) + f'Geometry objects from file {geofile} inserted at {xs * G.dx:g}m, {ys * G.dy:g}m, {zs * G.dz:g}m, with corresponding materials file {matfile}.')
+            grid.solid[xs:xs + data.shape[0], ys:ys + data.shape[1], zs:zs + data.shape[2]] = data + numexistmaterials
+            grid.rigidE[:, xs:xs + rigidE.shape[1], ys:ys + rigidE.shape[2], zs:zs + rigidE.shape[3]] = rigidE
+            grid.rigidH[:, xs:xs + rigidH.shape[1], ys:ys + rigidH.shape[2], zs:zs + rigidH.shape[3]] = rigidH
+            grid.ID[:, xs:xs + ID.shape[1], ys:ys + ID.shape[2], zs:zs + ID.shape[3]] = ID + numexistmaterials
+            logger.info(self.grid_name(grid) + f'Geometry objects from file {geofile} inserted at {xs * grid.dx:g}m, {ys * grid.dy:g}m, {zs * grid.dz:g}m, with corresponding materials file {matfile}.')
         except KeyError:
             averaging = False
-            build_voxels_from_array(xs, ys, zs, numexistmaterials, averaging, data, G.solid, G.rigidE, G.rigidH, G.ID)
-            logger.info(self.grid_name(G) + f'Geometry objects from file (voxels only){geofile} inserted at {xs * G.dx:g}m, {ys * G.dy:g}m, {zs * G.dz:g}m, with corresponding materials file {matfile}.')
+            build_voxels_from_array(xs, ys, zs, numexistmaterials, averaging, data, grid.solid, grid.rigidE, grid.rigidH, grid.ID)
+            logger.info(self.grid_name(grid) + f'Geometry objects from file (voxels only){geofile} inserted at {xs * grid.dx:g}m, {ys * grid.dy:g}m, {zs * grid.dz:g}m, with corresponding materials file {matfile}.')

user_libs/stltovoxel/__init__.py

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

user_libs/stltovoxel/convert.py

@@ -1,14 +1,18 @@
-from stl import mesh
 import numpy as np
+from stl import mesh
+
 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
 
 
@@ -18,12 +22,12 @@ def convert_file(input_file_path, discretization, pad=1, parallel=False):
 
 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)
+    vol = np.transpose(vol)
+    
     return vol
-
-

user_libs/stltovoxel/examples/bunny.py

@@ -0,0 +1,29 @@
+from pathlib import Path
+
+import gprMax
+
+fn = Path(__file__)
+
+title = gprMax.Title(name=fn.with_suffix('').name)
+domain = gprMax.Domain(p1=(0.148, 0.128, 0.148))
+dl = 0.001
+dxdydz = gprMax.Discretisation(p1=(dl, dl, dl))
+time_window = gprMax.TimeWindow(time=10e-9)
+
+g_read = gprMax.GeometryObjectsRead(p1=(0.020, 0.020, 0.020), geofile='bunny.h5', matfile= 'materials.txt')
+
+gv = gprMax.GeometryView(p1=(0, 0, 0),
+                         p2=domain.props.p1,
+                         dl=(dl, dl, dl),
+                         filename=fn.with_suffix('').name,
+                         output_type='n')
+                         
+scene = gprMax.Scene()
+scene.add(title)
+scene.add(domain)
+scene.add(dxdydz)
+scene.add(time_window)
+scene.add(g_read)
+scene.add(gv)
+
+gprMax.run(scenes=[scene], geometry_only=True, n=1, outputfile=fn)

user_libs/stltovoxel/perimeter.py

@@ -1,5 +1,6 @@
 from functools import reduce
 
+
 def lines_to_voxels(line_list, pixels):
     current_line_indices = set()
     x = 0
@@ -22,12 +23,14 @@ def slope_intercept(p1, p2):
     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
 
 
@@ -56,8 +59,10 @@ def paint_y_axis(lines, pixels, 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

@@ -1,8 +1,13 @@
-import numpy as np
 import multiprocessing as mp
+import sys
+
+import numpy as np
+from gprMax.utilities.utilities import get_terminal_width
 from tqdm import tqdm
+
 from . import perimeter
 
+
 def mesh_to_plane(mesh, bounding_box, parallel):
     if parallel:
         pool = mp.Pool(mp.cpu_count())
@@ -12,7 +17,7 @@ def mesh_to_plane(mesh, bounding_box, parallel):
 
     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:
+    with tqdm(total=bounding_box[2], desc="Processing Layers", ncols=get_terminal_width() - 1, file=sys.stdout) 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]
@@ -102,6 +107,7 @@ def where_line_crosses_z(p1, p2, z):
         distance = 0
     else:
         distance = (z - p1[2]) / (p2[2] - p1[2])
+    
     return linear_interpolation(p1, p2, distance)
 
 
@@ -121,8 +127,8 @@ def calculate_scale_shift(meshes, discretization):
     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
 
 
@@ -138,4 +144,5 @@ def generate_tri_events(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])

user_libs/stltovoxel/stltovoxel.py

@@ -0,0 +1,40 @@
+import argparse
+import logging
+from pathlib import Path
+
+import h5py
+import matplotlib.pyplot as plt
+
+from .convert import convert_file
+
+logger = logging.getLogger(__name__)
+logging.basicConfig(format='%(message)s', level=logging.INFO)
+
+if __name__ == "__main__":
+
+    # Parse command line arguments
+    parser = argparse.ArgumentParser(description='Allows the user to convert a STL files to voxelized mesh.', usage='cd gprMax; python -m user_libs.stltovoxel.stltovoxel stlfilename matindex dx_dy_dz')
+    parser.add_argument('stlfilename', help='name of STL file to convert including path')
+    parser.add_argument('-matindex', type=int, required=True,
+                        help='index of material to extract from STL file')
+    parser.add_argument('-dxdydz', nargs='+', type=float, required=True,
+                        help='discretisation to use in voxelisation process')
+    args = parser.parse_args()
+
+    filename_stl = Path(args.stlfilename)
+    dxdydz = tuple(args.dxdydz)
+
+    logger.info(f'\nConverting STL file: {filename_stl.name}')
+    model_array = convert_file(filename_stl, dxdydz)
+    model_array[model_array==0] = -1
+    model_array[model_array==1] = args.matindex
+    logger.info(f'Number of voxels: {model_array.shape[0]} x {model_array.shape[1]} x {model_array.shape[2]}')
+    logger.info(f'Spatial discretisation: {dxdydz[0]} x {dxdydz[1]} x {dxdydz[2]}m')
+    
+    # Write HDF5 file for gprMax using voxels
+    filename_hdf5 = filename_stl.with_suffix('.h5')
+    with h5py.File(filename_hdf5, 'w') as f:
+        f.create_dataset('data', data=model_array)
+        f.attrs['dx_dy_dz'] = (dxdydz[0], dxdydz[1], dxdydz[2])
+
+    logger.info(f'Written geometry object file: {filename_hdf5.name}')

user_models/bunny.py

@@ -1,34 +0,0 @@
-from pathlib import Path
-
-import gprMax
-
-fn = Path(__file__)
-
-title = gprMax.Title(name=fn.with_suffix('').name)
-domain = gprMax.Domain(p1=(0.12,0.1,0.12))
-dxdydz = gprMax.Discretisation(p1=(0.001, 0.001, 0.001))
-time_window = gprMax.TimeWindow(time=10e-9)
-
-g_stl = gprMax.GeometryObjectsReadSTL(stl_file='data_LasoMax/Stanford_Bunny.stl', mat_index=2, discretization=(0.001,0.001,0.001) )
-
-g_read = gprMax.GeometryObjectsRead(p1=(0 ,0 ,0), geofile=fn.with_suffix('.h5') , matfile= 'user_models/materials.txt')
-
-gv = gprMax.GeometryView(p1=(0, 0, 0),
-                         p2=(0.12, 0.1,0.12),
-                         dl=(0.001, 0.001, 0.001),
-                         filename=fn.with_suffix('').name,
-                         output_type='n')
-                         
-# create a scene
-scene = gprMax.Scene()
-# add the simulation objects to the scene
-scene.add(title)
-scene.add(domain)
-scene.add(dxdydz)
-scene.add(time_window)
-scene.add(g_stl)
-scene.add(g_read)
-scene.add(gv)
-
-# run the simulation
-gprMax.run(scenes=[scene], geometry_only=True, n=1, outputfile=fn)