Merge pull request #295 from karban8/devel

[GSoC 2021] Importing Geometrical Information from laser scanners

conda_env.yml

@@ -22,3 +22,4 @@ dependencies:
   - terminaltables
   - tqdm
   - git+https://github.com/craig-warren/PyEVTK.git
+  - numpy-stl

gprMax/__init__.py

@@ -19,6 +19,7 @@ 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,6 +29,58 @@ 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):
 
@@ -117,8 +169,8 @@ class GeometryObjectsRead(UserObjectGeometry):
             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(grid) + 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}.')
+            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}.')
         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(grid) + 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}.')
+            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}.')

user_libs/docs_lasoMax.ipynb

@@ -0,0 +1,163 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Importing Geometrical Information from Laser Scanners-\n",
+    "\n",
+    "The module focuses on the ability to directly model real objects without having to enter their geometries manually. A laser scanner takes distance measurements in every direction to rapidly capture the surface shape of objects, buildings and landscapes.The data acquired from the laser scanner is in the form of point cloud. This information is then used to construct a full 3D model of the object in the form of a STL mesh. The point cloud data can be processed through various free mesh processing softwares (for example, MeshLab) to be converted into STL format. The STL mesh is then voxelized. A voxel is a unit of graphic information that defines a point in three-dimensional space. The voxelized model array is further used to build FDTD cells upon which gprMax is based. "
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### GeometryObjectsReadSTL\n",
+    "\n",
+    "This function allows the user to convert the STL models to voxelized mesh. The input parameter requires the user to specify the path of the STL file, a material index (an integer corresponding to the materials given in the materials.txt file that you want the model to be built of) and also the spacial discretization/resolution of the model.\n",
+    "<br>\n",
+    "\n",
+    "Usage:\n",
+    "<br>\n",
+    "g_stl = gprMax.GeometryObjectsReadSTL(stl_file='your_model_file.stl', mat_index = n, discretization=(0.002,0.002,0.002))\n",
+    "<br>\n",
+    "where 'n' is an integer corresponding to the material\n",
+    "<br>\n",
+    " \n",
+    "The voxelized mesh is exported in the form of a numpy 3D array and further written to a HDF5 file which is outputted. \n",
+    "The HDF5 file along with the materials.txt file is given as input to the GeometryObjectsRead command which then builds the FDTD cells. The GeomteryView command is then used to output the model geometry which can be analysed/viewed in Paraview.\n",
+    "<br>\n",
+    "\n",
+    "#### NOTE:\n",
+    " - Before using this module make sure that your conda env has 'numpy-stl' pip dependency installed. You can check that by running 'conda list'. If not installed, then head over to the conda_env.yml file, under the pip section add 'numpy-stl' and rebuild the environment by using the following command\n",
+    " <br>\n",
+    " $ conda env update -f conda_env.yml\n",
+    "\n",
+    "#### Resolution\n",
+    "\n",
+    "Resolution is basically the no. of layers in which the model is divided. It can have different values for each axis, and these values are auto-calculated based on the spacial discretization entered by the user. A simple relation follows: smaller the spacial discretization, greater would be the resolution, thus greater would be the no. of voxels.\n",
+    "\n",
+    "\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "| STL Model of Mountain Topography |\n",
+    "|:----------------------------------|\n",
+    "|<img src='docs_laso_images/mont_blanc_stl.png' height=\"250\" width=\"350\">|"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "| Voxelized Model of Topography with resolution 1mm | Voxelized Model of Topography with resolution 2mm |\n",
+    "|:---------------------------------------------------|:---------------------------------------------------|\n",
+    "|<img src='docs_laso_images/mont_blanc.png' height=\"300\" width=\"370\">|<img src='docs_laso_images/mont_blanc_2.png' height=\"300\" width=\"370\">|"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "#Sample input file of Mont Blanc\n",
+    "from pathlib import Path\n",
+    "import gprMax\n",
+    "\n",
+    "basefile = Path().resolve().parents[0]\n",
+    "fn = basefile/'user_models'/'Mont_Blanc.py'\n",
+    "\n",
+    "title = gprMax.Title(name=fn.with_suffix('').name)\n",
+    "domain = gprMax.Domain(p1=(0.4,0.4 ,0.12))\n",
+    "dxdydz = gprMax.Discretisation(p1=(0.002, 0.002, 0.002))\n",
+    "time_window = gprMax.TimeWindow(time=10e-9)\n",
+    "\n",
+    "g_stl = gprMax.GeometryObjectsReadSTL(stl_file = basefile/'data_LasoMax'/'Mont_Blanc.stl', mat_index=2, discretization=(0.002,0.002,0.002))\n",
+    "\n",
+    "g_read = gprMax.GeometryObjectsRead(p1=(0.04,0.04,0.04), geofile=fn.with_suffix('.h5') , matfile=basefile/'user_models'/'materials.txt')\n",
+    "\n",
+    "gv = gprMax.GeometryView(p1=(0, 0, 0),\n",
+    "                         p2=(0.4, 0.4 , 0.12),\n",
+    "                         dl=(0.002, 0.002, 0.002),\n",
+    "                         filename=fn.with_suffix('').name,\n",
+    "                         output_type='n')\n",
+    "                         \n",
+    "# create a scene\n",
+    "scene = gprMax.Scene()\n",
+    "# add the simulation objects to the scene\n",
+    "scene.add(title)\n",
+    "scene.add(domain)\n",
+    "scene.add(dxdydz)\n",
+    "scene.add(time_window)\n",
+    "scene.add(g_stl)\n",
+    "scene.add(g_read)\n",
+    "scene.add(gv)\n",
+    "\n",
+    "# run the simulation\n",
+    "gprMax.run(scenes=[scene], geometry_only=True, n=1, outputfile=fn)\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "|Some example STL models | Corresponding Voxelized Models |\n",
+    "|:------------------------|:--------------------------------|\n",
+    "|<img src='docs_laso_images/bunny_stl.png' height=\"250\" width=\"300\">|<img src='docs_laso_images/bunny_voxel.png' height=\"250\" width=\"300\">|\n",
+    "|<img src='docs_laso_images/amphora_stl.png' height=\"300\" width=\"300\">|<img src='docs_laso_images/amphora.png' height=\"300\" width=\"300\">|\n",
+    "|<img src='docs_laso_images/house_stl.png' height=\"250\" width=\"300\">|<img src='docs_laso_images/house.png' height=\"250\" width=\"300\">|\n",
+    "|<img src='docs_laso_images/saint_michel_island_stl.png' height=\"250\" width=\"300\">|<img src='docs_laso_images/saint_michel_island.png' height=\"250\" width=\"300\">|"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Point Cloud to STL\n",
+    "\n",
+    "There are various free point cloud and mesh processing softwares available to convert laser point cloud data (in formats like .ply, .pcd, .xyz etc) to STL. Meshlab, which is one such software, is being used here. To convert your Point Cloud to STL,\n",
+    "firstly import your Point Cloud Data into Meshlab. If you feel the need you can scale your Point Cloud by heading over to \n",
+    "Filter -> Normal and Curvature -> Scale and scale it accordingly.\n",
+    "If you wish to reduce the density of points in the cloud then under Remeshing use 'Quadric Edge Collapse Decimation'.\n",
+    "Finally export the mesh as STL. Now you are ready to import this STL file into gprMax using the GeometryObjectsReadSTL function."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "| Point Cloud | Corresponding STL Mesh |\n",
+    "|:-------------|:------------------------|\n",
+    "|<img src='docs_laso_images/alps_point_cloud.png' height=\"250\" width=\"300\">|<img src='docs_laso_images/alp_stl.png' height=\"250\" width=\"300\">|"
+   ]
+  }
+ ],
+ "metadata": {
+  "interpreter": {
+   "hash": "65a73e13105b4b4b908542c2acdabd574f39be4b6464136c25622ad1e7e9b7ba"
+  },
+  "kernelspec": {
+   "display_name": "Python 3.9.6 64-bit ('gprMax': conda)",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.9.6"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 4
+}

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])

