Merge branch 'main' into cube-file-loading-and-viewing

README.md

@@ -24,6 +24,38 @@ To get new updates, simply run:
 ```
 git pull
 ```
+## Usage <a name=usage></a>
+TCviewer exports the `Screen` object which can be used to quickly build a scene. All the heavy work, such as setting up the camera, building and coloring meshes, will be done in the background. Please see the examples section for some simple example scripts for you to try out.
 
 ## Examples <a name=examples></a>
-[insert example outputs]
+To draw a molecule, simply run:
+
+```python
+from tcviewer import Screen
+from scm import plams
+
+molecule = 'path/to/molecule.xyz'  # molecules can be given as paths to xyz files
+molecule = plams.Molecule('path/to/molecule.xyz')  # or as plams.Molecule objects
+molecule.guess_bonds()
+with Screen() as scr:
+  scr.draw_molecule(molecule)
+```
+
+To load and draw the HOMO:
+
+```python
+from yutility import orbitals
+from tcviewer import Screen, materials
+import pathlib
+
+rkf_dir = pathlib.Path(__file__).parents[0]/'data'/'NH3BH3'
+
+# load orbitals and choose an MO to draw
+orbs = orbitals.Orbitals(rkf_dir/'adf.rkf')
+homo = orbs.mos['HOMO']
+# generate a cube file
+cub = homo.generate_orbital()
+
+with Screen() as scr:
+	scr.draw_cub(cub, 0.03, material=materials.orbital_matte)
+```

src/tcviewer/screen.py

@@ -153,16 +153,25 @@ def draw_isosurface(self, gridd, isovalue=0, color=None, material=None, with_pha
             with_phase = False
         if val.min() < isovalue < val.max():
             spacing = gridd.spacing if len(gridd.spacing) == 3 else gridd.spacing * 3
-            verts, faces, normals, values = skimage.measure.marching_cubes(val, isovalue, spacing=spacing, method='lorensen')
-            verts = o3d.cpu.pybind.utility.Vector3dVector(verts + gridd.origin)
-            triangles = o3d.cpu.pybind.utility.Vector3iVector(faces)
-
-            mesh = o3d.geometry.TriangleMesh(verts, triangles)
-            if color is not None:
-                mesh.paint_uniform_color(np.atleast_2d(color)[0])
-
-            mesh.compute_vertex_normals()
-            self.add_mesh(mesh, material=material)
+            verts, faces, normals, values = skimage.measure.marching_cubes(val, isovalue, spacing=spacing)
+
+            edges = []
+            for face in faces:
+                edges.extend(list(itertools.combinations(face, 2)))
+            g = nx.from_edgelist(edges)
+
+            # compute connected components and print results
+            components = list(nx.algorithms.components.connected_components(g))
+            # separate faces by component
+            for im, component in enumerate(components):
+                triangles = o3d.cpu.pybind.utility.Vector3iVector([face for face in faces if set(face) <= component])  # <= operator tests for subset relation
+                verts_ = o3d.cpu.pybind.utility.Vector3dVector(verts + gridd.origin)
+                mesh = o3d.geometry.TriangleMesh(verts_, triangles)
+                if color is not None:
+                    mesh.paint_uniform_color(np.atleast_2d(color)[0])
+
+                mesh.compute_vertex_normals()
+                self.add_mesh(mesh, material=material)
 
         if with_phase:
             self.draw_isosurface(gridd, isovalue=-isovalue, color=np.atleast_2d(color)[-1], material=material, with_phase=False)
@@ -262,4 +271,5 @@ def draw_axes(self, center=[0, 0, 0], length=1, width=.04, **kwargs):
         with Screen() as scr:
             scr.draw_molecule(mol)
             scr.draw_orbital(mo, material=materials.orbital_shiny, isovalue=.03)
+
             scr.draw_axes()