From 722e9d54282befe9eac75b0cd6879a336fe0e394 Mon Sep 17 00:00:00 2001 From: meta-androcto Date: Thu, 15 Jun 2017 22:13:42 +1000 Subject: io_convert_image_to_mesh_img: complete rewrite: T51754 --- io_convert_image_to_mesh_img/mesh/__init__.py | 25 +++ io_convert_image_to_mesh_img/mesh/dtm.py | 219 ++++++++++++++++++++ io_convert_image_to_mesh_img/mesh/terrain.py | 245 +++++++++++++++++++++++ io_convert_image_to_mesh_img/mesh/triangulate.py | 186 +++++++++++++++++ 4 files changed, 675 insertions(+) create mode 100644 io_convert_image_to_mesh_img/mesh/__init__.py create mode 100644 io_convert_image_to_mesh_img/mesh/dtm.py create mode 100644 io_convert_image_to_mesh_img/mesh/terrain.py create mode 100644 io_convert_image_to_mesh_img/mesh/triangulate.py (limited to 'io_convert_image_to_mesh_img/mesh') diff --git a/io_convert_image_to_mesh_img/mesh/__init__.py b/io_convert_image_to_mesh_img/mesh/__init__.py new file mode 100644 index 00000000..edbc8c88 --- /dev/null +++ b/io_convert_image_to_mesh_img/mesh/__init__.py @@ -0,0 +1,25 @@ +# This file is a part of the HiRISE DTM Importer for Blender +# +# Copyright (C) 2017 Arizona Board of Regents on behalf of the Planetary Image +# Research Laboratory, Lunar and Planetary Laboratory at the University of +# Arizona. +# +# This program is free software: you can redistribute it and/or modify it +# under the terms of the GNU General Public License as published by the Free +# Software Foundation, either version 3 of the License, or (at your option) +# any later version. +# +# This program is distributed in the hope that it will be useful, but +# WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY +# or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +# for more details. +# +# You should have received a copy of the GNU General Public License along +# with this program. If not, see . + +"""A sub-package for loading DTMs as 3D models""" + +from . import dtm +from . import terrain + +__all__ = ['dtm', 'terrain', ] diff --git a/io_convert_image_to_mesh_img/mesh/dtm.py b/io_convert_image_to_mesh_img/mesh/dtm.py new file mode 100644 index 00000000..a6ab6e30 --- /dev/null +++ b/io_convert_image_to_mesh_img/mesh/dtm.py @@ -0,0 +1,219 @@ +# This file is a part of the HiRISE DTM Importer for Blender +# +# Copyright (C) 2017 Arizona Board of Regents on behalf of the Planetary Image +# Research Laboratory, Lunar and Planetary Laboratory at the University of +# Arizona. +# +# This program is free software: you can redistribute it and/or modify it +# under the terms of the GNU General Public License as published by the Free +# Software Foundation, either version 3 of the License, or (at your option) +# any later version. +# +# This program is distributed in the hope that it will be useful, but +# WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY +# or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +# for more details. +# +# You should have received a copy of the GNU General Public License along +# with this program. If not, see . + +"""Objects for importing HiRISE DTMs.""" + +import numpy as np + +from .. import pvl + + +class DTM: + """ + HiRISE Digital Terrain Model + + This class imports a HiRISE DTM from a Planetary Data Systems (PDS) + compliant .IMG file. + + Parameters + ---------- + path : str + terrain_resolution : float, optional + Controls the resolution the DTM is read at. This should be a float + in the range [0.01, 1.0] (and will be constrained to this range). A + value of 1.0 will result in the DTM being read at full resolution. A + value of 0.01 will result in the DTM being read at 1/100th resolution. + Default is 1.0 (no downsampling). + + Todo + ---- + * Use GDAL for importing the DTM if it is installed for this Python + environment. If/when I have the time to do this, it probably + warrants breaking out separate importer classes. The benefits of + doing this are pretty substantial, though: + + + More reliable (doesn't rely on my PVL parser for finding the + valid values in the DTM, for locating the starting position of + the elevation data in the .IMG file) + + + Other, better, downsampling algorithms are already built in. + + + Would make this much better at general PDS DTM importing, + currently some of the import code is specific to HiRISE DTMs. + + """ + + # Special constants in our data: + # NULL : No data at this point. + # LRS : Low Representation Saturation + # LIS : Low Instrument Saturation + # HRS : High Representation Saturation + # HIS : High Insturment Saturation + SPECIAL_VALUES = { + "NULL": np.fromstring(b'\xFF\x7F\xFF\xFB', dtype='>f4')[0], + "LRS": np.fromstring(b'\xFF\x7F\xFF\xFC', dtype='>f4')[0], + "LIS": np.fromstring(b'\xFF\x7F\xFF\xFD', dtype='>f4')[0], + "HRS": np.fromstring(b'\xFF\x7F\xFF\xFE', dtype='>f4')[0], + "HIS": np.fromstring(b'\xFF\x7F\xFF\xFF', dtype='>f4')[0] + } + + def __init__(self, path, terrain_resolution=1.0): + self.path = path + self.terrain_resolution = terrain_resolution + self.label = self._read_label() + self.data = self._read_data() + + def _read_label(self): + """Returns a dict-like representation of a PVL label""" + return pvl.load(self.path) + + def _read_data(self): + """ + Reads elevation data from a PDS .IMG file. + + Notes + ----- + * Uses nearest-neighbor to downsample data. + + Todo + ---- + * Add other downsampling algorithms. + + """ + h, w = self.image_resolution + max_samples = int(w - w % self.bin_size) + + data = np.zeros(self.shape) + with open(self.path, 'rb') as f: + # Seek to the first byte of data + start_byte = self._get_data_start() + f.seek(start_byte) + # Iterate over each row of the data + for r in range(data.shape[0]): + # Each iteration, seek to the right location before + # reading a row. We determine this location as the + # first byte of data PLUS a offset which we calculate as the + # product of: + # + # 4, the number of bytes in a single record + # r, the current row index + # w, the number of records in a row of the DTM + # bin_size, the number of records in a bin + # + # This is where we account for skipping over rows. + offset = int(4 * r * w * self.bin_size) + f.seek(start_byte + offset) + # Read a row + row = np.fromfile(f, dtype=np.float32, count=max_samples) + # This is where we account for skipping over columns. + data[r] = row[::self.bin_size] + + data = self._process_invalid_data(data) + return data + + def _get_data_start(self): + """Gets the start position of the DTM data block""" + label_length = self.label['RECORD_BYTES'] + num_labels = self.label.get('LABEL_RECORDS', 1) + return int(label_length * num_labels) + + def _process_invalid_data(self, data): + """Sets any 'NULL' elevation values to np.NaN""" + invalid_data_mask = (data <= self.SPECIAL_VALUES['NULL']) + data[invalid_data_mask] = np.NaN + return data + + @property + def map_size(self): + """Geographic size of the bounding box around the DTM""" + scale = self.map_scale * self.unit_scale + w = self.image_resolution[0] * scale + h = self.image_resolution[1] * scale + return (w, h) + + @property + def mesh_scale(self): + """Geographic spacing between mesh vertices""" + return self.bin_size * self.map_scale * self.unit_scale + + @property + def map_info(self): + """Map Projection metadata""" + return self.label['IMAGE_MAP_PROJECTION'] + + @property + def map_scale(self): + """Geographic spacing between DTM posts""" + map_scale = self.map_info.get('MAP_SCALE', None) + return getattr(map_scale, 'value', 1.0) + + @property + def map_units(self): + """Geographic unit for spacing between DTM posts""" + map_scale = self.map_info.get('MAP_SCALE', None) + return getattr(map_scale, 'units', None) + + @property + def unit_scale(self): + """ + The function that creates a Blender mesh from this object will assume + that the height values passed into it are in meters --- this + property is a multiplier to convert DTM-units to meters. + """ + scaling_factors = { + 'KM/PIXEL': 1000, + 'METERS/PIXEL': 1 + } + return scaling_factors.get(self.map_units, 1.0) + + @property + def terrain_resolution(self): + """Vertex spacing, meters""" + return self._terrain_resolution + + @terrain_resolution.setter + def terrain_resolution(self, t): + self._terrain_resolution = np.clip(t, 0.01, 1.0) + + @property + def bin_size(self): + """The width of the (square) downsampling bin""" + return int(np.ceil(1 / self.terrain_resolution)) + + @property + def image_stats(self): + """Image statistics from the original DTM label""" + return self.label['IMAGE'] + + @property + def image_resolution(self): + """(Line, Sample) resolution of the original DTM""" + return (self.image_stats['LINES'], self.image_stats['LINE_SAMPLES']) + + @property + def size(self): + """Number of posts in our reduced DTM""" + return self.shape[0] * self.shape[1] + + @property + def shape(self): + """Shape of our reduced DTM""" + num_rows = self.image_resolution[0] // self.bin_size + num_cols = self.image_resolution[1] // self.bin_size + return (num_rows, num_cols) diff --git a/io_convert_image_to_mesh_img/mesh/terrain.py b/io_convert_image_to_mesh_img/mesh/terrain.py new file mode 100644 index 00000000..db866289 --- /dev/null +++ b/io_convert_image_to_mesh_img/mesh/terrain.py @@ -0,0 +1,245 @@ +# This file is a part of the HiRISE DTM Importer for Blender +# +# Copyright (C) 2017 Arizona Board of Regents on behalf of the Planetary Image +# Research Laboratory, Lunar and Planetary Laboratory at the University of +# Arizona. +# +# This program is free software: you can redistribute it and/or modify it +# under the terms of the GNU General Public License as published by the Free +# Software Foundation, either version 3 of the License, or (at your option) +# any later version. +# +# This program is distributed in the hope that it will be useful, but +# WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY +# or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +# for more details. +# +# You should have received a copy of the GNU General Public License along +# with this program. If not, see . + +"""Objects for creating 3D models in Blender""" + +import bpy +import bmesh + +import numpy as np + +from .triangulate import Triangulate + + +class BTerrain: + """ + Functions for creating Blender meshes from DTM objects + + This class contains functions that convert DTM objects to Blender meshes. + Its main responsiblity is to triangulate a mesh from the elevation data in + the DTM. Additionally, it attaches some metadata to the object and creates + a UV map for it so that companion ortho-images drape properly. + + This class provides two public methods: `new()` and `reload()`. + + `new()` creates a new object[1] and attaches a new mesh to it. + + `reload()` replaces the mesh that is attached to an already existing + object. This allows us to retain the location and orientation of the parent + object's coordinate system but to reload the terrain at a different + resolution. + + Notes + ---------- + [1] If you're unfamiliar with Blender, one thing that will help you in + reading this code is knowing the difference between 'meshes' and + 'objects'. A mesh is just a collection of vertices, edges and + faces. An object may have a mesh as a child data object and + contains additional information, e.g. the location and orientation + of the coordinate system its child-meshes are reckoned in terms of. + + """ + + @staticmethod + def new(dtm, name='Terrain'): + """ + Loads a new terrain + + Parameters + ---------- + dtm : DTM + name : str, optional + The name that will be assigned to the new object, defaults + to 'Terrain' (and, if an object named 'Terrain' already + exists, Blender will automatically extend the name of the + new object to something like 'Terrain.001') + + Returns + ---------- + obj : bpy_types.Object + + """ + bpy.ops.object.add(type="MESH") + obj = bpy.context.object + obj.name = name + + # Fill the object data with a Terrain mesh + obj.data = BTerrain._mesh_from_dtm(dtm) + + # Add some meta-information to the object + metadata = BTerrain._create_metadata(dtm) + BTerrain._setobjattrs(obj, **metadata) + + # Center the mesh to its origin and create a UV map for draping + # ortho images. + BTerrain._center(obj) + + return obj + + @staticmethod + def reload(obj, dtm): + """ + Replaces an exisiting object's terrain mesh + + This replaces an object's mesh with a new mesh, transferring old + materials over to the new mesh. This is useful for reloading DTMs + at different resolutions but maintaining textures/location/rotation. + + Parameters + ----------- + obj : bpy_types.Object + An already existing Blender object + dtm : DTM + + Returns + ---------- + obj : bpy_types.Object + + """ + old_mesh = obj.data + new_mesh = BTerrain._mesh_from_dtm(dtm) + + # Copy any old materials to the new mesh + for mat in old_mesh.materials: + new_mesh.materials.append(mat.copy()) + + # Swap out the old mesh for the new one + obj.data = new_mesh + + # Update out-dated meta-information + metadata = BTerrain._create_metadata(dtm) + BTerrain._setobjattrs(obj, **metadata) + + # Center the mesh to its origin and create a UV map for draping + # ortho images. + BTerrain._center(obj) + + return obj + + @staticmethod + def _mesh_from_dtm(dtm, name='Terrain'): + """ + Creates a Blender *mesh* from a DTM + + Parameters + ---------- + dtm : DTM + name : str, optional + The name that will be assigned to the new mesh, defaults + to 'Terrain' (and, if an object named 'Terrain' already + exists, Blender will automatically extend the name of the + new object to something like 'Terrain.001') + + Returns + ---------- + mesh : bpy_types.Mesh + + Notes + ---------- + * We are switching coordinate systems from the NumPy to Blender. + + Numpy: Blender: + + ----> (0, j) ^ (0, y) + | | + | | + v (i, 0) + ----> (x, 0) + + """ + # Create an empty mesh + mesh = bpy.data.meshes.new(name) + + # Get the xy-coordinates from the DTM, see docstring notes + y, x = np.indices(dtm.data.shape).astype('float64') + x *= dtm.mesh_scale + y *= -1 * dtm.mesh_scale + + # Create an array of 3D vertices + vertices = np.dstack([x, y, dtm.data]).reshape((-1, 3)) + + # Drop vertices with NaN values (used in the DTM to represent + # areas with no data) + vertices = vertices[~np.isnan(vertices).any(axis=1)] + + # Calculate the faces of the mesh + triangulation = Triangulate(dtm.data) + faces = triangulation.face_list() + + # Fill the mesh + mesh.from_pydata(vertices, [], faces) + mesh.update() + + # Create a new UV layer + mesh.uv_textures.new("HiRISE Generated UV Map") + # We'll use a bmesh to populate the UV map with values + bm = bmesh.new() + bm.from_mesh(mesh) + bm.faces.ensure_lookup_table() + uv_layer = bm.loops.layers.uv[0] + + # Iterate over each face in the bmesh + num_faces = len(bm.faces) + w = dtm.data.shape[1] + h = dtm.data.shape[0] + for face_index in range(num_faces): + # Iterate over each loop in the face + for loop in bm.faces[face_index].loops: + # Get this loop's vertex coordinates + vert_coords = loop.vert.co.xy + # And calculate it's uv coordinate. We do this by dividing the + # vertice's x and y coordinates by: + # + # d + 1, dimensions of DTM (in "posts") + # mesh_scale, meters/DTM "post" + # + # This has the effect of mapping the vertex to its + # corresponding "post" index in the DTM, and then mapping + # that value to the range [0, 1). + u = vert_coords.x / ((w + 1) * dtm.mesh_scale) + v = 1 + vert_coords.y / ((h + 1) * dtm.mesh_scale) + loop[uv_layer].uv = (u, v) + + bm.to_mesh(mesh) + + return mesh + + @staticmethod + def _center(obj): + """Move object geometry to object origin""" + bpy.context.scene.objects.active = obj + bpy.ops.object.origin_set(center='BOUNDS') + + @staticmethod + def _setobjattrs(obj, **attrs): + for key, value in attrs.items(): + obj[key] = value + + @staticmethod + def _create_metadata(dtm): + """Returns a dict containing meta-information about a DTM""" + return { + 'PATH': dtm.path, + 'MESH_SCALE': dtm.mesh_scale, + 'DTM_RESOLUTION': dtm.terrain_resolution, + 'BIN_SIZE': dtm.bin_size, + 'MAP_SIZE': dtm.map_size, + 'MAP_SCALE': dtm.map_scale * dtm.unit_scale, + 'UNIT_SCALE': dtm.unit_scale, + 'IS_TERRAIN': True, + 'HAS_UV_MAP': True + } diff --git a/io_convert_image_to_mesh_img/mesh/triangulate.py b/io_convert_image_to_mesh_img/mesh/triangulate.py new file mode 100644 index 00000000..bb0a50ac --- /dev/null +++ b/io_convert_image_to_mesh_img/mesh/triangulate.py @@ -0,0 +1,186 @@ +# This file is a part of the HiRISE DTM Importer for Blender +# +# Copyright (C) 2017 Arizona Board of Regents on behalf of the Planetary Image +# Research Laboratory, Lunar and Planetary Laboratory at the University of +# Arizona. +# +# This program is free software: you can redistribute it and/or modify it +# under the terms of the GNU General Public License as published by the Free +# Software Foundation, either version 3 of the License, or (at your option) +# any later version. +# +# This program is distributed in the hope that it will be useful, but +# WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY +# or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +# for more details. +# +# You should have received a copy of the GNU General Public License along +# with this program. If not, see . + +"""Triangulation algorithms""" + +import numpy as np + + +class Triangulate: + """ + A triangulation algorithm for creating a mesh from a DTM raster. + + I have been re-writing parts of the Blender HiRISE DTM importer in an + effort to cull its dependencies on external packages. Originally, the + add-on relied on SciPy's Delaunay triangulation (really a wrapper for + Qhull's Delaunay triangulation) to triangulate a mesh from a HiRISE DTM. + + This re-write is much better suited to the problem domain. The SciPy + Delaunay triangulation creates a mesh from any arbitrary point cloud and, + while robust, doesn't care about the fact that our HiRISE DTMs are + regularly gridded rasters. This triangulation algorithm is less robust + but much faster. Credit is due to Tim Spriggs for his work on the previous + Blender HiRISE DTM importer --- this triangulation algorithm largely + models the one in his add-on with a few changes (namely interfacing + with NumPy's API). + + Overview + ---------- + Suppose we have a DTM: + + .. code:: + + - - - - - - - - X X - - - - - + - - - - - - X X X X X - - - - + - - - - X X X X X X X X - - - + - - X X X X X X X X X X X - - + X X X X X X X X X X X X X X - + - X X X X X X X X X X X X X X + - - X X X X X X X X X X X - - + - - - X X X X X X X X - - - - + - - - - X X X X X - - - - - - + - - - - - X X - - - - - - - - + + where 'X' represents valid values and '-' represents invalid values. + Valid values should become vertices in the resulting mesh, invalid + values should be ignored. + + Our end goal is to supply Blender with: + + 1. an (n x 3) list of vertices + + 2. an (m x 3) list of faces. + + A vertex is a 3-tuple that we get from the DTM raster array. The + z-coordinate is whatever elevation value is in the DTM and the xy- + coordinates are the image indices multiplied by the resolution of the + DTM (e.g. if the DTM is at 5m/px, the first vertex is at (0m, 0m, + z_00) and the vertex to the right of it is at (5m, 0m, z_01)). + + A face is a 3-tuple (because we're using triangles) where each element + is an index of a vertex in the vertices list. Computing the faces is + tricky because we want to leverage the orthogonal structure of the DTM + raster for computational efficiency but we also need to reference + vertex indices in our faces, which don't observe any regular + structure. + + We take two rows at a time from the DTM raster and track the *raster + row* indices as well as well as the *vertex* indices. Raster row + indices are the distance of a pixel in the raster from the left-most + (valid *or* invalid) pixel of the row. The first vertex is index 0 and + corresponds to the upperleft-most valid pixel in the DTM raster. + Vertex indices increase to the right and then down. + + For example, the first two rows: + + .. code:: + + - - - - - - - - X X - - - - - + - - - - - - X X X X X - - - - + + in vertex indices: + + .. code:: + + - - - - - - - - 0 1 - - - - - + - - - - - - 2 3 4 5 6 - - - - + + and in raster row indices: + + .. code:: + + - - - - - - - - 9 10 - - - - - + - - - - - - 7 8 9 10 11 - - - - + + To simplify, we will only add valid square regions to our mesh. So, + for these first two rows the only region that will be added to our + mesh is the quadrilateral formed by vertices 0, 1, 4 and 5. We + further divide this area into 2 triangles and add the vertices to the + face list in CCW order (i.e. t1: (4, 1, 0), t2: (4, 5, 1)). + + After the triangulation between two rows is completed, the bottom + row is cached as the top row and the next row in the DTM raster is + read as the new bottom row. This process continues until the entire + raster has been triangulated. + + Todo + --------- + * It should be pretty trivial to add support for triangular + regions (i.e. in the example above, also adding the triangles + formed by (3, 4, 0) and (5, 6, 1)). + + """ + def __init__(self, array): + self.array = array + self.faces = self._triangulate() + + def _triangulate(self): + """Triangulate a mesh from a topography array.""" + # Allocate memory for the triangles array + max_tris = (self.array.shape[0] - 1) * (self.array.shape[1] - 1) * 2 + tris = np.zeros((max_tris, 3), dtype=int) + ntri = 0 + + # We initialize a vertex counter at 0 + prev_vtx_start = 0 + # We don't care about the values in the array, just whether or not + # they are valid. + prev = ~np.isnan(self.array[0]) + # We can sum this boolean array to count the number of valid entries + prev_num_valid = prev.sum() + # TODO: Probably a more clear (and faster) function than argmax for + # getting the first Truth-y value in a 1d array. + prev_img_start = np.argmax(prev) + + # Start quadrangulation + for i in range(1, self.array.shape[0]): + # Fetch this row, get our bearings in image *and* vertex space + curr = ~np.isnan(self.array[i]) + curr_vtx_start = prev_vtx_start + prev_num_valid + curr_img_start = np.argmax(curr) + curr_num_valid = curr.sum() + # Find the overlap between this row and the previous one + overlap = np.logical_and(prev, curr) + num_tris = overlap.sum() - 1 + overlap_start = np.argmax(overlap) + # Store triangles + for j in range(num_tris): + curr_pad = overlap_start - curr_img_start + j + prev_pad = overlap_start - prev_img_start + j + tris[ntri + 0] = [ + curr_vtx_start + curr_pad, + prev_vtx_start + prev_pad + 1, + prev_vtx_start + prev_pad + ] + tris[ntri + 1] = [ + curr_vtx_start + curr_pad, + curr_vtx_start + curr_pad + 1, + prev_vtx_start + prev_pad + 1 + ] + ntri += 2 + # Cache current row as previous row + prev = curr + prev_vtx_start = curr_vtx_start + prev_img_start = curr_img_start + prev_num_valid = curr_num_valid + + return tris[:ntri] + + def face_list(self): + return list(self.faces) -- cgit v1.2.3