5

I have a big FEM model from where I can get the "surface" of the model, say the elements and vertex that define the surface of that FEM model. For plotting purposes (nice plots are always a win!) I want to plot it nicely. My approach is just to use 

lungs.Vertex=vtx;
lungs.Faces=fcs;
patch(lungs,'facecolor','r','edgecolor','none')

NOTE: I need edgecolor none, as this is 4D data and different FEM have different triangulation, if edges are plotted user will be dizzy.

enter image description here

However this will output everything in a really plain red color, which is not nice (as it can not show the complexity of the figure, which are lungs, for the attentive to details).

therefore I decided to use ligthing:

camlight; camlight(-80,-10); lighting phong;

But again, this is not entirely correct. Actually it seems that the patch nromals are not computed correctly by Matlab.

enter image description here

My supposition is that maybe the patches are not always defined counter-clockwise and therefore some normals go to the wrong direction. However that is something its not straightforward to check.

Anyone has a similar problem, or ahint of how should I aproach this problem in order to have a nice surface ploted here?

EDIT

Just for the shake of plotting, here is the result obtained with @magnetometer answer:

enter image description here

Ander Biguri
  • 35,140
  • 11
  • 74
  • 120

1 Answers1

3

If your model gives you outward oriented normals, you can re-sort the faces of your model so that Matlab can properly calculate it's own normals. The following function works if you have triangular faces and outward oriented normals: 

function [FaceCor,nnew]=SortFaces(Faces,Normals,Vertices)
FaceCor=Faces;
nnew=Normals*0;
for jj=1:size(Faces,1)
    v1=Vertices(Faces(jj,3),:)-Vertices(Faces(jj,2),:);
    v2=Vertices(Faces(jj,2),:)-Vertices(Faces(jj,1),:);

    nvek=cross(v2,v1); %calculate normal vectors
    nvek=nvek/norm(nvek); 
    nnew(jj,:)=nvek;
    if dot(nvek,Normals(jj,:))<0
        FaceCor(jj,:)=[Faces(jj,3) Faces(jj,2) Faces(jj,1)]; 
        nnew(jj,:)=-nvek;
    end

end

If your FEM model doesn't give you outward directed normals, one way could be to reconstruct the surface using e.g. a crust algorithm that gives you outward directed normals or correctly oriented patches.

EDIT: As you don't have normals, the only solution that comes to my mind is to reconstruct the surface. This implementation of the crust algorithm has worked well for me in the past. All you need to do is:

[FacesNew,NormalsNew]=MyRobustCrust(Vertices);

If I remember correctly, FacesNew are not yet oriented counterclockwise, but you can use the SortFaces algorithm I posted above to correct for this, as you now have correctly oriented face normals, i.e. run:

[FaceCor,~]=SortFaces(FacesNew,NormalsNew,Vertices)

If you use Matlab's reducepatch (e.g. reducedmodel=reducepatch(fullmodel,reduction); ) to reduce the number of vertices, you will have to reconstruct the surface again, as reducepatch doesn't seem to keep the correct orientation of the patches.

Ander Biguri
  • 35,140
  • 11
  • 74
  • 120
magnetometer
  • 646
  • 4
  • 12
  • I definitely don't have normals, thats the problem. If all the normals were "outward oriented" or "inward oriented" then there wont be any problem neither, just reorient them and thts it. But it seems (from the picture) that some go to the inside and some to the outside, making ligthing wrong. – Ander Biguri Sep 23 '14 at 08:58
  • Very very interesting functions indeed. However Im still having problem. I have "face normals" but it seems matlab wants "vertex normals" for setting ligthing. How do you suggest I approach this? Vertexnormals=mean(facenormals_with_that_vertex)? – Ander Biguri Sep 23 '14 at 14:53
  • I did it like that and got amazing results. Thanks a lot, this is a very nice answer. – Ander Biguri Sep 23 '14 at 15:13
  • 1
    Averaging the face normals would be one possibility. Depending on which version of Matlab you are using, there might be other possibilities. At least Matlab 2014 seems to have built-in functions that could help: ``TR = triangulation(T,P);`` and ``VN = vertexNormal(TR)`` could be worth to look at, however I can test them right now. – magnetometer Sep 23 '14 at 18:21
  • Im still with 2013b, however averaging just did the work, as you can see in the edit of my question. Its not perfect, but perfect enough to impress a supervisor or two, so it does the job! thanks again ;) – Ander Biguri Sep 23 '14 at 18:38
  • I think it can can handle non-triangular faces if you replace `[Faces(jj,3) Faces(jj,2) Faces(jj,1)]` with `fliplr(Faces(jj,:))`. – saastn Jan 10 '17 at 20:24