Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Use flying edges for mesh generation #1741

Merged
merged 20 commits into from
Nov 16, 2024

Conversation

perminder-17
Copy link
Contributor

Resolves issue : #901
Developer Certificate of Origin
Version 1.1

Copyright (C) 2004, 2006 The Linux Foundation and its contributors.
1 Letterman Drive
Suite D4700
San Francisco, CA, 94129

Everyone is permitted to copy and distribute verbatim copies of this
license document, but changing it is not allowed.

Developer's Certificate of Origin 1.1

By making a contribution to this project, I certify that:

(a) The contribution was created in whole or in part by me and I
have the right to submit it under the open source license
indicated in the file; or

(b) The contribution is based upon previous work that, to the best
of my knowledge, is covered under an appropriate open source
license and I have the right under that license to submit that
work with modifications, whether created in whole or in part
by me, under the same open source license (unless I am
permitted to submit under a different license), as indicated
in the file; or

(c) The contribution was provided directly to me by some other
person who certified (a), (b) or (c) and I have not modified
it.

(d) I understand and agree that this project and the contribution
are public and that a record of the contribution (including all
personal information I submit with it, including my sign-off) is
maintained indefinitely and may be redistributed consistent with
this project or the open source license(s) involved.

@ghutchis ghutchis linked an issue Oct 17, 2024 that may be closed by this pull request
avogadro/core/cube.cpp Fixed Show fixed Hide fixed
avogadro/core/cube.cpp Fixed Show fixed Hide fixed
avogadro/qtgui/meshgenerator.cpp Fixed Show fixed Hide fixed
avogadro/qtgui/meshgenerator.cpp Fixed Show fixed Hide fixed
avogadro/qtgui/meshgenerator.cpp Fixed Show fixed Hide fixed
avogadro/qtgui/meshgenerator.cpp Fixed Show fixed Hide fixed
avogadro/qtgui/meshgenerator.cpp Fixed Show fixed Hide fixed
avogadro/qtgui/meshgenerator.cpp Fixed Show fixed Hide fixed
avogadro/qtgui/meshgenerator.cpp Fixed Show fixed Hide fixed
avogadro/qtgui/meshgenerator.cpp Fixed Show fixed Hide fixed
Signed-off-by: Perminder <[email protected]>
avogadro/core/cube.cpp Fixed Show fixed Hide fixed
avogadro/core/cube.cpp Fixed Show fixed Hide fixed
avogadro/qtgui/meshgenerator.cpp Fixed Show fixed Hide fixed
avogadro/qtgui/meshgenerator.cpp Fixed Show fixed Hide fixed
avogadro/qtgui/meshgenerator.cpp Fixed Show fixed Hide fixed
avogadro/qtgui/meshgenerator.cpp Fixed Show fixed Hide fixed
avogadro/qtgui/meshgenerator.cpp Fixed Show fixed Hide fixed
avogadro/qtgui/meshgenerator.cpp Fixed Show fixed Hide fixed
avogadro/qtgui/meshgenerator.cpp Fixed Show fixed Hide fixed
avogadro/qtgui/meshgenerator.cpp Fixed Show fixed Hide fixed
Copy link
Contributor

github-actions bot commented Nov 2, 2024

Here are the build results
macOS.dmg
Artifacts will only be retained for 90 days.

