Use flying edges for mesh generation

avogadro/core/cube.cpp

@@ -7,6 +7,7 @@
 
 #include "molecule.h"
 #include "mutex.h"
+#include <iostream>
 
 namespace Avogadro::Core {
 
@@ -23,21 +24,32 @@ Cube::~Cube()
   m_lock = nullptr;
 }
 
+std::vector<float>::const_iterator Cube::getRowIter(size_t j, size_t k) const
+{
+    return m_data.cbegin() + m_points.x() * (k * m_points.y() + j);
+}
+
+
 bool Cube::setLimits(const Vector3& min_, const Vector3& max_,
                      const Vector3i& points)
 {
-  // We can calculate all necessary properties and initialise our data
+  // Calculate the maximum dimension
+  int maxPoints = std::max({points.x(), points.y(), points.z()});
+  Vector3i equalPoints(maxPoints, maxPoints, maxPoints);
+
+  // Recalculate spacing based on equal points
   Vector3 delta = max_ - min_;
-  m_spacing =
-    Vector3(delta.x() / (points.x() - 1), delta.y() / (points.y() - 1),
-            delta.z() / (points.z() - 1));
+  m_spacing = Vector3(delta.x() / (equalPoints.x() - 1),
+                      delta.y() / (equalPoints.y() - 1),
+                      delta.z() / (equalPoints.z() - 1));
   m_min = min_;
   m_max = max_;
-  m_points = points;
+  m_points = equalPoints;
   m_data.resize(m_points.x() * m_points.y() * m_points.z());
   return true;
 }
 
+
 bool Cube::setLimits(const Vector3& min_, const Vector3& max_, float spacing_)
 {
   Vector3 delta = max_ - min_;
@@ -194,6 +206,168 @@ float Cube::value(int i, int j, int k) const
     return 0.0;
 }
 
