diff --git a/avogadro/core/cube.cpp b/avogadro/core/cube.cpp index b9f8dabf11..868cdaea35 100644 --- a/avogadro/core/cube.cpp +++ b/avogadro/core/cube.cpp @@ -24,27 +24,23 @@ Cube::~Cube() m_lock = nullptr; } -std::vector::const_iterator Cube::getRowIter(int j, int k) const -{ - return m_data.cbegin() + m_points.x() * (k * m_points.y() + j); -} +// std::vector::const_iterator Cube::getRowIter(int j, int k) const +// { +// return m_data.cbegin() + m_points.z() * k * m_points.y() + m_points.z() * j; +// } bool Cube::setLimits(const Vector3& min_, const Vector3& max_, const Vector3i& points) { - // 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 + // We can calculate all necessary properties and initialise our data Vector3 delta = max_ - min_; - m_spacing = Vector3(delta.x() / (equalPoints.x() - 1), - delta.y() / (equalPoints.y() - 1), - delta.z() / (equalPoints.z() - 1)); + m_spacing = + Vector3(delta.x() / (points.x() - 1), delta.y() / (points.y() - 1), + delta.z() / (points.z() - 1)); m_min = min_; m_max = max_; - m_points = equalPoints; + m_points = points; m_data.resize(m_points.x() * m_points.y() * m_points.z()); return true; } @@ -205,69 +201,58 @@ float Cube::value(int i, int j, int k) const else return 0.0; } - std::array 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, 3> x; std::array run; - int 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]; + // X-direction + if (i == 0) { + x[0][0] = m_data[dataIdx + ny * nz]; // (i+1, j, k) + x[0][1] = m_data[dataIdx]; // (i, j, k) 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]; + } else if (i == (nx - 1)) { + x[0][0] = m_data[dataIdx]; // (i, j, k) + x[0][1] = m_data[dataIdx - ny * nz]; // (i-1, j, k) run[0] = m_spacing.x(); - } - else - { - x[0][0] = m_data[dataIdx + 1]; - x[0][1] = m_data[dataIdx - 1]; + } else { + x[0][0] = m_data[dataIdx + ny * nz]; // (i+1, j, k) + x[0][1] = m_data[dataIdx - ny * nz]; // (i-1, j, k) run[0] = 2 * m_spacing.x(); } - if (j == 0) - { - x[1][0] = m_data[dataIdx + m_points.x()]; - x[1][1] = m_data[dataIdx]; + // Y-direction + if (j == 0) { + x[1][0] = m_data[dataIdx + nz]; // (i, j+1, k) + x[1][1] = m_data[dataIdx]; // (i, j, k) 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()]; + } else if (j == (ny - 1)) { + x[1][0] = m_data[dataIdx]; // (i, j, k) + x[1][1] = m_data[dataIdx - nz]; // (i, j-1, k) run[1] = m_spacing.y(); - } - else - { - x[1][0] = m_data[dataIdx + m_points.x()]; - x[1][1] = m_data[dataIdx - m_points.x()]; + } else { + x[1][0] = m_data[dataIdx + nz]; // (i, j+1, k) + x[1][1] = m_data[dataIdx - nz]; // (i, j-1, k) 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]; + // Z-direction + if (k == 0) { + x[2][0] = m_data[dataIdx + 1]; // (i, j, k+1) + x[2][1] = m_data[dataIdx]; // (i, j, k) 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())]; + } else if (k == (nz - 1)) { + x[2][0] = m_data[dataIdx]; // (i, j, k) + x[2][1] = m_data[dataIdx - 1]; // (i, j, k-1) 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())]; + } else { + x[2][0] = m_data[dataIdx + 1]; // (i, j, k+1) + x[2][1] = m_data[dataIdx - 1]; // (i, j, k-1) run[2] = 2 * m_spacing.z(); } @@ -314,10 +299,12 @@ std::array Cube::getValsCube(int i, int j, int k) const } -inline float -Cube::getData(int i, int j, int k) const +inline float Cube::getData(int i, int j, int k) const { - return m_data[(k * m_points.x() * m_points.y()) + (j * m_points.x()) + i]; + int nx = m_points.x(); + int ny = m_points.y(); + int nz = m_points.z(); + return m_data[(i * ny * nz) + (j * nz) + k]; } @@ -366,7 +353,7 @@ std::array, 8> Cube::getPosCube(int i, int j, int k) const } 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()) @@ -465,4 +452,4 @@ bool Cube::fillStripe( return true; } -} // End Avogadro namespace +} // End Avogadro namespace \ No newline at end of file diff --git a/avogadro/qtgui/meshgenerator.cpp b/avogadro/qtgui/meshgenerator.cpp index b6959b6fd9..37813b1bd8 100644 --- a/avogadro/qtgui/meshgenerator.cpp +++ b/avogadro/qtgui/meshgenerator.cpp @@ -69,29 +69,19 @@ void MeshGenerator::FlyingEdgesAlgorithmPass1() { // Loop through z-dimension - for(int k = 0; k != m_dim.z(); ++k) { - for(int j = 0; j != m_dim.y(); ++j) - { - - auto curEdgeCases = edgeCases.begin() + (m_dim.x() - 1) * (k * m_dim.y() + j); - - auto curPointValues = m_cube->getRowIter(j, k); - - // Initialize the isGE array to track isosurface crossings - std::array isGE; - isGE[0] = (curPointValues[0] >= m_iso); // first element initialization - - // loop through x-dimension - for(int i = 1; i != m_dim.x(); ++i) - { - // update isGE for the current point + for (int k = 0; k != m_dim.z(); ++k) { + for (int j = 0; j != m_dim.y(); ++j) { - isGE[i % 2] = (curPointValues[i] >= m_iso); + auto curEdgeCases = edgeCases.begin() + (m_dim.x() - 1) * (k * m_dim.y() + j); + std::array isGE; + isGE[0] = (m_cube->getData(0, j, k) >= m_iso); - curEdgeCases[i-1] = calcCaseEdge(isGE[(i+1)%2], isGE[i%2]); + for (int i = 1; i != m_dim.x(); ++i) { + isGE[i % 2] = (m_cube->getData(i, j, k) >= m_iso); + curEdgeCases[i - 1] = calcCaseEdge(isGE[(i + 1) % 2], isGE[i % 2]); + } + } } - } -} for(int k = 0; k != m_dim.z(); ++k){ for(int j = 0; j != m_dim.y(); ++j) @@ -471,6 +461,7 @@ void MeshGenerator::FlyingEdgesAlgorithmPass4() const char* caseTri = caseTriangles[caseId]; // size 16 for(int idx = 0; caseTri[idx] != -1; idx += 3) { + tris[triIdx][0] = globalIdxs[caseTri[idx]]; tris[triIdx][1] = globalIdxs[caseTri[idx+1]]; tris[triIdx][2] = globalIdxs[caseTri[idx+2]]; @@ -488,7 +479,6 @@ void MeshGenerator::run() qDebug() << "No mesh or cube set - nothing to find isosurface of…"; return; } - m_mesh->setStable(false); m_mesh->clear();