Implement Laplacian mesh smoothing (#975)

* Implement Laplacian mesh smoothing

Signed-off-by: Aritz Erkiaga <[email protected]>

avogadro/core/mesh.cpp

@@ -6,6 +6,7 @@
 #include "mesh.h"
 
 #include "mutex.h"
+#include "neighborperceiver.h"
 
 using std::vector;
 
@@ -177,5 +178,95 @@ Mesh& Mesh::operator=(const Mesh& other)
   return *this;
 }
 
+void Mesh::smooth(int iterationCount)
+{
+  if (m_vertices.size() == 0)
+    return;
+  if (iterationCount <= 0)
+    return;
+
+  // Map vertices to a line and pass them to NeighborPerceiver
+  Array<Vector3> linearList(m_vertices.size());
+  for (size_t i = 0; i < m_vertices.size(); i++)
+    // Empirical constants to make the distribution more homogeneous
+    linearList[i] = Vector3(
+      double(m_vertices[i](0) + 1.31*m_vertices[i](1) + 0.97*m_vertices[i](2)),
+    0.0, 0.0);
+  NeighborPerceiver perceiver(linearList, 0.01);
+
+  // Identify degenerate vertices
+  std::vector<int> indexToVertexID(m_vertices.size(), -1);
+  std::vector<std::vector<size_t>> vertexIDToIndices;
+  Array<size_t> neighbors;
+  for (size_t i = 0; i < m_vertices.size(); i++) {
+    if (indexToVertexID[i] != -1)
+      continue;
+    perceiver.getNeighborsInclusiveInPlace(neighbors, linearList[i]);
+    size_t vertexID = vertexIDToIndices.size();
+    for (size_t n: neighbors) {
+      if ((m_vertices[n] - m_vertices[i]).norm() < 0.0001) {
+        if (vertexID == vertexIDToIndices.size())
+          vertexIDToIndices.emplace_back();
+        indexToVertexID[n] = vertexID;
+        vertexIDToIndices[vertexID].push_back(n);
+      }
+    }
+  }
+
+  // Compute 1-ring
+  std::vector<std::vector<size_t>> vertexIDTo1Ring(vertexIDToIndices.size());
+  for (size_t id = 0; id < vertexIDToIndices.size(); id++) {
+    for (size_t v: vertexIDToIndices[id]) {
+      size_t relative = v % 3;
+      size_t triangle = v - relative;
+      std::array<size_t, 2> candidates{{
+        triangle + (relative + 1) % 3,
+        triangle + (relative + 2) % 3
+      }};
+      for (size_t candidate: candidates) {
+        size_t newID = indexToVertexID[candidate];
+        if (std::find(vertexIDToIndices[id].begin(), vertexIDToIndices[id].end(), newID)
+        == vertexIDToIndices[id].end())
+          vertexIDTo1Ring[id].push_back(newID);
+      }
+    }
+  }
+
+  float weight = 1.0f;
+  for (int iteration = 0; iteration < iterationCount; iteration++) {
+    // Copy vertices by ID into source array
+    std::vector<Vector3f> inputVertices(vertexIDToIndices.size());
+    for (size_t id = 0; id < vertexIDToIndices.size(); id++)
+      inputVertices[id] = m_vertices[vertexIDToIndices[id][0]];
+
+    // Apply Laplacian smoothing
+    for (size_t id = 0; id < inputVertices.size(); id++) {
+      Vector3f output(0.0f, 0.0f, 0.0f);
+      for (size_t neighbor: vertexIDTo1Ring[id])
+        output += inputVertices[neighbor];
+      output += weight * inputVertices[id];
+      output *= 1.0f / (weight + vertexIDTo1Ring[id].size());
+      for (size_t i: vertexIDToIndices[id])
+        m_vertices[i] = output;
+    }
+  }
+
+  // Recompute normals
+  for (size_t id = 0; id < vertexIDToIndices.size(); id++) {
+    Vector3f normal(0.0f, 0.0f, 0.0f);
+    for (size_t v: vertexIDToIndices[id]) {
+      size_t relative = v % 3;
+      size_t triangle = v - relative;
+      Vector3f &a = m_vertices[v];
+      Vector3f &b = m_vertices[triangle + (relative + 1) % 3];
+      Vector3f &c = m_vertices[triangle + (relative + 2) % 3];
+      Vector3f triangleNormal = (b - a).cross(c - a);
+      normal += triangleNormal.normalized();
+    }
+    for (size_t i: vertexIDToIndices[id])
+      m_normals[i] = normal.normalized();
+  }
+}
+
 } // End namespace QtGui
 } // End namespace Avogadro

avogadro/core/mesh.h

@@ -214,6 +214,12 @@ class AVOGADROCORE_EXPORT Mesh
    */
   Mutex* lock() const { return m_lock; }
 
+  /**
+   * Applies Laplacian smoothing.
+   * @param iterationCount number of smoothing passes to make.
+   */
+  void smooth(int iterationCount = 6);
+
   friend class Molecule;
 
 private:

avogadro/qtgui/meshgenerator.cpp

@@ -19,33 +19,34 @@ using Core::Cube;
 using Core::Mesh;
 
 MeshGenerator::MeshGenerator(QObject* p)
