Use flying edges for mesh generation

avogadro/core/cube.cpp

@@ -37,7 +37,6 @@
   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 @@
     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 @@
   return true;
 }
 
-} // End Avogadro namespace
+} // End Avogadro namespace

avogadro/core/cube.h

@@ -13,6 +13,7 @@
 #include "vector.h"
 
 #include <vector>
+#include <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(int j, int 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,6 +239,51 @@ 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;
+
+  /**
+   * Get the positions 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 positions at the eight 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;