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,6 +24,17 @@ Cube::~Cube()
   m_lock = nullptr;
 }
 
+std::vector<float>::const_iterator Cube::getRowIter(size_t j, size_t k) const
+{
+    size_t startIdx = k * m_points.x() * m_points.y() + j * m_points.x();
+    if (startIdx >= m_data.size()) {
+        std::cout << "Error: startIdx out of bounds in getRowIter.\n";
+        // return m_data.cend();
+    }
+    return m_data.cbegin() + startIdx;
+}
+
+
 bool Cube::setLimits(const Vector3& min_, const Vector3& max_,
                      const Vector3i& points)
 {
@@ -194,6 +206,165 @@ 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();
+
+    if (i == 0)
+    {
+        x[0][0] = m_data[dataIdx + 1];
+        x[0][1] = m_data[dataIdx];
+        run[0] = m_spacing[0];
+    }
+    else if (i == (m_points.x() - 1))
+    {
+        x[0][0] = m_data[dataIdx];
+        x[0][1] = m_data[dataIdx - 1];
+        run[0] = m_spacing[0];
+    }
+    else
+    {
+        x[0][0] = m_data[dataIdx + 1];
+        x[0][1] = m_data[dataIdx - 1];
+        run[0] = 2 * m_spacing[0];
+    }
+
+    if (j == 0)
+    {
+        x[1][0] = m_data[dataIdx + m_points.x()];
+        x[1][1] = m_data[dataIdx];
+        run[1] = m_spacing[1];
+    }
+    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[1];
+    }
+    else
+    {
+        x[1][0] = m_data[dataIdx + m_points.x()];
+        x[1][1] = m_data[dataIdx - m_points.x()];
+        run[1] = 2 * m_spacing[1];
+    }
+
+    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[2];
+    }
+    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[2];
+    }
+    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[2];
+    }
+
+    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;
+
+  float xpos = m_min.x() + i * m_spacing.x();
+  float ypos = m_min.y() + j * m_spacing.y();
+  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();
+
+  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;