Issue connectivity faces with isosurface

Hi,

Say I have a gaussian-like function from which I extract the 0.5 isosurface. When I look at the list of faces, I see that the following faces are not necessarily connected i.e. the face(i,:) does not necessarily have two indices in common with the face(i+1,:). Does anyone know an efficient way to order that list in order to get this kind of connectivity property ? Thanks !

Code to see the issue:

x=linspace(-1,1,100);
[X,Y,Z]=meshgrid(x,x,x);
gaus = exp(-(X.^2+Y.^2+Z.^2)/2/0.1^2);
fv = isosurface(X,Y,Z,gaus,0.5);
hold on
for i=1:1:100
    p1 = fv.vertices(fv.faces(i,1),:);
    p2 = fv.vertices(fv.faces(i,2),:);
    p3 = fv.vertices(fv.faces(i,3),:);
    plot3([p1(1) p2(1)],[p1(2) p2(2)],[p1(3) p2(3)])
    plot3([p1(1) p3(1)],[p1(2) p3(2)],[p1(3) p3(3)])
    plot3([p3(1) p2(1)],[p3(2) p2(2)],[p3(3) p2(3)])
    pause
end

[Ouatehaouks]

Anyone who's an expert on faces and vertices could help ? I feel I could use ismember and filter like this but it's going to take ages before re-ordering everything....maybe there's some magic trick function available ?

[Walter Roberson]

Being able to do what you are asking for is equivalent to saying: if you create a Graph of the face connectivity, then is there a Hamiltonian path on that graph?

There is for at least some simple regular solids, but I do not know if there is in general.

I suspect that if you were to start with a graph that did not have a Hamiltonian that you could work back to a solid, but at the moment I do not know if you could find an example in 3d.

[Ouatehaouks]

The hamiltonian path connects all vertices in a unique path but it does not connect the triangles. I realize my code example was wrong so I edited it. In this example all triangles share at least one edge with another so all triangles share two points. This means that there must be a way to sorte the faces so that faces(i,:) and faces(i+1,:) have two indices in common. If I think of a dirty way to sort faces it would be

sorted_faces = ones(size(fv.faces))
sorted_faces(1,:) = fv.faces(1,:);
for i=2:size(fv.faces,1)
    if length(find(ismember(fv.faces(i,:),fv.faces(i-1,:))))==2
        sorted_faces(i,:) = fv.faces(i,:);
    else
        for j=i:size(fv.faces,2)
           if length(find(ismember(fv.faces(i,:),fv.faces(i-1,:))))==2
               sorted_faces(i,:) = fv.faces(i,:);
           end
        end
    end
end

This code will definitly not work but the idea is to loop, find the triangle that share one vertex with i and if not then find another one that shares an edges with triangle i. You would then to remove tha triangle j from the list to search for the next triangle etc...very heavy and slow.

[Walter Roberson]

You missed what I said about taking the graph of the face connectivity.

For example for a cube with faces L R U D F B then U (up) is connected to L (left), R (right), F (front), B (back) but not D (down)

%      U D L R F B
adj = [0 0 1 1 1 1;
       0 0 1 1 1 1;
       1 1 0 0 1 1;
       1 1 0 0 1 1;
       1 1 1 1 0 0;
       1 1 1 1 0 0];
G = graph(adj);
plot(G)

and that graph does have hamiltonian paths, including Up Left Front Down Right Back .

Up shares an edge with Left, Left shares an edge with Front, Front shares an edge with Down, Down shares an edge with Right, Right shares an edge with Back .

So in the case of a simple cube it is possible to order the faces such that each face shares an edge with the next face.

Hmmm... I guess this might be different from what you asked; I will think more about this.

Ah yes, each face is described in full each time, so if it shares an edge with another face then the shared edge must be in the list of edges for the face that was constructed in step faces(i,:) and so Yes, it would appear in the list of edges for the next step faces(i+1,:) .

Therefore if you find a hamiltonian path on the face connection graph then you can construct the face list with the property you asked for.

Likewise, if you have a list of faces vertex indices that has the property you asked for of there being a shared edge between adjacent faces, then provided that no face is duplicated in the faces vertex list, then following the list of shared edges would correspond to traversing a hamiltonian on the face connection graph -- the hamiltonian would correspond to not visiting the same face twice (not drawing the same face twice). So the ordering you want to achieve is possible if and only if there is a hamiltonian for the face connection graph -- the graph of which face is connected to which face.

[Ouatehaouks]

oh indeed sorry I missed that. And that's a great idea ! I'll try the hamiltonian path on the face connectivity then. I can then make this more flexible perhaps and authorize the path to go over the same node twice if necessary. I'll try if that works and update the post. thx.

[Walter Roberson]

It appears that isosurface() only produces triangle meshes. If it does not duplicate faces, then possibly that has implications for a traversal order.

Pick a vertex and a starting face. Traverse the three triangles in (say) counter-clockwise order, and then flip "down" to the next level, and traverse in the same order along the bottom of what you had already traversed, then move to the next level and so on.

[Ouatehaouks]

what do you mean by level ? If it is the row in faces, I'm not sure this will be sufficient. When one launches my exemple, one can see that triangles appear connected from one another and then at some point it jumps to another region of my sphere which is completely disconnected, so I need to find another level that stays connected and it does not seem to be a matter of counterclockwise to me. Or maybe I misunderstood your point.

Answer 1

[Fabio Freschi]

Sorting of the faces does not guaratee the property you think exist. If you check, your isosurface is a closed surface: all edges are shared by two triangles

clear all, close all

% your code
x=linspace(-1,1,100);
[X,Y,Z]=meshgrid(x,x,x);
gaus = exp(-(X.^2+Y.^2+Z.^2)/2/0.1^2);
fv = isosurface(X,Y,Z,gaus,0.5);

% plot isosurface (graphic check)
figure, axis equal, view([1 1 1]);
p = patch(fv,'FaceColor','r');

% load data on triangulation class
DT = triangulation(fv.faces,fv.vertices);
% mesh edges
E = edges(DT);

% find triangles attached to edges
tri = edgeAttachments(DT,E);

% check the number of triangels attached to edges
nTri = cellfun(@length,tri);

% check if the all edges are shared by two triangles
all(nTri == 2)

ans = logical
   1

[Ouatehaouks]

Thank you for your help. The triangulation seems to be a good way, I did not know you could triangulate a 3D surface, I thought it would connect points in 3D a well...I'm not sure though how to use E and tri to sort fv.faces so that fv.faces(i,:) and fv.faces(i+1,:) have two indices in common ?

[Fabio Freschi]

I don't think it is possible to have that property for a generic mesh, if you want to label all triangles