Posts Tagged ‘graphics’

Determine boundary faces from tetrahedral mesh

Wednesday, April 20th, 2011

Here’s a matlab function that takes a list of tetrahedron indices (4 indices to a row) and finds the triangles that are on the surface of the volume. It does this simply by finding all of the faces that only occur once.


function F = boundary_faces(T)
  % BOUNDARY_FACES
  % F = boundary_faces(T)
  % Determine boundary faces of tetrahedra stored in T
  %
  % Input:
  %  T  tetrahedron index list, m by 4, where m is the number of tetrahedra
  %
  % Output:
  %  F  list of boundary faces, n by 3, where n is the number of boundary faces
  %

  % get all faces
  allF = [ ...
    T(:,1) T(:,2) T(:,3); ...
    T(:,1) T(:,3) T(:,4); ...
    T(:,1) T(:,4) T(:,2); ...
    T(:,2) T(:,4) T(:,3)];
  % sort rows so that faces are reorder in ascending order of indices
  sortedF = sort(allF,2);
  % determine uniqueness of faces
  [u,m,n] = unique(sortedF,'rows');
  % determine counts for each unique face
  counts = accumarray(n(:), 1);
  % extract faces that only occurred once
  sorted_exteriorF = u(counts == 1,:);
  % find in original faces so that ordering of indices is correct
  F = allF(ismember(sortedF,sorted_exteriorF,'rows'),:);
end

With this you can easily determine the vertices of a tetmesh that are on the boundary:


% V are your vertex positions, T are your tet indices
% get boundary faces
F = boundary_faces(T);
% get boundary vertices
b = unique(F(:));
subplot(1,2,1);
% plot boundary positions
plot3(V(b,1),V(b,2),V(b,3),'.');
subplot(1,2,2);
% plot just  boundary faces
trisurf(F,V(:,1),V(:,2),V(:,3),'FaceAlpha',0.3)

Rotate a point around another point

Tuesday, February 22nd, 2011

Composing rotations and translations is one of the most important operations in computer graphics. One useful application is the ability to compose rotations and translations to rotate a point around another point.

rotate around other point

Here’s how I learned to do this. I find this method the most intuitive.

  1. Translate so that the other point is at the origin
  2. Rotate about the origin by however much you want
  3. Translate back so that the other point is back to its original position

rotate around other point

This works out nicely with matrix math since rotations around the origin are easily stored as 2×2 matrices:


                               / cos θ  -sin θ\
Rotation around origin by θ:  |                |
                               \ sin θ   cos θ/

If we call that matrix, R, then we can write the whole operation that rotates a point, a, around another point, b,, as:
R*(a-b) + b. Be careful to note the order of operations: (a-b) corresponds to step 1, then left multiply with R to step 2, and finally adding back b is step 3.

One convenient fact is that when we look at the transformation as a whole where T(a): a → R*(a-b) + b. Because T is just the composition of rotations and translations it can be decomposed into a single rotation followed by a single translation. Namely, T(a): a → R*(a-b) + b may be re-written as T(a): a → S(a) + t, for some rotation matrix S and some translation vector t. To see this just use the distributive law: R*(a-b) + b = R*a – R*b + b, then S = R and t = R*b + b.

So that gives a new derivation of the transformation that rotates a point, a, around another point, b,: R*a – R*b + b. Building an intuition as to why this works is a little tricky. We have just seen how it can be derived using linear algebra but actually seeing this version as each step is elusive. Turns out it helps to make a subtle change. Reverse the last two terms, so that you have: R*a + (b – R*b). Now intuitively we can follow the order of operations and build an intuition:

  1. Rotate first the point about the origin (since the other point is not the origin we’ve “missed” where we should have been rotated by a certain error amount)
  2. Rotate the other point about the origin by the same amout
  3. Translate by the difference between the original other point and the rotated other point (this is the error amount, because we know that rotating other pointabout itself shouldn’t change its position)

rotate around other point

Update: Here’s an image comparing these two compositions:

rotate around other point

Barycentric coordinates and point-triangle queries

Friday, February 4th, 2011

Recently I needed to test whether a bunch of points where on a given triangle in 3D. There are many ways of querying for this, but I found one that was for me both intuitive and convenient. The method uses barycentric coordinates.