user_models/Mont_Blanc.py

@@ -0,0 +1,34 @@
+from pathlib import Path
+
+import gprMax
+
+fn = Path(__file__)
+
+title = gprMax.Title(name=fn.with_suffix('').name)
+domain = gprMax.Domain(p1=(0.4,0.4 ,0.12))
+dxdydz = gprMax.Discretisation(p1=(0.002, 0.002, 0.002))
+time_window = gprMax.TimeWindow(time=10e-9)
+
+g_stl = gprMax.GeometryObjectsReadSTL(stl_file='data_LasoMax/Mont_Blanc.stl', mat_index=2, discretization=(0.002,0.002,0.002))
+
+g_read = gprMax.GeometryObjectsRead(p1=(0.04,0.04,0.04), geofile=fn.with_suffix('.h5') , matfile='user_models/materials.txt')
+
+gv = gprMax.GeometryView(p1=(0, 0, 0),
+                         p2=(0.4, 0.4 , 0.12),
+                         dl=(0.002, 0.002, 0.002),
+                         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)

user_models/amphora.py

@@ -0,0 +1,34 @@
+from pathlib import Path
+
+import gprMax
+
+fn = Path(__file__)
+
+title = gprMax.Title(name=fn.with_suffix('').name)
+domain = gprMax.Domain(p1=(0.08,0.07,0.11))
+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/amphora.stl', mat_index=2, discretization =(0.001,0.001,0.001))
+
+g_read = gprMax.GeometryObjectsRead(p1=(0.001,0.001,0.0), geofile=fn.with_suffix('.h5') , matfile= 'user_models/materials.txt')
+
+gv = gprMax.GeometryView(p1=(0, 0, 0),
+                         p2=(0.08, 0.07,0.11),
+                         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)

