Use flying edges for mesh generation

avogadro/core/cube.cpp

@@ -37,7 +37,6 @@ bool Cube::setLimits(const Vector3& min_, const Vector3& max_,
   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,8 +193,159 @@ float Cube::value(int i, int j, int k) const
     return 0.0;
 }
 
-float Cube::value(const Vector3i& pos) const
+std::array<float, 3> Cube::computeGradient(int i, int j, int k) const
+{
+    int nx = m_points.x();
+    int ny = m_points.y();
+    int nz = m_points.z();
+    int dataIdx = (i * ny * nz) + (j * nz) + k;
+
+    std::array<std::array<float, 2>, 3> x;
+    std::array<float, 3> run;
+
+    // X-direction
+    if (i == 0) {
+        x[0][0] = m_data[dataIdx + ny * nz];
+        x[0][1] = m_data[dataIdx];
+        run[0] = m_spacing.x();
+    } else if (i == (nx - 1)) {
+        x[0][0] = m_data[dataIdx];
+        x[0][1] = m_data[dataIdx - ny * nz];
+        run[0] = m_spacing.x();
+    } else {
+        x[0][0] = m_data[dataIdx + ny * nz];
+        x[0][1] = m_data[dataIdx - ny * nz];
+        run[0] = 2 * m_spacing.x();
+    }
+
+    // Y-direction
+    if (j == 0) {
+        x[1][0] = m_data[dataIdx + nz];
+        x[1][1] = m_data[dataIdx];
+        run[1] = m_spacing.y();
+    } else if (j == (ny - 1)) {
+        x[1][0] = m_data[dataIdx];
+        x[1][1] = m_data[dataIdx - nz];
+        run[1] = m_spacing.y();
+    } else {
+        x[1][0] = m_data[dataIdx + nz];
+        x[1][1] = m_data[dataIdx - nz];
+        run[1] = 2 * m_spacing.y();
+    }
+
+    // Z-direction
+    if (k == 0) {
+        x[2][0] = m_data[dataIdx + 1];
+        x[2][1] = m_data[dataIdx];         
+        run[2] = m_spacing.z();
+    } else if (k == (nz - 1)) {
+        x[2][0] = m_data[dataIdx];        
+        x[2][1] = m_data[dataIdx - 1];     
+        run[2] = m_spacing.z();
+    } else {
+        x[2][0] = m_data[dataIdx + 1];    
+        x[2][1] = m_data[dataIdx - 1];     
+        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(int i, int j, int 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(int i, int j, int k) const
+{
+    std::array<float, 8> vals;
+
+    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;
+}
+
+
+float Cube::getData(int i, int j, int k) const
+{
+    int nx = m_points.x();
+    int ny = m_points.y();
+    int nz = m_points.z();
+    return m_data[(i * ny * nz) + (j * nz) + k];
+}
+
+
+std::array<std::array<float, 3>, 8> Cube::getPosCube(int i, int j, int 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 =
     pos.x() * m_points.y() * m_points.z() + pos.y() * m_points.z() + pos.z();
   if (index < m_data.size())
@@ -294,4 +444,4 @@ bool Cube::fillStripe(
   return true;
 }
 
-} // End Avogadro namespace
+} // End Avogadro namespace

avogadro/core/cube.h

@@ -11,6 +11,7 @@
 #include "vector.h"
 
 #include <vector>
+#include <array>
 
 namespace Avogadro::Core {
 
@@ -21,6 +22,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
+ * @author Perminder Singh
  */
 
 class AVOGADROCORE_EXPORT Cube
@@ -115,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;
@@ -130,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(int j, int k) const;
+
   /**
    * @return Index of the point closest to the position supplied.
    * @param pos Position to get closest index for.
@@ -198,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.
@@ -212,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; }
 
@@ -232,6 +239,56 @@ class AVOGADROCORE_EXPORT Cube
    */
   Mutex* lock() const { return m_lock; }
 
+  /**
+   * 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(int i, int j, int k) const;
+
+  /**
+   * Get the values of the eight corners of a cube defined by the indices (i, j, k).
+   * @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(int i, int j, int k) const;
+
+  /**
+  * Get the gradients at the eight corners of the cube defined by the indices (i, j, k).
+  * @param i x index
+  * @param j y index
+  * @param k z index
+  * @return Array of gradients at the eight corners. Each gradient is a 3D vector.
+  */
+  std::array<std::array<float, 3>, 8> getGradCube(int i, int j, int 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(int i, int j, int k) const;
+
+/**
+ * Retrieves the positions of the eight corners of a cube at grid indices (i, j, k).
+ *
+ * The indices (i, j, k) are converted to real-space positions (xpos, ypos, zpos), mapping grid indices to physical coordinates.
+ * The method returns a cube that spans one step in each of the x, y, and z directions, with step sizes defined by `m_spacing`.
+ *
+ * @param i X-index.
+ * @param j Y-index.
+ * @param k Z-index.
+ * @return A `std::array` of eight `(x, y, z)` coordinates representing the cube's corners.
+ */
+
+std::array<std::array<float, 3>, 8> getPosCube(int i, int j, int k) const;
+
 protected:
   std::vector<float> m_data;
   Vector3 m_min, m_max, m_spacing;