Given a triangle made up of points A, B, and C you can describe any point, p, inside the triangle as a linear combination of A, B, and C: p = λa*A + λb*B + λc*C. That is to say, p can be reproduced as a weighted average of A, B and C. Where the λs are the weights.

One way to realize this is true is to notice that a triangle can be defined by two vectors, say A→B and A→C, and as long as these two vectors are not parallel we know from Linear Algebra that we can span the whole plane that they define (including the triangle A, B,C which lies on that plane).

abc) are called the “barycentric coordinates” of p. Another (more or less unused) name for these are the “area coordinates” of p. This is because λa, λb, and λc are very conveniently defined as:
λa = area(B,C,P)/area(A,B,C)
λb = area(C,A,P)/area(A,B,C)
λc = area(A,B,P)/area(A,B,C)

An important quality of barycentric coordinates is that they partition unity, that is: λabc = 1. To see this is true, notice that these sub triangles exactly cover the area of the whole triangle and since we divide be that whole triangle’s area we are normalizing the sub triangle areas such that they must sum to one: the barycentric coordinates are the percentages of the area covered by the sub triangles and the whole triangle area. Each sub triangle corresponds to the original triangle corner that is not a corner of the sub triangle. Thus, A corresponds to the sub triangle B, C, p and so forth.


As we move p inside the triangle, ABC, the barycentric coordinates shift. Notice that all coordinates are positive as long as p is inside the triangle.


If we move p onto an edge of the triangle then the area of the sub triangle corresponding to the corner opposite that edge becomes zero, thus the barycentric coordinate for p corresponding to that corner is zero.


If we move p onto a corner of the triangle then the sub triangle corresponding to that corner is the same (and thus has the as area) as the original triangle, so its corresponding barycentric coordinate to p is 1, the other areas and coordinates being 0.


If we move p outside of the triangle the total area covered by the sub triangles will be greater than the original triangle. This is easy to see as the sub triangles always cover the original triangle. We could use this to test whether p is inside the triangle ABC, but there actually a more elegant way that will be handy later.

If instead of computing regular (non-negative) area, we compute signed area then even when we move p outside of ABC the areas will sum to the original area of the triangle ABC. One way to think of signed area is by looking at the orientation of the triangles. Start with p inside ABC.


Imagine coloring the side of each triangle facing you green and the back side red. When we move p the triangles change shape, but as long as p stays inside the triangle we still see the same, green side.


When we pull p outside the triangle the triangle on the edge we crossed gets flipped, and now we see the back, red, side. Now we declare that if we can see the green side of the triangle, its area is positive and if we see the red side of the triangle its area is negative. Now, finally notice that the red side of the flipped triangle when p is outside ABC is always covered exactly by the other two, green triangles. Thus the negative area is cancelled out and we are only left with exactly the original ABC area.

So if we use signed area to compute barycentric coordinates of p, we can determine if p is inside, on, or outside a triangle by looking to see if all coordinates are positive, if any are 0 or if any are negative.

If our triangle is on the 2D plane then this sums up everything we might want to know about p and a given triangle ABC. But if ABC is an arbitrary triangle in 3D and p is also an arbitrary point, then we may also want to know if p whether p is on the plane of the triangle ABC, or if it is above or below the triangle.

We showed above that if p is on the plane of the triangle then the total un-signed area of the sub triangles is greater than the original triangle area and the total signed area of the sub triangles is always exactly equal to the area of the original triangle.

However, if we pull p off the plane of the sum of the total un-signed area of the sub triangles will necessarily be greater than the area of the original triangle. To see this notice that by pulling p off the triangle we are making a tetrahedron pABC with a certain amount of volume, where before pABC was flattened onto the base triangle ABC with no volume. Where ever p is in space it makes such a tetrahedron with ABC and we can always flatten (project) it to the plane of the triangle where we know the total area of the flattened triangles must cover (be greater than) the original triangle ABC. Now finally notice that “flattening” a triangle can never increase the triangle’s area (imagine holding a paper triangle in front of your face, your eyes project it onto the plane which you see. There’s no way of rotating the triangle so that the projected area (the area you see) is bigger than the actual area of the triangle)).

