From b273acba928270961c125284674f2280eb4a2527 Mon Sep 17 00:00:00 2001
From: David Minton <daminton@purdue.edu>
Date: Sat, 19 Feb 2022 14:23:06 -0500
Subject: [PATCH] Updated the viewer to plot elevation points as square pixels
 as it is characterized internally

---
 python/ctem/ctem/ctem_viewer_3d.py     | 78 ++++++++++++++++++++------
 python/ctem/tests/viewer3d/polytest.py | 13 +++--
 2 files changed, 69 insertions(+), 22 deletions(-)

diff --git a/python/ctem/ctem/ctem_viewer_3d.py b/python/ctem/ctem/ctem_viewer_3d.py
index 4d545b69..093ae1e3 100644
--- a/python/ctem/ctem/ctem_viewer_3d.py
+++ b/python/ctem/ctem/ctem_viewer_3d.py
@@ -1,7 +1,7 @@
 import numpy as np
 import os
 import open3d
-import ctem_io_readers
+from ctem import ctem_io_readers
 import os
 
 
@@ -57,31 +57,76 @@ def ctem2mesh(self, surface_file):
         # Build mesh grid
         s = parameters['gridsize']
         pix = parameters['pix']
-        ygrid, xgrid = np.mgrid[-s / 2:s / 2, -s / 2:s / 2] * pix
+        ygrid, xgrid = np.mgrid[-s/2:s/2, -s/2:s/2] * pix
+        y0, x0 = ygrid - pix/2, xgrid - pix/2
+        y1, x1 = ygrid - pix/2, xgrid + pix/2
+        y2, x2 = ygrid + pix/2, xgrid + pix/2
+        y3, x3 = ygrid + pix/2, xgrid - pix/2
+
     
-        xvals = xgrid.flatten()
-        yvals = ygrid.flatten()
-        zvals = self.dem.flatten()
+        xvals = np.append(
+                np.append(
+                np.append(x0.flatten(),
+                          x1.flatten()),
+                          x2.flatten()),
+                          x3.flatten())
+        yvals = np.append(
+                np.append(
+                np.append(y0.flatten(),
+                          y1.flatten()),
+                          y2.flatten()),
+                          y3.flatten())
+        zvals = np.append(
+                np.append(
+                np.append(self.dem.flatten(),
+                          self.dem.flatten()),
+                          self.dem.flatten()),
+                          self.dem.flatten())
         verts = np.array((xvals, yvals, zvals)).T
-        faces = np.empty([2 * (s - 1) ** 2, 3], dtype=np.int64)
-        for j in range(s - 1):
-            for i in range(s - 1):
-                idx = (s - 1) * j + i
-                # faces[idx,:] = [ j*s+i, j*s+i+1, (j+1)*s+i+1, (j+1)*s+i ]
-                faces[idx, :] = [j * s + i, j * s + i + 1, (j + 1) * s + i + 1]
-                idx += (s - 1) ** 2
-                faces[idx, :] = [(j + 1) * s + i + 1, (j + 1) * s + i, j * s + i ]
+        nface_triangles = 10
+        faces = np.full([nface_triangles*s**2, 3], -1, dtype=np.int64)
+        for j in range(s):
+            for i in range(s):
+                i0 = s*j + i
+                i1 = i0 + s**2
+                i2 = i1 + s**2
+                i3 = i2 + s**2
+
+                fidx = np.arange(nface_triangles*i0,nface_triangles*(i0+1))
+                # Make the two top triangles
+                faces[fidx[0],:] = [i0, i1, i2]
+                faces[fidx[1],:] = [i0, i2, i3]
+                # Make the two west side triangles
+                if i > 0:
+                    faces[fidx[0],:] = [i0, i3, i3-1]
+                    faces[fidx[1],:] = [i0, i3-1, i0-1]
+                # Make the two south side triangles
+                if j > 0:
+                    faces[fidx[4],:] = [i1, i0, i3-s ]
+                    faces[fidx[5],:] = [i1, i3-s, i2-s]
+                # Make the two east side triangles
+                if i < (s - 1):
+                    faces[fidx[6],:] = [i2, i1, i0+1]
+                    faces[fidx[7],:] = [i2, i0+1, i3+1]
+                #Make the two north side triangles
+                if j < (s -1):
+                    faces[fidx[8],:] = [i3, i2, i1+s]
+                    faces[fidx[9],:] = [i3, i1+s, i0+s]
+        nz = faces[:,0] != -1
+        f2 = faces[nz,:]
         self.mesh.vertices = open3d.utility.Vector3dVector(verts)
-        self.mesh.triangles = open3d.utility.Vector3iVector(faces)
+        self.mesh.triangles = open3d.utility.Vector3iVector(f2)
         self.mesh.compute_vertex_normals()
-        
+        self.mesh.compute_triangle_normals()
+
         return
     
     def visualize(self):
         vis = open3d.visualization.Visualizer()
         vis.create_window()
         vis.add_geometry(self.mesh)
-        
+        opt = vis.get_render_option()
+        opt.background_color = np.asarray([0, 0, 0])
 
         # zmax = np.amax(surf.dem)
         # zmin = np.amin(surf.dem)
@@ -99,6 +144,7 @@ def visualize(self):
     surf = polysurface()
     surf.ctem2mesh(surface_file="surface_dem.dat")
     surf.visualize()
+    #open3d.visualization.draw_geometries_with_editing([surf.mesh])
 
     
 
diff --git a/python/ctem/tests/viewer3d/polytest.py b/python/ctem/tests/viewer3d/polytest.py
index 13cee521..c528374c 100644
--- a/python/ctem/tests/viewer3d/polytest.py
+++ b/python/ctem/tests/viewer3d/polytest.py
@@ -1,8 +1,9 @@
 import numpy as np
 from ctem import ctem_viewer_3d
-
-surf = ctem_viewer_3d.surf()
-surf.io_read_ply("216kleopatra-mesh.ply")
-v = np.asarray(surf.mesh)
-t = np.asarray(surf.mesh)
-surf.render()
\ No newline at end of file
+import open3d
+surf = ctem_viewer_3d.polysurface()
+#surf.io_read_ply("216kleopatra-mesh.ply")
+surf.mesh = open3d.geometry.TriangleMesh.create_box()
+v = np.asarray(surf.mesh.vertices)
+t = np.asarray(surf.mesh.triangles)
+surf.visualize()
\ No newline at end of file