Add automatic calculation of multiple isosurfaces for e.g. deformatio…

…n density (requires occ update) Bump OCC_VERSION string for updates
peterspackman · Oct 31, 2024 · 60ea5a1 · 60ea5a1
1 parent 528aa96
commit 60ea5a1
Show file tree

Hide file tree

Showing 16 changed files with 564 additions and 408 deletions.
diff --git a/resources/surface_description.json b/resources/surface_description.json
diff --git a/src/core/isosurface_parameters.cpp b/src/core/isosurface_parameters.cpp
@@ -10,7 +10,6 @@ void to_json(nlohmann::json &j,
        {"occName", f.occName},
        {"displayName", f.displayName},
        {"units", f.needsWavefunction},
-       {"needsIsovalue", f.needsIsovalue},
        {"needsOrbital", f.needsOrbital},
        {"description", f.description}};
 }
@@ -22,19 +21,15 @@ void from_json(const nlohmann::json &j,
   j.at("displayName").get_to(f.displayName);
   j.at("description").get_to(f.description);
 
-  if(j.contains("units")) {
+  if (j.contains("units")) {
     j.at("units").get_to(f.units);
   }
-  if(j.contains("needsIsovalue")) {
-    j.at("needsIsovalue").get_to(f.needsIsovalue);
-  }
-  if(j.contains("needsOrbital")) {
+  if (j.contains("needsOrbital")) {
     j.at("needsOrbital").get_to(f.needsOrbital);
   }
-  if(j.contains("needsWavefunction")) {
+  if (j.contains("needsWavefunction")) {
     j.at("needsWavefunction").get_to(f.needsWavefunction);
   }
-
 }
 
 void to_json(nlohmann::json &j, const isosurface::SurfaceDescription &f) {
@@ -47,37 +42,41 @@ void to_json(nlohmann::json &j, const isosurface::SurfaceDescription &f) {
        {"needsCluster", f.needsCluster},
        {"periodic", f.periodic},
        {"units", f.units},
-       {"requestableProperties", f.requestableProperties}};
+       {"requestableProperties", f.requestableProperties},
+       {"computeNegativeIsovalue", f.computeNegativeIsovalue}};
 }
 
 void from_json(const nlohmann::json &j, isosurface::SurfaceDescription &s) {
   j.at("displayName").get_to(s.displayName);
   j.at("occName").get_to(s.occName);
   j.at("description").get_to(s.description);
-  if(j.contains("needsIsovalue")) {
+  if (j.contains("needsIsovalue")) {
     j.at("needsIsovalue").get_to(s.needsIsovalue);
   }
-  if(j.contains("defaultIsovalue")) {
+  if (j.contains("defaultIsovalue")) {
     j.at("defaultIsovalue").get_to(s.defaultIsovalue);
   }
-  if(j.contains("needsWavefunction")) {
+  if (j.contains("needsWavefunction")) {
     j.at("needsWavefunction").get_to(s.needsWavefunction);
   }
-  if(j.contains("needsOrbital")) {
+  if (j.contains("needsOrbital")) {
     j.at("needsOrbital").get_to(s.needsOrbital);
   }
-  if(j.contains("needsCluster")) {
+  if (j.contains("needsCluster")) {
     j.at("needsCluster").get_to(s.needsCluster);
   }
-  if(j.contains("periodic")) {
+  if (j.contains("periodic")) {
     j.at("periodic").get_to(s.periodic);
   }
-  if(j.contains("units")) {
+  if (j.contains("units")) {
     j.at("units").get_to(s.units);
   }
-  if(j.contains("requestableProperties")) {
+  if (j.contains("requestableProperties")) {
     j.at("requestableProperties").get_to(s.requestableProperties);
   }
+  if (j.contains("computeNegativeIsovalue")) {
+    j.at("computeNegativeIsovalue").get_to(s.computeNegativeIsovalue);
+  }
 }
 
 namespace isosurface {
@@ -181,8 +180,9 @@ loadSurfaceDescriptions(const nlohmann::json &json) {
       SurfaceDescription sd;
       from_json(item.value(), sd);
       surfaces.insert(s(item.key()), sd);
-    } catch (nlohmann::json::exception& e) {
-      qWarning() << "Failed to parse surface" << s(item.key()) << ":" << e.what();
+    } catch (nlohmann::json::exception &e) {
+      qWarning() << "Failed to parse surface" << s(item.key()) << ":"
+                 << e.what();
     }
   }
   return surfaces;
@@ -192,7 +192,8 @@ QMap<QString, double> loadResolutionLevels(const nlohmann::json &json) {
   QMap<QString, double> resolutions;
   auto s = [](const std::string &str) { return QString::fromStdString(str); };
 
-  if (!json.contains("resolutionLevels") || !json["resolutionLevels"].is_object()) {
+  if (!json.contains("resolutionLevels") ||
+      !json["resolutionLevels"].is_object()) {
     qWarning() << "JSON does not contain a 'resolutions' object";
     return resolutions;
   }
@@ -201,8 +202,9 @@ QMap<QString, double> loadResolutionLevels(const nlohmann::json &json) {
     try {
       double value = item.value().get<double>();
       resolutions.insert(s(item.key()), value);
-    } catch (nlohmann::json::exception& e) {
-      qWarning() << "Failed to parse resolution" << s(item.key()) << ":" << e.what();
+    } catch (nlohmann::json::exception &e) {
+      qWarning() << "Failed to parse resolution" << s(item.key()) << ":"
+                 << e.what();
     }
   }
   return resolutions;

diff --git a/src/core/isosurface_parameters.h b/src/core/isosurface_parameters.h
@@ -70,6 +70,7 @@ struct Parameters {
   Kind kind;
   float isovalue{0.0};
   float separation{0.2};
+  bool computeNegativeIsovalue{false};
   ChemicalStructure *structure{nullptr};
   MolecularWavefunction *wfn{nullptr};
   Eigen::Isometry3d wfn_transform{Eigen::Isometry3d::Identity()};
@@ -91,7 +92,6 @@ struct SurfacePropertyDescription {
   QString displayName;
   QString units;
   bool needsWavefunction{false};
-  bool needsIsovalue{false};
   bool needsOrbital{false};
   QString description;
 };
@@ -105,6 +105,7 @@ struct SurfaceDescription {
   bool needsOrbital{false};
   bool needsCluster{false};
   bool periodic{false};
+  bool computeNegativeIsovalue{false};
   QString units{""};
   QString description{"Unknown"};
   QStringList requestableProperties{"none"};

diff --git a/src/core/mesh.cpp b/src/core/mesh.cpp
@@ -1,7 +1,7 @@
 #include "mesh.h"
-#include "meshinstance.h"
-#include "json.h"
 #include "eigen_json.h"
+#include "json.h"
+#include "meshinstance.h"
 #include <QFile>
 #include <QSignalBlocker>
 #include <fmt/os.h>
@@ -22,7 +22,6 @@ Mesh::Mesh(Eigen::Ref<const VertexList> vertices,
   m_vertexAreas = computeVertexAreas();
   m_vertexMask = Eigen::Matrix<bool, Eigen::Dynamic, 1>(m_vertices.cols());
   resetVertexMask(true);
-
 }
 
 Mesh::Mesh(Eigen::Ref<const VertexList> vertices, QObject *parent)
@@ -340,11 +339,12 @@ Mesh *Mesh::newFromJson(const nlohmann::json &object, QObject *parent) {
 
   // Parse properties
   const auto propertiesObj = object["properties"];
-  for (const auto &item: propertiesObj.items()) {
+  for (const auto &item : propertiesObj.items()) {
     const auto propertyValuesArray = item.value();
     ScalarPropertyValues propertyValues(propertyValuesArray.size());
     propertyValuesArray.get_to(propertyValues);
-    pointCloud->setVertexProperty(QString::fromStdString(item.key()), propertyValues);
+    pointCloud->setVertexProperty(QString::fromStdString(item.key()),
+                                  propertyValues);
   }
 
   return pointCloud;
@@ -374,12 +374,12 @@ bool Mesh::haveChildMatchingTransform(
     const Eigen::Isometry3d &transform) const {
   // TODO do this in bulk
 
-  auto check = [](occ::Vec3 a, occ::Vec3 b,
-                  double tolerance = 1e-1) {
-      return (a - b).norm() < tolerance;
+  auto check = [](occ::Vec3 a, occ::Vec3 b, double tolerance = 1e-1) {
+    return (a - b).norm() < tolerance;
   };
 
-  auto candidateCentroid = transform.rotation() * m_centroid + transform.translation();
+  auto candidateCentroid =
+      transform.rotation() * m_centroid + transform.translation();
 
   for (auto *child : children()) {
     auto *instance = qobject_cast<MeshInstance *>(child);
@@ -391,22 +391,99 @@ bool Mesh::haveChildMatchingTransform(
   return false;
 }
 
-void Mesh::resetFaceMask(bool value) {
-    m_faceMask.array() = value;
-}
+void Mesh::resetFaceMask(bool value) { m_faceMask.array() = value; }
 
-void Mesh::resetVertexMask(bool value) {
-    m_vertexMask.array() = value;
-}
+void Mesh::resetVertexMask(bool value) { m_vertexMask.array() = value; }
 
 ScalarPropertyValues Mesh::computeVertexAreas() const {
-  ScalarPropertyValues vertexAreas = ScalarPropertyValues::Zero(m_vertices.cols());
+  ScalarPropertyValues vertexAreas =
+      ScalarPropertyValues::Zero(m_vertices.cols());
 
   for (int i = 0; i < m_faces.cols(); ++i) {
-        float v = m_faceAreas(i) / 3.0f;
-        vertexAreas(m_faces(0, i)) += v;
-        vertexAreas(m_faces(1, i)) += v;
-        vertexAreas(m_faces(2, i)) += v;
+    float v = m_faceAreas(i) / 3.0f;
+    vertexAreas(m_faces(0, i)) += v;
+    vertexAreas(m_faces(1, i)) += v;
+    vertexAreas(m_faces(2, i)) += v;
   }
   return vertexAreas;
 }
+
+Mesh *Mesh::combine(const QList<Mesh *> &meshes) {
+  if (meshes.isEmpty()) {
+    return nullptr;
+  }
+
+  // Validate meshes are compatible
+  const auto &first = meshes.first();
+  for (const auto *mesh : meshes) {
+    if (!mesh)
+      continue;
+
+    if (mesh->kind() != first->kind()) {
+      qWarning() << "Cannot combine meshes of different kinds";
+      return nullptr;
+    }
+
+    if (mesh->atomsInside() != first->atomsInside() ||
+        mesh->atomsOutside() != first->atomsOutside()) {
+      qWarning() << "Cannot combine meshes with different atom configurations";
+      return nullptr;
+    }
+  }
+
+  int totalVertices = 0;
+  int totalFaces = 0;
+  for (const auto *mesh : meshes) {
+    if (!mesh)
+      continue;
+    totalVertices += mesh->numberOfVertices();
+    totalFaces += mesh->numberOfFaces();
+  }
+
+  // Prepare combined vertex and face matrices
+  VertexList combinedVertices(3, totalVertices);
+  FaceList combinedFaces(3, totalFaces);
+  ScalarPropertyValues isovalues(totalVertices);
+
+  // Copy data with offset tracking
+  int vertexOffset = 0;
+  int faceOffset = 0;
+
+  for (const auto *mesh : meshes) {
+    if (!mesh)
+      continue;
+
+    const int nVerts = mesh->numberOfVertices();
+    const int nFaces = mesh->numberOfFaces();
+
+    // Copy vertices
+    combinedVertices.block(0, vertexOffset, 3, nVerts) = mesh->vertices();
+
+    // Copy faces with offset
+    auto faces = mesh->faces();
+    faces.array() += vertexOffset;
+    combinedFaces.block(0, faceOffset, 3, nFaces) = faces;
+
+    // Set isovalues for this mesh's vertices
+    isovalues.segment(vertexOffset, nVerts).array() =
+        mesh->parameters().isovalue;
+
+    vertexOffset += nVerts;
+    faceOffset += nFaces;
+  }
+
+  // Create new mesh
+  auto *combinedMesh = new Mesh(combinedVertices, combinedFaces);
+
+  // Copy parameters from first mesh
+  combinedMesh->setParameters(first->parameters());
+
+  // Set the isovalue property
+  combinedMesh->setVertexProperty("isovalue", isovalues);
+
+  // Copy description
+  combinedMesh->setDescription(
+      QString("Combined mesh from %1 isosurfaces").arg(meshes.size()));
+
+  return combinedMesh;
+}
diff --git a/src/core/mesh.h b/src/core/mesh.h
@@ -140,6 +140,8 @@ class Mesh : public QObject {
   inline void highlightVertex(int v) { m_vertexHighlights.insert(v); }
   [[nodiscard]] const auto &vertexHighlights() const { return m_vertexHighlights; }
 
+  static Mesh * combine(const QList<Mesh*> &meshes);
+
 signals:
   void visibilityChanged();
   void transparencyChanged();
@@ -187,3 +189,4 @@ class Mesh : public QObject {
   ScalarPropertyValues m_emptyProperty;
   isosurface::Parameters m_params;
 };
+