Finally, we just showed that pulling p away from the plane of the triangle increases the magnitude of each of the sub triangles’ area. We may also use this fact to see that the total signed area of these sub triangles is only ever equal to the area of the original triangle if we keep p on the triangle’s plane. If we pull p above the triangle (off its front side) then the total signed area is always greater than the original area, and if we pull p below the triangle (off its back side) then the total signed area is always less than the original area. This is harder to see and hopefully I will find a nice way of explaining it.

It’s important to note that the areas we get when we pull p off of the plane of the triangle ABC can no longer partition unity, in the sense that if we use these areas as above to define λa, λb, and λc, p will still equal λa*A + λb*B + λc*C, but λa + λb + λc will not equal 1.

Triangle mesh filled contour in MATLAB

Monday, December 6th, 2010

Here is matlab snippet that takes a 2D triangle mesh and a function and produces a pretty, filled contour plot in MATLAB. If your mesh vertices are stored in V, and your triangles in F, then use my previously posted code to get the boundary/outline edges of your mesh in O. Then you can create a contour plot of some function W defined over V using:


% resample W with a grid, THIS IS VERY SLOW if V is large
[Xr,Yr,Wr] = griddata(V(:,1),V(:,2),W,unique(V(:,1)),unique(V(:,2))');
% find all points inside the mesh, THIS IS VERY SLOW if V is large
IN = inpolygon(Xr,Yr,V(O,1),V(O,2));
% don't draw points outside the mesh
Wr(~IN) = NaN;
contourf(Xr,Yr,Wr)

And for a mesh like this one:
gingerbread man mesh plot

You get something like this:
gingerbread man contour 10 level sets

You can also approach a continuous contour with this:


contourf(Xr,Yr,Wr,200)
shading flat

gingerbread man contour 200 level sets

Compare this to what you can get from the much faster:


trisurf(F,V(:,1),V(:,2),W,'FaceColor','interp','EdgeColor','interp')
view(2)

which produces:
gingerbread man colored plot

The only problem with this method is that the boundary looks a little nasty because of the resampling. Nothing that can’t be fixed quickly in photoshop… Or even in MATLAB with something like:


line(V([O(:);O(1)],1),V([O(:);O(1)],2),'Color','k','LineWidth',6)

Which puts an ugly thick outline around the contour, but makes it look a little better…
gingerbread man contour 10 level sets with outline

Source

Find outline of triangle mesh and plot in MATLAB

Monday, December 6th, 2010

Here’s a little snippet to determine the edges along the boundary of a triangle mesh in matlab. Given that your vertex values are in V and your triangle indices are in F, this will fill O with a list of edges make up the outline of your mesh:


% Find all edges in mesh, note internal edges are repeated
E = sort([F(:,1) F(:,2); F(:,2) F(:,3); F(:,3) F(:,1)]')';
% determine uniqueness of edges
[u,m,n] = unique(E,'rows');
% determine counts for each unique edge
counts = accumarray(n(:), 1);
% extract edges that only occurred once
O = u(counts==1,:);

Then you can view the outline with, in 2D:


plot([V(O(:,1),1) V(O(:,2),1)]',[V(O(:,1),2) V(O(:,2),2)]','-')

and in 3D:


plot3([V(O(:,1),1) V(O(:,2),1)]',[V(O(:,1),2) V(O(:,2),2)]',[V(O(:,1),3) V(O(:,2),3)]','-')

See also: Display wireframe mesh in matlab and save as vector graphics

“Mixed Finite Elements for Variational Surface Modeling” project page

Monday, October 25th, 2010

point boundary conditions used to make bumps in flat domain

I’ve finally put a project page for the paper we presented this summer at SGP: “Mixed Finite Elements for Variational Surface Modeling” by Alec Jacobson, Elif Tosun, Olga Sorkine and Denis Zorin. So far I have the paper, slides from SGP, slides from my recent Disney Tech Talk at ETH and some videos and images. Hopefully some I will post some of the working MATLAB code base.

OpenGL render GL_POINTS as circles not squares

Wednesday, September 1st, 2010

Bugged me into staying late today. To make OpenGL render points (GL_POINTS) as circles instead of squares be sure to add:


glEnable(GL_POINT_SMOOTH);

aka point antialiasing.

Real time tester

Thursday, July 22nd, 2010


I’ve adapted an old applet to function as a real time tester. The idea being that you hear a new graphics algorithm is “real-time” or “interactive” because it runs in only “59ms” per frame, but your internal clock isn’t accurate enough to really know what that feels like. Here you can specify the simulated “solve” time and drag around the bezier curve as if hard-core math is going on between frames.

OpenGL texture map all white

Monday, June 28th, 2010

After a long battle today, I found a bug that was making all my textures render white. Here’s pseudo code of what I had:

Bugged code

GLwidget.h


class GLwidget : QGLWidget {
  GLWidget();
  ...
  void initializeGL();
  void paintGL();
  ...
}

GLwidget.cpp


GLWidget::GLWidget(){
  // code that generates texture
  makeMyTexture(...);
  // opengl calls that setup texture in GL land
  glGenTextures(1, &texName);
  glBindTexture(GL_TEXTURE_2D, texName);
  glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_REPEAT);
  ...
  glTexImage2D(GL_TEXTURE_2D, ...) ;
}
  ...
void GLWidget ::initializeGL(){
  // opengl initialization code
}
void GLWidget :: paintGL(){
  // code that actually uses the texture
  glEnable(GL_TEXTURE_2D)
  ...
  glDisable(GL_TEXTURE_2D)
}

My GLWidget class inherited from a class that called initializeGL and paintGL at the appropriate times. The only problem was that I was building and setting my texture in the class constructor before initializeGL() was ever called. Then when it was, it wiped out my texture.

The correct order should be:

Correct

GLwidget.cpp


GLWidget::GLWidget(){
  ...
}
void GLWidget ::initializeGL(){
  // opengl initialization code
  ...
  // code that generates texture
  makeMyTexture(...);
  // opengl calls that setup texture in GL land
  glGenTextures(1, &texName);
  glBindTexture(GL_TEXTURE_2D, texName);
  glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_REPEAT);
  ...
  glTexImage2D(GL_TEXTURE_2D, ...) ;
}
void GLWidget :: paintGL(){
  // code that actually uses the texture
  glEnable(GL_TEXTURE_2D)
  ...
  glDisable(GL_TEXTURE_2D)
}

Initialize your texture after you have initialized opengl.

Simple MATLAB slerp

Monday, April 12th, 2010

Interpolating between two vectors can catch you off guard. If you just linearly interpolate (x,y,z) values you could end up with a vanishing interpolated vector. Just imagine interpolating a 2-D vector from [-1,0] to [1,0]. Halfway you’ve got [0,0].

vector interpolation

One way to do this a little better is to interpolate a rotation between the vectors. The natural path to take is the shortest one which happens to also be the shortest path (of two unit vectors) of their respect points along the unit sphere. This is called spherical linear interpolation or slerp.

slerp

There’s more complicated ways to do this but in MATLAB I can simply do:


function [c] = slerp(a,b,t)
  angle = acos( dot(a,b) );
  % easy degenerate case
  if(0==angle)
    c = a;
  % hard case
  elseif(pi==angle) 
    error('SLERP: angle between vectors cannot be exactly PI...');
  else
    c = (sin((1.0-t)*angle)/sin(angle))*a + (sin(t*angle)/sin(angle))*b; 
  end
end

Note: that slerping is not well defined for opposite parallel vectors. This makes sense geometrically because their respective points on the unit sphere have infinitely many, all equally minimal paths between them. E.g. the north pole has infinitely many shortest paths along the earth to the south pole.

In the end, I lerp the magnitudes of the two vectors then use that to scale the slerp of their unit vectors:


c = ( (1.0-t) * norm(a) + t* norm(b) ) * slerp( a./norm(a) , b./norm(b) , t);

where t is your interpolation parameter.

Slightly more efficiently:


mag_a = norm(a); mag_b = norm(b); c = ((1.0-t)*mag_a + t*mag_b) * slerp(a./mag_a,b./mag_b,t);

Note: Careful using norm it is not vectorized as you might think. I.e. it doesn’t return a row vector of scalars if you pass it a list of vectors (instead the matrix norm of that list treated as a matrix).