Hi, I am trying to implement flying-edges. Can somebody help me? I am working on a software which previously using marching cubes and I am shifting it flying-edges. (3d). Is this the right place to ask a question here?
1 Like
Look in the VTK Examples, there are arouund 10 C++ examples (with corresponding Python ones) using vtkFlyingEdges3D, may of these examples also use vtkMarchingCubes as an option. E.g. (MedicalDemo1)[https://examples.vtk.org/site/Cxx/Medical/MedicalDemo1/], also here is a list of classes demonstraten in the examples: VTK Classes Used
Hopefully you’ve read the paper:
https://www.researchgate.net/publication/282975362_Flying_Edges_A_High-Performance_Scalable_Isocontouring_Algorithm.
thankyou so much WILL. I have read that and I feel I have implemented the same way as they have described. But My output comes out to be something like this:
It should basically cover the whole molecule.
Also, I have a bit confusion, can you confirm if I am pushing m_vertices and m_normals in a right way at the end?
void MeshGenerator::FlyingEdgesAlgorithmPass4()
{
size_t nx = m_dim.x();
size_t ny = m_dim.y();
size_t nz = m_dim.z();
for(size_t k = 0; k != nz -1; ++k) {
for(size_t j = 0; j != ny-1; ++j)
{
// find adjusted trim values
size_t xl, xr;
calcTrimValues(xl, xr, j, k); // xl, xr set in this function
if(xl == xr)
continue;
size_t triIdx = triCounter[k*(ny-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];
gridEdge const& ge1 = gridEdges[k* m_dim.y() + j + 1];
gridEdge const& ge2 = gridEdges[(k+1)* m_dim.y() + j];
gridEdge const& ge3 = gridEdges[(k+1)* m_dim.y() + j + 1];
size_t x0counter = 0;
size_t y0counter = 0;
size_t z0counter = 0;
size_t x1counter = 0;
size_t z1counter = 0;
size_t x2counter = 0;
size_t y2counter = 0;
size_t x3counter = 0;
bool isYEnd = (j == ny-2);
bool isZEnd = (k == nz-2);
for(size_t i = xl; i != xr; ++i)
{
bool isXEnd = (i == nx-2);
unsigned char caseId = curCubeCaseIds[i];
if(caseId == 0 || caseId == 255)
{
continue;
}
const bool* isCutCase = isCut[caseId]; // size 12
std::array<std::array<float, 3>, 8> pointCube = m_cube->getPosCube(i, j, k);
std::array<float, 8> isovalCube = m_cube->getValsCube(i, j, k);
std::array<std::array<float, 3>, 8> gradCube = m_cube->getGradCube(i, j, k);
// Add Points and normals.
// Calculate global indices for triangles
std::array<size_t, 12> globalIdxs;
if(isCutCase[0])
{
// points-> array<vector3f>
size_t idx = ge0.xstart + x0counter;
std::array<float, 3> interpolatedPoint = interpolateOnCube(pointCube, isovalCube, 0);
std::array<float, 3> interpolatedNormal = interpolateOnCube(gradCube, isovalCube, 0);
points[idx] = Vector3f(interpolatedPoint[0], interpolatedPoint[1], interpolatedPoint[2]);
normals[idx] = Vector3f(interpolatedNormal[0], interpolatedNormal[1], interpolatedNormal[2]);
globalIdxs[0] = idx;
++x0counter;
}
if(isCutCase[3])
{
size_t idx = ge0.ystart + y0counter;
std::array<float, 3> interpolatedPoint = interpolateOnCube(pointCube, isovalCube, 3);
std::array<float, 3> interpolatedNormal = interpolateOnCube(gradCube, isovalCube, 3);
points[idx] = Vector3f(interpolatedPoint[0], interpolatedPoint[1], interpolatedPoint[2]);
normals[idx] = Vector3f(interpolatedNormal[0], interpolatedNormal[1], interpolatedNormal[2]);
globalIdxs[3] = idx;
++y0counter;
}
if(isCutCase[8])
{
size_t idx = ge0.zstart + z0counter;
std::array<float, 3> interpolatedPoint = interpolateOnCube(pointCube, isovalCube, 8);
std::array<float, 3> interpolatedNormal = interpolateOnCube(gradCube, isovalCube, 8);
points[idx] = Vector3f(interpolatedPoint[0], interpolatedPoint[1], interpolatedPoint[2]);
normals[idx] = Vector3f(interpolatedNormal[0], interpolatedNormal[1], interpolatedNormal[2]);
globalIdxs[8] = idx;
++z0counter;
}
if(isCutCase[1])
{
size_t idx = ge0.ystart + y0counter;
if(isXEnd)
{
std::array<float, 3> interpolatedPoint = interpolateOnCube(pointCube, isovalCube, 1);
std::array<float, 3> interpolatedNormal = interpolateOnCube(gradCube, isovalCube, 1);
points[idx] = Vector3f(interpolatedPoint[0], interpolatedPoint[1], interpolatedPoint[2]);
normals[idx] = Vector3f(interpolatedNormal[0], interpolatedNormal[1], interpolatedNormal[2]);
// y0counter counter doesn't need to be incremented
// because it won't be used again.
}
globalIdxs[1] = idx;
}
if(isCutCase[9])
{
size_t idx = ge0.zstart + z0counter;
if(isXEnd)
{
std::array<float, 3> interpolatedPoint = interpolateOnCube(pointCube, isovalCube, 9);
std::array<float, 3> interpolatedNormal = interpolateOnCube(gradCube, isovalCube, 9);
points[idx] = Vector3f(interpolatedPoint[0], interpolatedPoint[1], interpolatedPoint[2]);
normals[idx] = Vector3f(interpolatedNormal[0], interpolatedNormal[1], interpolatedNormal[2]);
// z0counter doesn't need to in incremented.
}
globalIdxs[9] = idx;
}
if(isCutCase[2])
{
size_t idx = ge1.xstart + x1counter;
if(isYEnd)
{
std::array<float, 3> interpolatedPoint = interpolateOnCube(pointCube, isovalCube, 2);
std::array<float, 3> interpolatedNormal = interpolateOnCube(gradCube, isovalCube, 2);
points[idx] = Vector3f(interpolatedPoint[0], interpolatedPoint[1], interpolatedPoint[2]);
normals[idx] = Vector3f(interpolatedNormal[0], interpolatedNormal[1], interpolatedNormal[2]);
}
globalIdxs[2] = idx;
++x1counter;
}
if(isCutCase[10])
{
size_t idx = ge1.zstart + z1counter;
if(isYEnd)
{
std::array<float, 3> interpolatedPoint = interpolateOnCube(pointCube, isovalCube, 10);
std::array<float, 3> interpolatedNormal = interpolateOnCube(gradCube, isovalCube, 10);
points[idx] = Vector3f(interpolatedPoint[0], interpolatedPoint[1], interpolatedPoint[2]);
normals[idx] = Vector3f(interpolatedNormal[0], interpolatedNormal[1], interpolatedNormal[2]);
}
globalIdxs[10] = idx;
++z1counter;
}
if(isCutCase[4])
{
size_t idx = ge2.xstart + x2counter;
if(isZEnd)
{
std::array<float, 3> interpolatedPoint = interpolateOnCube(pointCube, isovalCube, 4);
std::array<float, 3> interpolatedNormal = interpolateOnCube(gradCube, isovalCube, 4);
points[idx] = Vector3f(interpolatedPoint[0], interpolatedPoint[1], interpolatedPoint[2]);
normals[idx] = Vector3f(interpolatedNormal[0], interpolatedNormal[1], interpolatedNormal[2]);
}
globalIdxs[4] = idx;
++x2counter;
}
if(isCutCase[7])
{
size_t idx = ge2.ystart + y2counter;
if(isZEnd)
{
std::array<float, 3> interpolatedPoint = interpolateOnCube(pointCube, isovalCube, 7);
std::array<float, 3> interpolatedNormal = interpolateOnCube(gradCube, isovalCube, 7);
points[idx] = Vector3f(interpolatedPoint[0], interpolatedPoint[1], interpolatedPoint[2]);
normals[idx] = Vector3f(interpolatedNormal[0], interpolatedNormal[1], interpolatedNormal[2]);
}
globalIdxs[7] = idx;
++y2counter;
}
if(isCutCase[11])
{
size_t idx = ge1.zstart + z1counter;
if(isXEnd and isYEnd)
{
std::array<float, 3> interpolatedPoint = interpolateOnCube(pointCube, isovalCube, 11);
std::array<float, 3> interpolatedNormal = interpolateOnCube(gradCube, isovalCube, 11);
points[idx] = Vector3f(interpolatedPoint[0], interpolatedPoint[1], interpolatedPoint[2]);
normals[idx] = Vector3f(interpolatedNormal[0], interpolatedNormal[1], interpolatedNormal[2]);
}
globalIdxs[11] = idx;
}
if(isCutCase[5])
{
size_t idx = ge2.ystart + y2counter;
if(isXEnd and isZEnd)
{
std::array<float, 3> interpolatedPoint = interpolateOnCube(pointCube, isovalCube, 5);
std::array<float, 3> interpolatedNormal = interpolateOnCube(gradCube, isovalCube, 5);
points[idx] = Vector3f(interpolatedPoint[0], interpolatedPoint[1], interpolatedPoint[2]);
normals[idx] = Vector3f(interpolatedNormal[0], interpolatedNormal[1], interpolatedNormal[2]);
// y2 counter does not need to be incremented.
}
globalIdxs[5] = idx;
}
if(isCutCase[6])
{
size_t idx = ge3.xstart + x3counter;
if(isYEnd and isZEnd)
{
std::array<float, 3> interpolatedPoint = interpolateOnCube(pointCube, isovalCube, 6);
std::array<float, 3> interpolatedNormal = interpolateOnCube(gradCube, isovalCube, 6);
points[idx] = Vector3f(interpolatedPoint[0], interpolatedPoint[1], interpolatedPoint[2]);
normals[idx] = Vector3f(interpolatedNormal[0], interpolatedNormal[1], interpolatedNormal[2]);
}
globalIdxs[6] = idx;
++x3counter;
}
// 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]);
}
}
}}
// qDebug() << "pass 4 ended";
}
Please if you have any suggestions do let me know.
@will.schroeder Any help?
It’s almost certainly an indexing problem. Managing all of the offsets and indices is what makes the algorithm tricky to implement and why my hair turned gray 
I recommend that you start with a simple volume 3x3x3 with just the single center point val==1, and all other points ==0. Make sure that is working. Then build up from there: 5x5x5 volumes, and also move the point with val==1 towards the various boundaries to make sure all of that works.
1 Like
Thanks so much @will.schroeder sir. I wasn’t aware that you only have wrote the research papers. Really a pleasure of mine that you replied. I have completed the algorithm and it’s working perfectly.
One last doubt, Do you think we can remove the artifacts by laplacian smoothing? I used to remove the artifacts and make it smooth via laplacian smoothing when i used marching cubes.
Can we do the same with flying-edges as well?
I prefer vtkWindowedSincPolyDataFilter. The laplacian smoother tends to “shrink” the surface.