-  : QThread(p), m_iso(0.0), m_reverseWinding(false), m_cube(nullptr),
+  : QThread(p), m_iso(0.0), m_passes(6), m_reverseWinding(false), m_cube(nullptr),
     m_mesh(nullptr), m_stepSize(0.0, 0.0, 0.0), m_min(0.0, 0.0, 0.0),
     m_dim(0, 0, 0), m_progmin(0), m_progmax(0)
 {
 }
 
 MeshGenerator::MeshGenerator(const Cube* cube_, Mesh* mesh_, float iso,
-                             bool reverse, QObject* p)
-  : QThread(p), m_iso(0.0), m_reverseWinding(reverse), m_cube(nullptr),
+                             int passes, bool reverse, QObject* p)
+  : QThread(p), m_iso(0.0), m_passes(6), m_reverseWinding(reverse), m_cube(nullptr),
     m_mesh(nullptr), m_stepSize(0.0, 0.0, 0.0), m_min(0.0, 0.0, 0.0),
     m_dim(0, 0, 0), m_progmin(0), m_progmax(0)
 {
-  initialize(cube_, mesh_, iso);
+  initialize(cube_, mesh_, iso, passes);
 }
 
 MeshGenerator::~MeshGenerator()
 {
 }
 
 bool MeshGenerator::initialize(const Cube* cube_, Mesh* mesh_, float iso,
-                               bool reverse)
+                               int passes, bool reverse)
 {
   if (!cube_ || !mesh_)
     return false;
   m_cube = cube_;
   m_mesh = mesh_;
   m_iso = iso;
+  m_passes = passes;
   m_reverseWinding = reverse;
   if (!m_cube->lock()->tryLock()) {
     qDebug() << "Cannot get a read lock…";
@@ -102,11 +103,15 @@ void MeshGenerator::run()
   // Now we are done give all that memory back
   m_vertices.resize(0);
   m_normals.resize(0);
+
+  // Smooth out the mesh
+  m_mesh->smooth(m_passes);
 }
 
 void MeshGenerator::clear()
 {
   m_iso = 0.0;
+  m_passes = 6;
   m_cube = nullptr;
   m_mesh = nullptr;
   m_stepSize.setZero();

avogadro/qtgui/meshgenerator.h

@@ -53,10 +53,11 @@ class AVOGADROQTGUI_EXPORT MeshGenerator : public QThread
    * @param cube The source Cube with the volumetric data.
    * @param mesh The Mesh that will hold the isosurface.
    * @param iso The iso value of the surface.
+   * @param passes Number of smoothing passes to perform.
    * @return True if the MeshGenerator was successfully initialized.
    */
   MeshGenerator(const Core::Cube* cube, Core::Mesh* mesh, float iso,
-                bool reverse = false, QObject* parent = nullptr);
+                int passes = 6, bool reverse = false, QObject* parent = nullptr);
 
   /**
    * Destructor.
@@ -69,9 +70,10 @@ class AVOGADROQTGUI_EXPORT MeshGenerator : public QThread
    * @param cube The source Cube with the volumetric data.
    * @param mesh The Mesh that will hold the isosurface.
    * @param iso The iso value of the surface.
+   * @param passes Number of smoothing passes to perform.
    */
   bool initialize(const Core::Cube* cube, Core::Mesh* mesh, float iso,
-                  bool reverse = false);
+                  int passes = 6, bool reverse = false);
 
   /**
    * Use this function to begin Mesh generation. Uses an asynchronous thread,
@@ -135,10 +137,11 @@ class AVOGADROQTGUI_EXPORT MeshGenerator : public QThread
   bool marchingCube(const Vector3i& pos);
 
   float m_iso;              /** The value of the isosurface. */
-  bool m_reverseWinding;    /** Whether the winding and normals are reversed */
+  int m_passes;             /** Number of smoothing passes to perform. */
+  bool m_reverseWinding;    /** Whether the winding and normals are reversed. */
   const Core::Cube* m_cube; /** The cube that we are generating a Mesh from. */
   Core::Mesh* m_mesh;       /** The mesh that is being generated. */
-  Vector3f m_stepSize;      /** The step size vector for cube */
+  Vector3f m_stepSize;      /** The step size vector for cube. */
   Vector3f m_min;           /** The minimum point in the cube. */
   Vector3i m_dim;           /** The dimensions of the cube. */
   Core::Array<Vector3f> m_vertices, m_normals;

avogadro/qtplugins/surfaces/surfaces.cpp

@@ -471,7 +471,7 @@ void Surfaces::displayMesh()
     m_meshGenerator1 = new QtGui::MeshGenerator;
     connect(m_meshGenerator1, SIGNAL(finished()), SLOT(meshFinished()));
   }
-  m_meshGenerator1->initialize(m_cube, m_mesh1, -m_isoValue);
+  m_meshGenerator1->initialize(m_cube, m_mesh1, -m_isoValue, m_smoothingPasses);
 
   // TODO - only do this if we're generating an orbital
   //    and we need two meshes
@@ -482,7 +482,7 @@ void Surfaces::displayMesh()
     m_meshGenerator2 = new QtGui::MeshGenerator;
     connect(m_meshGenerator2, SIGNAL(finished()), SLOT(meshFinished()));
   }
-  m_meshGenerator2->initialize(m_cube, m_mesh2, m_isoValue, true);
+  m_meshGenerator2->initialize(m_cube, m_mesh2, m_isoValue, m_smoothingPasses, true);
 
   // Start the mesh generation - this needs an improved mutex with a read lock
   // to function as expected. Write locks are exclusive, read locks can have

avogadro/qtplugins/surfaces/surfaces.h

@@ -106,6 +106,7 @@ private slots:
   QtGui::MeshGenerator* m_meshGenerator2 = nullptr;
 
   float m_isoValue = 0.01;
+  int m_smoothingPasses = 6;
   int m_meshesLeft = 0;
 
   bool m_recordingMovie = false;