Signed-off-by: Perminder <[email protected]>
Comment on lines 510 to 529
// Add triangles
const char* caseTri = caseTriangles[caseId]; // size 16
for(int idx = 0; caseTri[idx] != -1; idx += 3)
{


size_t globalIdx = globalIdxs[caseTri[idx]];
size_t globalIdx1 = globalIdxs[caseTri[idx + 1]];
size_t globalIdx2 = globalIdxs[caseTri[idx + 2]];



m_vertices.push_back(points[globalIdx]);
m_vertices.push_back(points[globalIdx1]);
m_vertices.push_back(points[globalIdx2]);

m_normals.push_back(normals[globalIdx]);
m_normals.push_back(normals[globalIdx1]);
m_normals.push_back(normals[globalIdx2]);

Copy link
Contributor Author

@perminder-17 perminder-17 Nov 3, 2024

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Hi @cryos and @ghutchis , does this part looks okay to you guys. I think I have done all things right but my output comes out to be something like. Any suggestions? Do you think I am pushing for the vertices and normals in a right way?

image

@perminder-17
Copy link
Contributor Author

perminder-17 commented Nov 3, 2024

I can do cleanups in the code once I start getting the expected output.

Also tagging @matterhorn103 and @aerkiaga for their inputs.

Let me know if you guys have any feedbacks/suggestions?

avogadro/core/cube.cpp Fixed Show fixed Hide fixed
avogadro/qtgui/meshgenerator.cpp Fixed Show fixed Hide fixed
avogadro/qtgui/meshgenerator.cpp Fixed Show fixed Hide fixed
avogadro/qtgui/meshgenerator.cpp Fixed Show fixed Hide fixed
avogadro/qtgui/meshgenerator.cpp Fixed Show fixed Hide fixed
Signed-off-by: Perminder <[email protected]>
avogadro/core/cube.cpp Fixed Show fixed Hide fixed
avogadro/core/cube.cpp Fixed Show fixed Hide fixed
avogadro/core/cube.cpp Fixed Show fixed Hide fixed
avogadro/core/mesh.h Fixed Show resolved Hide resolved
avogadro/qtgui/meshgenerator.cpp Fixed Show fixed Hide fixed
avogadro/qtgui/meshgenerator.cpp Fixed Show fixed Hide fixed
avogadro/qtgui/meshgenerator.cpp Fixed Show fixed Hide fixed
avogadro/qtgui/meshgenerator.cpp Fixed Show fixed Hide fixed
avogadro/qtgui/meshgenerator.cpp Fixed Show fixed Hide fixed
avogadro/qtplugins/meshes/meshes.cpp Fixed Show fixed Hide fixed
@matterhorn103
Copy link
Contributor

@perminder-17 I'm afraid rendering is something I know absolutely nothing about

@perminder-17
Copy link
Contributor Author

@perminder-17 I'm afraid rendering is something I know absolutely nothing about

Aah... Thanks for your reply. No problem, I figured that out. It was like, we need to pass the triangles as well.

@ghutchis
Copy link
Member

ghutchis commented Nov 6, 2024

Great, I'll start a review tomorrow

avogadro/qtgui/meshgenerator.cpp Fixed Show fixed Hide fixed
avogadro/qtgui/meshgenerator.cpp Fixed Show fixed Hide fixed
unsigned char edgeCaseZ = edgeCases[edgeCaseIdxZ];

// If the sum is odd, the edge along the z-axis is cut
if ((edgeCase + edgeCaseZ) % 2 == 1)

Check warning

Code scanning / CodeQL

Bad check for oddness Warning

Possibly invalid test for oddness. This will fail for negative numbers.
Copy link
Contributor

github-actions bot commented Nov 7, 2024

Here are the build results
macOS.dmg
Artifacts will only be retained for 90 days.

@ghutchis ghutchis changed the title WIP : Consider switching to flying edges from parallel marching cubes Use flying edges for mesh generation Nov 8, 2024
@avo-bot
Copy link

avo-bot commented Nov 8, 2024

This pull request has been mentioned on Avogadro Discussion. There might be relevant details there:

https://discuss.avogadro.cc/t/november-2024-live-updates/6446/1

avogadro/core/cube.cpp Fixed Show fixed Hide fixed
avogadro/core/cube.cpp Fixed Show fixed Hide fixed
@perminder-17 perminder-17 marked this pull request as ready for review November 10, 2024 01:54
@perminder-17
Copy link
Contributor Author

perminder-17 commented Nov 10, 2024

Hi everyone, Here's the results after the implementation of flying-edges. requesting @ghutchis for the review, still needs some cleanups in the code, removal of unnecessary codeblocks and variables.

Marching cube takes 1.5 minutes to render the mesh

Screenshot from 2024-11-07 02-04-36

Flying edges takes only 3 seconds to render it. (could vary with different systems)

Screenshot from 2024-11-07 02-06-40

Also it produces a little smoother mesh (left->flying edges, right->marching cubes)

Screenshot from 2024-11-07 18-56-29

After Laplacian smoothing we have:

Screenshot from 2024-11-10 07-17-54
Screenshot from 2024-11-10 07-18-29

@ghutchis
Copy link
Member

I'm happy to help do some performance tests - I'm not seeing the "Time taken" code at the moment in your pull request.

I have a range of molecule sizes and resolutions to give this a good benchmark.

Doing a rough stopwatch count, I'm seeing ~2-2.5x performance improvement, but that's including the whole thing, not just mesh generation.

I can say that with caffeine as an example, the flying edges version appears instantly on my MacBook even at "Very High" resolution. With a larger molecule, it's going from ~10s to display the VdW mesh to ~4.5s.

If you want to send me the patch you used for timing, I'd be happy to perform some real benchmarks tonight in release mode.

@perminder-17
Copy link
Contributor Author

perminder-17 commented Nov 12, 2024

All I see is debug mode in what you posted, if that is the case my point stands - benchmarks in debug mode are not very useful and we know the STL etc is very slow in debug. Comparing release numbers (which I also expect to be good) would be more informative.

Thanks @cryos for your insights. Really helps a lot and also thanks @ghutchis for the help in performance testing.

If I could say what code I used for checking the time, I used the chrono library https://en.cppreference.com/w/cpp/chrono/high_resolution_clock/now

If you include header #include<chrono> and replace your run() function with this, you can see the time taken.

#include <chrono> // at the top (header file)



// Rest code remains the same


// Replace the run function to test
 
void MeshGenerator::run() // you can directly copy the run() function and paste it with the original ones
{

  if (!m_cube || !m_mesh) {
    qDebug() << "No mesh or cube set - nothing to find isosurface of…";
    return;
  }
  m_mesh->setStable(false);
  m_mesh->clear();

  auto start = std::chrono::high_resolution_clock::now();
  FlyingEdgesAlgorithmPass1();
  FlyingEdgesAlgorithmPass2();
  FlyingEdgesAlgorithmPass3();
  FlyingEdgesAlgorithmPass4();

  auto end = std::chrono::high_resolution_clock::now();
  std::chrono::duration<double> elapsed = end - start;
  qDebug() << "Flying Edges Time: " << elapsed.count() << " s";    

  m_mesh->setVertices(points);
  m_mesh->setNormals(normals);
  m_mesh->setTriangles(tris);
  m_mesh->smooth(m_passes);
  m_mesh->setStable(true);

  points.resize(0);
  normals.resize(0);
  tris.resize(0);

}

Screenshot from 2024-11-13 02-37-19

also, if you want to test and compare the time with marching cubes.

here how you can check the time taken for marching cubes as well.

OLD CODE

#include <chrono> // similarly include this

// just a little changes in the run function 

void MeshGenerator::run()
{

// rest code remains the same.

 auto start = std::chrono::high_resolution_clock::now(); // add this

  for (int i = 0; i < m_dim.x() - 1; ++i) {
    for (int j = 0; j < m_dim.y() - 1; ++j) {
      for (int k = 0; k < m_dim.z() - 1; ++k) {
        marchingCube(Vector3i(i, j, k));
      }
    }
    if (m_vertices.capacity() < m_vertices.size() + m_dim.y() * m_dim.x() * 3) {
      m_vertices.reserve(m_vertices.capacity() * 2);
      m_normals.reserve(m_normals.capacity() * 2);
    }
    emit progressValueChanged(i);
  }

// add this
  auto end = std::chrono::high_resolution_clock::now();
  std::chrono::duration<double> elapsed = end - start;
  qDebug() << "Marching Cube Timing: " << elapsed.count() << " s";    

// rest code remains the same
}

@perminder-17
Copy link
Contributor Author

With a larger molecule, it's going from ~10s to display the VdW mesh to ~4.5s.

Also, after the merge of this, I can handle the merge conflicts (if any) in this PR #1692. If we add a progressive refinement, it won't take even 4.5 seconds to render and will generate the mesh in no time.

Do let me know if everything looks fine in the results, I could do many cleanups in the code (added too many debugging peices, unnecessary functions, variables haha). Still needs some formatting and cleanups.

@perminder-17
Copy link
Contributor Author

Btw @ghutchis I just noticed one bug (when you click on the calculate button twice for the same mesh generation it takes more time). It's because we were not clearing the memory. In the above commit I have fixed that as well. Sorry for the oversight.

@ghutchis
Copy link
Member

I'm seeing ~2x - 2.2x speedup for mesh generation on a range of 7 molecules with 0.05Å resolution. It depends a lot on the size of the molecule and the shape (e.g., smaller molecules and more spherical / cubic see the greatest speedup).

I'll do a bit of profiling, but of course this is also the single-core implementation (vs. parallel)

@perminder-17
Copy link
Contributor Author

perminder-17 commented Nov 13, 2024

IMG-20241113-WA0002

Correct! I looked at the results mentioned on the paper. They also have 2-2.2x speedups.

I'll do some cleanups in the code tomorrow.

@ghutchis
Copy link
Member

I also checked the number of vertices for marching cubes vs. flying edges. It's getting exactly 6.0 times fewer vertices, so the mesh is much smaller at the same resolution (and thus faster to render, less memory, etc.)

Even on a really big DNA chunk, I'm seeing ~0.1s for mesh creation with the "automatic" resolution. (And yes, the progressive refinement could be a help here too.)

Signed-off-by: Perminder <[email protected]>

float Cube::getData(int i, int j, int k) const
{
int nx = m_points.x();

Check notice

Code scanning / CodeQL

Unused local variable Note

Variable nx is not used.
Copy link
Contributor

Here are the build results
Avogadro2.AppImage
Ubuntu-Latest.tar.gz
macOS.dmg
Artifacts will only be retained for 90 days.

Copy link
Member

@ghutchis ghutchis left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Mostly some small things with the comments and documentation consistency.

std::array<float, 3> computeGradient(int i, int j, int k) const;

/**
* Get the values of the eight corners of a cube cell.
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm not quite sure I understand from the docs here. I've got some index (i, j, k) into the cube. What's the "cell" of the cube? Is it just going to get me the corners of the entire cube? If so, why do I need i, j, k parameters?

avogadro/core/cube.h Show resolved Hide resolved
float getData(int i, int j, int k) const;

/**
* Get the positions of the eight corners of a cube cell.
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Again, what's the "cell" here? Are these the min, max of the corners of the cube? If so, why do you need the i, j, k?

avogadro/core/mesh.h Fixed Show resolved Hide resolved
@@ -228,6 +232,7 @@ class AVOGADROCORE_EXPORT Mesh
Core::Array<Vector3f> m_vertices;
Core::Array<Vector3f> m_normals;
Core::Array<Color3f> m_colors;
Core::Array<Vector3f> tris;
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Please use m_triangles to be consistent with the other variables.

avogadro/qtgui/meshgenerator.h Show resolved Hide resolved
@@ -107,34 +173,44 @@ class AVOGADROQTGUI_EXPORT MeshGenerator : public QThread
int progressMaximum() { return m_progmax; }

signals:

/**
* The current value of the calculation's progress.
*/
void progressValueChanged(int);

protected:
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Can we keep this? Sometimes it's helpful to get normals for arbitrary points.

*/
Vector3f normal(const Vector3f& pos);

unsigned long duplicate(const Vector3i& c, const Vector3f& pos);
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Needs documentation.

std::vector<unsigned char> cubeCases; //size ((m_dim.x() - 1) * (m_dim.y() - 1) * (m_dim.z() - 1))
std::vector<int> triCounter; // size ((m_dim.y() - 1) * (m_dim.z() - 1))
std::vector<unsigned char> edgeCases; // size ((m_dim.x() - 1) * (m_dim.y()) * (m_dim.z()))
Core::Array<Vector3f> tris; // triangles of a mesh
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Please remove any un-used variables (from Marching Cubes)

Please also use m_triangles (or similar) to keep consistent with member variables through the rest of the code.

static const unsigned char numTris[256];
static const bool isCut[256][12];
static const char caseTriangles[256][16];
static const unsigned char edgeVertices[12][2];
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The static analysis says these class members are never used?
If they're used, please keep to m_edgeVertices etc. for consistency.

Signed-off-by: Perminder <[email protected]>
Signed-off-by: Perminder <[email protected]>
Copy link
Contributor

Here are the build results
Avogadro2.AppImage
Ubuntu-Latest.tar.gz
macOS.dmg
Artifacts will only be retained for 90 days.

1 similar comment
Copy link
Contributor

Here are the build results
Avogadro2.AppImage
Ubuntu-Latest.tar.gz
macOS.dmg
Artifacts will only be retained for 90 days.


std::array<std::array<float, 3>, 8> pos;

float xpos = m_min.x() + (i * m_spacing.x());
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm sorry - I still don't understand what this method is supposed to do.

xpos, ypos, zpos are the real-space versions of i, j, k indices
.. so the method gives me a cube that's one step in each of the x, y, z directions?

If so, shouldn't that be the doc string?

@@ -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
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@author Marcus D. Hanwell
@author Perminder Singh

@@ -138,14 +142,6 @@ class AVOGADROCORE_EXPORT Mesh
*/
const Core::Array<Vector3f>& normals() const;

/**
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I still don't understand why you're removing numNormals() which seems like a useful method.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

But i was unable to find where we were using numnormals() in our actual code?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We're not, currently. But since this is a core class, it's also API for future development.

Copy link
Member

@ghutchis ghutchis left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Mostly resolved - please note the other comments, thanks.

Signed-off-by: Perminder <[email protected]>
Copy link
Contributor

Here are the build results
Ubuntu-Latest.tar.gz
Avogadro2.AppImage
Artifacts will only be retained for 90 days.

1 similar comment
Copy link
Contributor

Here are the build results
Ubuntu-Latest.tar.gz
Avogadro2.AppImage
Artifacts will only be retained for 90 days.

Copy link
Member

@ghutchis ghutchis left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Sorry, I missed these last time. I think this fixes the compile on Windows.

ge2.xstart += isCutCase[4];
ge2.ystart += isCutCase[7];
}
if(isXEnd and isYEnd)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

&& not and - this is the build failure on Windows

{
ge1.zstart += isCutCase[11];
}
if(isXEnd and isZEnd)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Same - use &&

{
ge2.ystart += isCutCase[5];
}
if(isYEnd and isZEnd)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

&&

if(isCutCase[11])
{
int idx = ge1.zstart + z1counter;
if(isXEnd and isYEnd)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

&& not and

if(isCutCase[5])
{
int idx = ge2.ystart + y2counter;
if(isXEnd and isZEnd)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

&& not and

if(isCutCase[6])
{
int idx = ge3.xstart + x3counter;
if(isYEnd and isZEnd)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

&& not and

Signed-off-by: Perminder Singh <[email protected]>
@perminder-17
Copy link
Contributor Author

Done @ghutchis

Copy link
Contributor

Here are the build results
Ubuntu-Latest.tar.gz
Avogadro2.AppImage
macOS.dmg
Win64.exe
Artifacts will only be retained for 90 days.

@ghutchis ghutchis merged commit 141745d into OpenChemistry:master Nov 16, 2024
22 of 23 checks passed
Copy link
Contributor

Here are the build results
Avogadro2.AppImage
Ubuntu-Latest.tar.gz
Win64.exe
Artifacts will only be retained for 90 days.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

Consider switching to flying edges from parallel marching cubes
5 participants