user_models/bunny.py

@@ -0,0 +1,34 @@
+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)

user_models/dubai.py

@@ -0,0 +1,34 @@
+from pathlib import Path
+
+import gprMax
+
+fn = Path(__file__)
+
+title = gprMax.Title(name=fn.with_suffix('').name)
+domain = gprMax.Domain(p1=(0.14,0.12,0.12))
+dxdydz = gprMax.Discretisation(p1=(0.0008, 0.0008, 0.0008))
+time_window = gprMax.TimeWindow(time=10e-9)
+
+g_stl = gprMax.GeometryObjectsReadSTL(stl_file='data_LasoMax/Dubai.stl', mat_index=1, discretization=(0.0008,0.0008,0.0008) )
+
+g_read = gprMax.GeometryObjectsRead(p1=(0.01 ,0.01,0.01), geofile=fn.with_suffix('.h5') , matfile= 'user_models/materials.txt')
+
+gv = gprMax.GeometryView(p1=(0, 0, 0),
+                         p2=(0.14, 0.12,0.12),
+                         dl=(0.0008, 0.0008, 0.0008),
+                         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)

user_models/house.py

@@ -0,0 +1,34 @@
+from pathlib import Path
+
+import gprMax
+
+fn = Path(__file__)
+
+title = gprMax.Title(name=fn.with_suffix('').name)
+domain = gprMax.Domain(p1=(0.112,0.104,0.082))
+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/house.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.112, 0.104,0.082),
+                         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)

user_models/materials.txt

@@ -0,0 +1,4 @@
+#material: 8 0 1 0 pseudo_mat
+#material: 9 0 1 0 mat1
+#material: 8.5 0 1 0 mat2
+#material: 7 0 1 0 mat3

user_models/newyorkcity.py

@@ -0,0 +1,34 @@
+from pathlib import Path
+
+import gprMax
+
+fn = Path(__file__)
+
+title = gprMax.Title(name=fn.with_suffix('').name)
+domain = gprMax.Domain(p1=(0.11,0.115,0.07))
+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/newyorkcity.stl', mat_index=2, discretization=(0.001,0.001,0.001))
+
+g_read = gprMax.GeometryObjectsRead(p1=(0.001,0.001,0.0), geofile=fn.with_suffix('.h5') , matfile= 'user_models/materials.txt')
+
+gv = gprMax.GeometryView(p1=(0, 0, 0),
+                         p2=(0.11, 0.115,0.07),
+                         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)