+std::array<float, 3> Cube::computeGradient(size_t i, size_t j, size_t k) const
+{
+    std::array<std::array<float, 2>, 3> x;
+    std::array<float, 3> run;
+
+    size_t dataIdx = i + (j * m_points.x()) + (k * m_points.x() * m_points.y());
+
+        // std::cout << m_spacing[0];
+    if (i == 0)
+    {
+        x[0][0] = m_data[dataIdx + 1];
+        x[0][1] = m_data[dataIdx];
+        run[0] = m_spacing.x();
+    }
+    else if (i == (m_points.x() - 1))
+    {
+        x[0][0] = m_data[dataIdx];
+        x[0][1] = m_data[dataIdx - 1];
+        run[0] = m_spacing.x();
+    }
+    else
+    {
+        x[0][0] = m_data[dataIdx + 1];
+        x[0][1] = m_data[dataIdx - 1];
+        run[0] = 2 * m_spacing.x();
+    }
+
+    if (j == 0)
+    {
+        x[1][0] = m_data[dataIdx + m_points.x()];
+        x[1][1] = m_data[dataIdx];
+        run[1] = m_spacing.y();
+    }
+    else if (j == (m_points.y() - 1))
+    {
+        x[1][0] = m_data[dataIdx];
+        x[1][1] = m_data[dataIdx - m_points.x()];
+        run[1] = m_spacing.y();
+    }
+    else
+    {
+        x[1][0] = m_data[dataIdx + m_points.x()];
+        x[1][1] = m_data[dataIdx - m_points.x()];
+        run[1] = 2 * m_spacing.y();
+    }
+
+    if (k == 0)
+    {
+        x[2][0] = m_data[dataIdx + (m_points.x() * m_points.y())];
+        x[2][1] = m_data[dataIdx];
+        run[2] = m_spacing.z();
+    }
+    else if (k == (m_points.z() - 1))
+    {
+        x[2][0] = m_data[dataIdx];
+        x[2][1] = m_data[dataIdx - (m_points.x() * m_points.y())];
+        run[2] = m_spacing.z();
+    }
+    else
+    {
+        x[2][0] = m_data[dataIdx + (m_points.x() * m_points.y())];
+        x[2][1] = m_data[dataIdx - (m_points.x() * m_points.y())];
+        run[2] = 2 * m_spacing.z();
+    }
+
+    std::array<float, 3> ret;
+
+    ret[0] = (x[0][1] - x[0][0]) / run[0];
+    ret[1] = (x[1][1] - x[1][0]) / run[1];
+    ret[2] = (x[2][1] - x[2][0]) / run[2];
+
+    return ret;
+}
+
+std::array<std::array<float, 3>, 8>
+Cube::getGradCube(size_t i, size_t j, size_t k) const
+{
+    std::array<std::array<float, 3>, 8> grad;
+
+    grad[0] = computeGradient(i, j, k);
+    grad[1] = computeGradient(i + 1, j, k);
+    grad[2] = computeGradient(i + 1, j + 1, k);
+    grad[3] = computeGradient(i, j + 1, k);
+    grad[4] = computeGradient(i, j, k + 1);
+    grad[5] = computeGradient(i + 1, j, k + 1);
+    grad[6] = computeGradient(i + 1, j + 1, k + 1);
+    grad[7] = computeGradient(i, j + 1, k + 1);
+
+    return grad;
+}
+
+std::array<float, 8> Cube::getValsCube(size_t i, size_t j, size_t k) const
+{
+    std::array<float, 8> vals;
+
+    size_t idx = i + (j * m_points.x()) + (k * m_points.x() * m_points.y());
+
+    vals[0] = getData(i, j, k);
+    vals[1] = getData(i + 1, j, k);
+    vals[2] = getData(i + 1, j + 1, k);
+    vals[3] = getData(i, j + 1, k);
+    vals[4] = getData(i, j, k + 1);
+    vals[5] = getData(i + 1, j, k + 1);
+    vals[6] = getData(i + 1, j + 1, k + 1);
+    vals[7] = getData(i, j + 1, k + 1);
+
+    return vals;
+}
+
+
+inline float
+Cube::getData(size_t i, size_t j, size_t k) const
+{
+    return m_data[(k * m_points.x() * m_points.y()) + (j * m_points.x()) + i];
+}
+
+
+std::array<std::array<float, 3>, 8> Cube::getPosCube(size_t i, size_t j, size_t k) const
+{
+  std::array<std::array<float, 3>, 8> pos;
+
+  // std::cout<< m_spacing.x();
+  float xpos = m_min.x() + i * m_spacing.y();
+  float ypos = m_min.y() + j * m_spacing.x();
+  float zpos = m_min.z() + k * m_spacing.z();
+
+  pos[0][0] = xpos;
+  pos[0][1] = ypos;
+  pos[0][2] = zpos;
+
+  pos[1][0] = xpos + m_spacing.x();
+  pos[1][1] = ypos;
+  pos[1][2] = zpos;
+
+  pos[2][0] = xpos + m_spacing.x();
+  pos[2][1] = ypos + m_spacing.y();
+  pos[2][2] = zpos;
+
+  pos[3][0] = xpos;
+  pos[3][1] = ypos + m_spacing.y();
+  pos[3][2] = zpos;
+
+  pos[4][0] = xpos;
+  pos[4][1] = ypos;
+  pos[4][2] = zpos + m_spacing.z();
+
+  pos[5][0] = xpos + m_spacing.x();
+  pos[5][1] = ypos;
+  pos[5][2] = zpos + m_spacing.z();
+
+  pos[6][0] = xpos + m_spacing.x();
+  pos[6][1] = ypos + m_spacing.y();
+  pos[6][2] = zpos + m_spacing.z();
+
+  pos[7][0] = xpos;
+  pos[7][1] = ypos + m_spacing.y();
+  pos[7][2] = zpos + m_spacing.z();
+
+  // std::cout << "pos" << " " << pos[0][0] << std::endl;
+  return pos;
+}
+
 float Cube::value(const Vector3i& pos) const
 {
   unsigned int index =

avogadro/core/cube.h

@@ -13,6 +13,7 @@
 #include "vector.h"
 
 #include <vector>
+#include <array>  // Include this for std::array
 
 namespace Avogadro {
 namespace Core {
@@ -23,9 +24,7 @@ class Mutex;
 /**
  * @class Cube cube.h <avogadro/core/cube.h>
  * @brief Provide a data structure for regularly spaced 3D grids.
- * @author Marcus D. Hanwell
  */
-
 class AVOGADROCORE_EXPORT Cube
 {
 public:
@@ -118,7 +117,7 @@ class AVOGADROCORE_EXPORT Cube
   bool setLimits(const Molecule& mol, float spacing, float padding);
 
   /**
-   * @return Vector containing all the data in a one-dimensional array.
+   * @return Pointer to the vector containing all the data in a one-dimensional array.
    */
   std::vector<float>* data();
   const std::vector<float>* data() const;
@@ -133,6 +132,11 @@ class AVOGADROCORE_EXPORT Cube
    */
   bool addData(const std::vector<float>& values);
 
+  /**
+   * @return Const iterator to the beginning of a specific row in the cube data.
+   */
+  std::vector<float>::const_iterator getRowIter(size_t j, size_t k) const;
+
   /**
    * @return Index of the point closest to the position supplied.
    * @param pos Position to get closest index for.
@@ -201,7 +205,7 @@ class AVOGADROCORE_EXPORT Cube
    * @param value Value to fill the cube with.
    */
   void fill(float value);
-  
+
   /**
    * Sets all indices in a Z stripe of the cube to the specified value.
    * @param i x component of the position.
@@ -215,12 +219,12 @@ class AVOGADROCORE_EXPORT Cube
   );
 
   /**
-   * @return The minimum  value at any point in the Cube.
+   * @return The minimum value at any point in the Cube.
    */
   float minValue() const { return m_minValue; }
 
   /**
-   * @return The maximum  value at any point in the Cube.
+   * @return The maximum value at any point in the Cube.
    */
   float maxValue() const { return m_maxValue; }
 
@@ -235,11 +239,51 @@ class AVOGADROCORE_EXPORT Cube
    */
   Mutex* lock() const { return m_lock; }
 
+  // Add declarations for the functions defined in cube.cpp
+  /**
+   * Compute the gradient at a specific point in the cube.
+   * @param i x index
+   * @param j y index
+   * @param k z index
+   * @return Gradient vector at the specified point.
+   */
+  std::array<float, 3> computeGradient(size_t i, size_t j, size_t k) const;
+
+  /**
+   * Get the values of the eight corners of a cube cell.
+   * @param i x index
+   * @param j y index
+   * @param k z index
+   * @return Array of values at the eight corners.
+   */
+  std::array<float, 8> getValsCube(size_t i, size_t j, size_t k) const;
+
+
+  std::array<std::array<float, 3>, 8> getGradCube(size_t i, size_t j, size_t k) const;        
+  /**
+   * Get the data value at the specified indices.
+   * @param i x index
+   * @param j y index
+   * @param k z index
+   * @return Value at the specified indices.
+   */
+  float getData(size_t i, size_t j, size_t k) const;
+
+  /**
+   * Get the positions of the eight corners of a cube cell.
+   * @param i x index
+   * @param j y index
+   * @param k z index
+   * @return Array of positions at the eight corners.
+   */
+std::array<std::array<float, 3>, 8> getPosCube(size_t i, size_t j, size_t k) const;
+
 protected:
   std::vector<float> m_data;
   Vector3 m_min, m_max, m_spacing;
   Vector3i m_points;
   float m_minValue, m_maxValue;
+  // Removed conflicting 'data' member variable
   std::string m_name;
   Type m_cubeType;
   Mutex* m_lock;

avogadro/core/mesh.cpp

@@ -62,6 +62,20 @@ const Vector3f* Mesh::vertex(int n) const
   return &(m_vertices[n * 3]);
 }
 
+
+bool Mesh::setTriangles(const std::vector<std::array<size_t, 3> >& values)
+{
+  tris.clear();
+  tris = values;
+  return true;
+}
+
+const std::vector<std::array<size_t, 3> >& Mesh::triangles() const
+{
+  return tris;
+}
+
+
 bool Mesh::setVertices(const Core::Array<Vector3f>& values)
 {
   m_vertices.clear();

avogadro/core/mesh.h

@@ -122,6 +122,10 @@ class AVOGADROCORE_EXPORT Mesh
    */
   const Vector3f* vertex(int n) const;
 
+  bool setTriangles(const std::vector<std::array<size_t, 3> >& values);
+  const std::vector<std::array<size_t, 3> >& triangles() const;
+
+
   /**
    * Clear the vertices vector and assign new values.
    */
@@ -141,10 +145,10 @@ class AVOGADROCORE_EXPORT Mesh
   /**
    * @return The number of normals.
    */
-  unsigned int numNormals() const
-  {
-    return static_cast<unsigned int>(m_normals.size());
-  }
+  // unsigned int numNormals() const
+  // {
+  //   return static_cast<unsigned int>(m_normals.size());
+  // }
 
   /**
    * @return Pointer to the first normal of the specified triangle.
@@ -228,6 +232,7 @@ class AVOGADROCORE_EXPORT Mesh
   Core::Array<Vector3f> m_vertices;
   Core::Array<Vector3f> m_normals;
   Core::Array<Color3f> m_colors;
+  std::vector<std::array<size_t, 3> > tris;
   std::string m_name;
   bool m_stable;
   float m_isoValue;