diff --git a/avogadro/core/cube.h b/avogadro/core/cube.h index 359a43d5da..afafec6267 100644 --- a/avogadro/core/cube.h +++ b/avogadro/core/cube.h @@ -259,7 +259,8 @@ class AVOGADROCORE_EXPORT Cube std::array getValsCube(int i, int j, int k) const; - std::array, 8> getGradCube(int i, int j, int k) const; + std::array, 8> getGradCube(int i, int j, int k) const; + /** * Get the data value at the specified indices. * @param i x index diff --git a/avogadro/qtgui/meshgenerator.cpp b/avogadro/qtgui/meshgenerator.cpp index 796acd98eb..b6959b6fd9 100644 --- a/avogadro/qtgui/meshgenerator.cpp +++ b/avogadro/qtgui/meshgenerator.cpp @@ -68,15 +68,12 @@ bool MeshGenerator::initialize(const Cube* cube_, Mesh* mesh_, float iso, void MeshGenerator::FlyingEdgesAlgorithmPass1() { // Loop through z-dimension - int nx = m_dim.x(); - int ny = m_dim.y(); - int nz = m_dim.z(); - for(int k = 0; k != nz; ++k) { - for(int j = 0; j != ny; ++j) + for(int k = 0; k != m_dim.z(); ++k) { + for(int j = 0; j != m_dim.y(); ++j) { - auto curEdgeCases = edgeCases.begin() + (nx - 1) * (k * ny + j); + auto curEdgeCases = edgeCases.begin() + (m_dim.x() - 1) * (k * m_dim.y() + j); auto curPointValues = m_cube->getRowIter(j, k); @@ -96,10 +93,10 @@ void MeshGenerator::FlyingEdgesAlgorithmPass1() } } - for(int k = 0; k != nz; ++k){ - for(int j = 0; j != ny; ++j) + for(int k = 0; k != m_dim.z(); ++k){ + for(int j = 0; j != m_dim.y(); ++j) { - gridEdge& curGridEdge = gridEdges[k * ny + j]; + gridEdge& curGridEdge = gridEdges[k * m_dim.y() + j]; curGridEdge.xl = m_dim.x(); @@ -121,19 +118,16 @@ void MeshGenerator::FlyingEdgesAlgorithmPass1() void MeshGenerator::FlyingEdgesAlgorithmPass2() { - int nx = m_dim.x(); - int ny = m_dim.y(); - int nz = m_dim.z(); - for(int k = 0; k != nz - 1; ++k){ - for(int j = 0; j != ny - 1; ++j) + for(int k = 0; k != m_dim.z() - 1; ++k){ + for(int j = 0; j != m_dim.y() - 1; ++j) { int xl, xr; calcTrimValues(xl, xr, j, k); // xl, xr set in this function - gridEdge& ge0 = gridEdges[k * ny + j]; - gridEdge& ge1 = gridEdges[k* ny +j + 1]; - gridEdge& ge2 = gridEdges[(k+1) * ny + j]; - gridEdge& ge3 = gridEdges[(k+1) * ny + j + 1]; + gridEdge& ge0 = gridEdges[k * m_dim.y() + j]; + gridEdge& ge1 = gridEdges[k* m_dim.y() +j + 1]; + gridEdge& ge2 = gridEdges[(k+1) * m_dim.y() + j]; + gridEdge& ge3 = gridEdges[(k+1) * m_dim.y() + j + 1]; auto const& ec0 = edgeCases.begin() + (m_dim.x()-1) * (k * m_dim.y() + j); auto const& ec1 = edgeCases.begin() + (m_dim.x()-1) * (k * m_dim.y() + j + 1); @@ -141,19 +135,19 @@ void MeshGenerator::FlyingEdgesAlgorithmPass2() auto const& ec3 = edgeCases.begin() + (m_dim.x()-1) * ((k+1) * m_dim.y() + j + 1); // Count the number of triangles along this row of cubes - int& curTriCounter = *(triCounter.begin() + k * (ny - 1) + j); + int& curTriCounter = *(triCounter.begin() + k * (m_dim.y() - 1) + j); auto curCubeCaseIds = cubeCases.begin() + (m_dim.x() - 1) * (k * (m_dim.y() - 1) + j); - bool isYEnd = (j == ny - 2); - bool isZEnd = (k == nz - 2); + bool isYEnd = (j == m_dim.y() - 2); + bool isZEnd = (k == m_dim.z() - 2); for(int i = xl; i != xr; ++i) { - bool isXEnd = (i == nx - 2); + bool isXEnd = (i == m_dim.x() - 2); unsigned char caseId = calcCubeCase(ec0[i], ec1[i], ec2[i], ec3[i]); // todo cubeCase not decleared curCubeCaseIds[i] = caseId; @@ -210,13 +204,10 @@ void MeshGenerator::FlyingEdgesAlgorithmPass3() { int tmp; int triAccum = 0; - int nx = m_dim.x(); - int ny = m_dim.y(); - int nz = m_dim.z(); - for(int k = 0; k != nz -1; ++k) { - for(int j = 0; j != ny-1; ++j) + for(int k = 0; k != m_dim.z() -1; ++k) { + for(int j = 0; j != m_dim.y()-1; ++j) { - int& curTriCounter = triCounter[k*(ny-1)+j]; + int& curTriCounter = triCounter[k*(m_dim.y()-1)+j]; tmp = curTriCounter; curTriCounter = triAccum; @@ -224,10 +215,10 @@ void MeshGenerator::FlyingEdgesAlgorithmPass3() }} int pointAccum = 0; - for(int k = 0; k != nz; ++k) { - for(int j = 0; j != ny; ++j) + for(int k = 0; k != m_dim.z(); ++k) { + for(int j = 0; j != m_dim.y(); ++j) { - gridEdge& curGridEdge = gridEdges[(k * ny) + j]; + gridEdge& curGridEdge = gridEdges[(k * m_dim.y()) + j]; tmp = curGridEdge.xstart; curGridEdge.xstart = pointAccum; @@ -250,12 +241,9 @@ void MeshGenerator::FlyingEdgesAlgorithmPass3() void MeshGenerator::FlyingEdgesAlgorithmPass4() { - int nx = m_dim.x(); - int ny = m_dim.y(); - int nz = m_dim.z(); - for(int k = 0; k != nz -1; ++k) { - for(int j = 0; j != ny-1; ++j) + for(int k = 0; k != m_dim.z() -1; ++k) { + for(int j = 0; j != m_dim.y()-1; ++j) { // find adjusted trim values int xl, xr; @@ -264,7 +252,7 @@ void MeshGenerator::FlyingEdgesAlgorithmPass4() if(xl == xr) continue; - int triIdx = triCounter[(k*(ny-1)) + j]; + int triIdx = triCounter[(k*(m_dim.y()-1)) + j]; auto curCubeCaseIds = cubeCases.begin() + (m_dim.x()-1)*(k*(m_dim.y()-1) + j); gridEdge const& ge0 = gridEdges[k* m_dim.y() + j]; @@ -284,12 +272,12 @@ void MeshGenerator::FlyingEdgesAlgorithmPass4() int x3counter = 0; - bool isYEnd = (j == ny-2); - bool isZEnd = (k == nz-2); + bool isYEnd = (j == m_dim.y()-2); + bool isZEnd = (k == m_dim.z()-2); for(int i = xl; i != xr; ++i) { - bool isXEnd = (i == nx-2); + bool isXEnd = (i == m_dim.x()-2); unsigned char caseId = curCubeCaseIds[i]; @@ -589,12 +577,9 @@ unsigned char MeshGenerator::calcCubeCase( bool MeshGenerator::isCutEdge(int const& i, int const& j, int const& k) const { - int nx = m_dim.x(); - int ny = m_dim.y(); - int nz = m_dim.z(); // Assuming edgeCases are all set - int edgeCaseIdx = k * ((nx - 1) * ny) + (j * (nx - 1)) + i; + int edgeCaseIdx = k * ((m_dim.x() - 1) * m_dim.y()) + (j * (m_dim.x() - 1)) + i; unsigned char edgeCase = edgeCases[edgeCaseIdx]; if (edgeCase == 1 || edgeCase == 2) @@ -602,9 +587,9 @@ bool MeshGenerator::isCutEdge(int const& i, int const& j, int const& k) const return true; } - if (j != ny - 1) + if (j != m_dim.y() - 1) { - int edgeCaseIdxY = (k * (nx - 1) * ny) + ((j + 1) * (nx - 1)) + i; + int edgeCaseIdxY = (k * (m_dim.x() - 1) * m_dim.y()) + ((j + 1) * (m_dim.x() - 1)) + i; unsigned char edgeCaseY = edgeCases[edgeCaseIdxY]; // If the sum is odd, the edge along the y-axis is cut @@ -614,9 +599,9 @@ bool MeshGenerator::isCutEdge(int const& i, int const& j, int const& k) const } } - if (k != nz - 1) + if (k != m_dim.z() - 1) { - int edgeCaseIdxZ = ((k + 1) * (nx - 1) * ny) + (j * (nx - 1)) + i; + int edgeCaseIdxZ = ((k + 1) * (m_dim.x() - 1) * m_dim.y()) + (j * (m_dim.x() - 1)) + i; unsigned char edgeCaseZ = edgeCases[edgeCaseIdxZ]; // If the sum is odd, the edge along the z-axis is cut @@ -648,12 +633,11 @@ unsigned char MeshGenerator::calcCaseEdge(bool const& prevEdge, bool const& curr // 2000 void MeshGenerator::calcTrimValues(int& xl, int& xr, int const& j, int const& k) const { - int ny = m_dim.y(); - const gridEdge& ge0 = gridEdges[k * ny + j]; - const gridEdge& ge1 = gridEdges[k * ny + j + 1]; - const gridEdge& ge2 = gridEdges[(k + 1) * ny + j]; - const gridEdge& ge3 = gridEdges[(k + 1) * ny + j + 1]; + const gridEdge& ge0 = gridEdges[k * m_dim.y() + j]; + const gridEdge& ge1 = gridEdges[k * m_dim.y() + j + 1]; + const gridEdge& ge2 = gridEdges[(k + 1) * m_dim.y() + j]; + const gridEdge& ge3 = gridEdges[(k + 1) * m_dim.y() + j + 1]; xl = std::min({ge0.xl, ge1.xl, ge2.xl, ge3.xl}); xr = std::max({ge0.xr, ge1.xr, ge2.xr, ge3.xr}); diff --git a/avogadro/qtgui/meshgenerator.h b/avogadro/qtgui/meshgenerator.h index fb38199151..32f8157410 100644 --- a/avogadro/qtgui/meshgenerator.h +++ b/avogadro/qtgui/meshgenerator.h @@ -81,6 +81,10 @@ class AVOGADROQTGUI_EXPORT MeshGenerator : public QThread */ void run() override; + /** + * It holds the range and starting offsets of isosurface-intersected + * edges along the x, y, and z axes for each grid cell. + */ struct gridEdge { gridEdge() @@ -104,10 +108,29 @@ class AVOGADROQTGUI_EXPORT MeshGenerator : public QThread }; - + /** + * Pass 1 for flying edges. Pass1 detects and records + * where the isosurface intersects each row of grid edges + * along the x-axis. + */ void FlyingEdgesAlgorithmPass1(); + + /** + * Pass2 assigns case identifiers to each grid cell based on + * intersected edges and tallies the number of triangles needed + * for mesh construction. + */ void FlyingEdgesAlgorithmPass2(); + + /** + * Pass3 computes cumulative offsets for triangles + * and vertices and allocates memory for the mesh structures. + */ void FlyingEdgesAlgorithmPass3(); + + /** + * Calculates normals, triangles and vertices. + */ void FlyingEdgesAlgorithmPass4(); /** @@ -115,12 +138,16 @@ class AVOGADROQTGUI_EXPORT MeshGenerator : public QThread */ const Core::Cube* cube() const { return m_cube; } + /** + * Determines the x-range (xl to xr) where isosurface intersections + * occur, optimizing calculations within this range. + */ inline void calcTrimValues( int& xl, int& xr, int const& j, int const& k) const; - - void calculateTotalPoints(); - + /** + * Indicates which edges intersects the isosurface. + */ inline unsigned char calcCubeCase(unsigned char const& ec0, unsigned char const& ec1, unsigned char const& ec2, unsigned char const& ec3) const; @@ -146,6 +173,7 @@ class AVOGADROQTGUI_EXPORT MeshGenerator : public QThread int progressMaximum() { return m_progmax; } signals: + /** * The current value of the calculation's progress. */ @@ -171,22 +199,32 @@ class AVOGADROQTGUI_EXPORT MeshGenerator : public QThread unsigned long duplicate(const Vector3i& c, const Vector3f& pos); /** - * Perform a marching cubes step on a single cube. - */ - bool marchingCube(const Vector3i& pos); - + * isCutEdge checks whether the grid edge at position (i, j, k) is + * intersected by the isosurface based on edge case conditions. + * @return Boolean if it's intersected or not. + */ bool isCutEdge(int const& i, int const& j, int const& k) const; + /** + * It computes the 3D intersection point on a cube edge via interpolation. + */ inline std::array interpolateOnCube( std::array, 8> const& pts, std::array const& isovals, unsigned char const& edge) const; + /** + * It linearly interpolates between two 3D points, a and b, using + * the given weight to determine the intermediate position. + */ inline std::array interpolate( std::array const& a, std::array const& b, float const& weight) const; + /** + * calcCaseEdge determines an edge case code (0–3) based on two boolean edge comparisons. + */ inline unsigned char calcCaseEdge(bool const& prevEdge, bool const& currEdge) const; float m_iso; /** The value of the isosurface. */ @@ -198,26 +236,18 @@ class AVOGADROQTGUI_EXPORT MeshGenerator : public QThread Vector3f m_min; /** The minimum point in the cube. */ Vector3i m_dim; /** The dimensions of the cube. */ - std::array, 8> cube_t; - std::array scalarCube_t; - - Core::Array m_vertices, m_normals; - Core::Array m_indices; - std::vector gridEdges; - std::vector cubeCases; // size (nx-1)*(ny-1)*(nz-1) - - std::vector triCounter; // size of (ny-1)*(nz-1) - std::vector edgeCases; - Core::Array points; // - Core::Array normals; // The output - Core::Array tris; + Core::Array points, normals; + std::vector gridEdges; // size (m_dim.y() * m_dim.z()) + std::vector cubeCases; //size ((m_dim.x() - 1) * (m_dim.y() - 1) * (m_dim.z() - 1)) + std::vector triCounter; // size ((m_dim.y() - 1) * (m_dim.z() - 1)) + std::vector edgeCases; // size ((m_dim.x() - 1) * (m_dim.y()) * (m_dim.z())) + Core::Array tris; // triangles of a mesh int m_progmin; int m_progmax; /** - * These are the tables of constants for the marching cubes and tetrahedra - * algorithms. They are taken from the public domain source at - * http://local.wasp.uwa.edu.au/~pbourke/geometry/polygonise/ + * These are the lookup tables for flying edges. + * Reference : https://github.com/sandialabs/miniIsosurface/blob/master/flyingEdges/util/MarchingCubesTables.h */ static const unsigned char numTris[256]; static const bool isCut[256][12];