Posts Tagged ‘mesh’

Technical report: A Cotangent Laplacian for Images as Surfaces

Monday, April 23rd, 2012

A Cotangent Laplacian for Images as Surfaces
We decided to write up a quick, 2 page technical report about some ideas we’ve had in the last year.

By embedding images as surfaces in a high dimensional coordinate space defined by each pixel’s Cartesian coordinates and color values, we directly define and employ cotangent-based, discrete differential-geometry operators. These operators define discrete energies useful for image segmentation and colorization.

Isoline plots on triangle meshes in Matlab

Friday, March 9th, 2012

I recently posted how to plot nice looking iso intervals of scalar fields on triangle meshes using matlab. Along with the intervals I gave an ad hoc way of also showing isolines. This works ok if your data has fairly uniform gradient but they’re not really isolines.

Here’s a little matlab function to output real isolines. Save it in a file called isoline.m:

function [LS,LD,I] = isolines(V,F,S,iso)
  % ISOLINES compute a list of isolines for a scalar field S defined on the
  % mesh (V,F)
  % [LS,LD,I] = isolines(V,F,S,iso)
  % Inputs:
  %   V  #V by dim list of vertex positions
  %   F  #F by 3 list of triangle indices
  %   S  #V list of scalar values defined on V
  %   iso #iso list of iso values
  % Outputs:
  %   LS  #L by dim list of isoline start positions
  %   LD  #L by dim list of isoline destination positions
  %   I  #L list of indices into iso revealing corresponding iso value
  % alternative: tracing isolines should result in smaller plot data (every
  % vertex only appears once

  % make sure iso is a ROW vector
  iso = iso(:)';
  % number of isolines
  niso = numel(iso);

  % number of domain positions
  n = size(V,1);
  % number of dimensions
  dim = size(V,2);

  % number of faces
  m = size(F,1);

  % Rename for convenience
  S1 = S(F(:,1),:);
  S2 = S(F(:,2),:);
  S3 = S(F(:,3),:);

  % t12(i,j) reveals parameter t where isovalue j falls on the line from
  % corner 1 to corner 2 of face i
  t12 = bsxfun(@rdivide,bsxfun(@minus,iso,S1),S2-S1);
  t23 = bsxfun(@rdivide,bsxfun(@minus,iso,S2),S3-S2);
  t31 = bsxfun(@rdivide,bsxfun(@minus,iso,S3),S1-S3);

  % replace values outside [0,1] with NaNs
  t12( (t12<-eps)|(t12>(1+eps)) ) = NaN;
  t23( (t23<-eps)|(t23>(1+eps)) ) = NaN;
  t31( (t31<-eps)|(t31>(1+eps)) ) = NaN;

  % masks for line "parallel" to 12
  l12 = ~isnan(t23) & ~isnan(t31);
  l23 = ~isnan(t31) & ~isnan(t12);
  l31 = ~isnan(t12) & ~isnan(t23);

  % find non-zeros (lines) from t23 to t31
  [F12,I12,~] = find(l12);
  [F23,I23,~] = find(l23);
  [F31,I31,~] = find(l31);
  % indices directly into t23 and t31 corresponding to F12 and I12
  %ti12 = sub2ind(size(l12),F12,I12);
  %ti23 = sub2ind(size(l23),F23,I23);
  %ti31 = sub2ind(size(l31),F31,I31);
  % faster sub2ind
  ti12 = F12+(I12-1)*size(l12,1);
  ti23 = F23+(I23-1)*size(l23,1);
  ti31 = F31+(I31-1)*size(l31,1);

  % compute actual position values
  LS = [ ...
    ... % average of vertex positions between 2 and 3
    bsxfun(@times,1-t23(ti12), V(F(F12,2),:)) + ...
    bsxfun(@times,  t23(ti12), V(F(F12,3),:)) ...
    ; ... % average of vertex positions between 2 and 3
    bsxfun(@times,1-t31(ti23), V(F(F23,3),:)) + ...
    bsxfun(@times,  t31(ti23), V(F(F23,1),:)) ...
    ;... % average of vertex positions between 2 and 3
    bsxfun(@times,1-t12(ti31), V(F(F31,1),:)) + ...
    bsxfun(@times,  t12(ti31), V(F(F31,2),:)) ...
  LD = [...
    ... % hcat with average of vertex positions between 3 and 1
    bsxfun(@times,1-t31(ti12), V(F(F12,3),:)) + ...
    bsxfun(@times,  t31(ti12), V(F(F12,1),:)) ...
    ;... % hcat with average of vertex positions between 3 and 1
    bsxfun(@times,1-t12(ti23), V(F(F23,1),:)) + ...
    bsxfun(@times,  t12(ti23), V(F(F23,2),:)) ...
    ;... % hcat with average of vertex positions between 3 and 1
    bsxfun(@times,1-t23(ti31), V(F(F31,2),:)) + ...
    bsxfun(@times,  t23(ti31), V(F(F31,3),:)) ...

  % I is just concatenation of each I12
  I = [I12;I23;I31];


Then if you have a mesh in (V,F) and so scalar data in S. You can display isolines using:

axis equal;
[LS,LD,I] = isolines(V,F,S,linspace(min(S),max(S),nin+1));
hold on;
plot([LS(:,1) LD(:,1)]',[LS(:,2) LD(:,2)]','k','LineWidth',3);
hold off
set(gcf, 'Color', [1,1,1]);
set(gca, 'visible', 'off');
%O = outline(F);
%hold on;
%plot([V(O(:,1),1) V(O(:,2),1)]',[V(O(:,1),2) V(O(:,2),2)]','-k','LineWidth',2);
%hold off

Uncomment the last lines to display the outline of the mesh. This produces:
woody matlab isolines
If you use the my anti-alias script you can get pretty good results, almost good enough for a camera-ready paper:
woody matlab isolines anti-aliased

Get curve sketch (pen tool) from user in MATLAB figure

Thursday, March 1st, 2012

Here’s a little function that asks the user for a curve drawn on the current figure’s current axis. Similar to a pen/pencil tool in Photoshop etc.

Save this in a file called get_pencil_curve.m:

function P = get_pencil_curve(f)
  % GET_PENCIL_CURVE Get a curve (sequence of points) from the user by dragging
  % on the current plot window
  % P = get_pencil_curve()
  % P = get_pencil_curve(f)
  % Inputs:
  %   f  figure id
  % Outputs:
  %   P  #P by 2 list of point positions

  % Get the input figure or get current one (creates new one if none exist)
  if nargin == 0 || isempty(f)
    f = gcf;

  % get axes of current figure (creates on if doesn't exist)
  a = gca;

  % set equal axis 
  axis equal;
  % freeze axis
  axis manual;
  % set view to XY plane

  set(gcf,'keypressfcn',        @onkeypress);
  % living variables
  P = [];
  p = [];

  % loop until mouse up or ESC is pressed
  done = false;

  % We've been also gathering Z coordinate, drop it
  P = P(:,1:2);

  % Callback for mouse press
  function ondown(src,ev)
    % Tell window that we'll handle drag and up events
    set(gcf,'windowbuttonmotionfcn', @ondrag);
    set(gcf,'windowbuttonupfcn',     @onup);

  % Callback for mouse drag
  function ondrag(src,ev)

  % Callback for mouse release
  function onup(src,ev)
    % Tell window to handle down, drag and up events itself

  function onkeypress(src,ev)
    % escape character id
    ESC = char(27);
    switch ev.Character
    case ESC
      error(['Unknown key: ' ev.Character]);

  function append_current_point()
    % get current mouse position
    cp = get(gca,'currentpoint');
    % append to running list
    P = [P;cp(1,:)];
    if isempty(p)
      % init plot
      hold on;
      p = plot(P(:,1),P(:,2));
      hold off;
      % update plot

  function finish()
    done = true;


Here I’m using the function ask the user to draw a curve (blue) and then immediately simplify it (red), then mesh the interior (green).
matlab mesh sketch

Display wireframe mesh in matlab and save as vector graphics

Saturday, January 14th, 2012

For our siggraph submission I wanted to overlay a 2D image with the wireframe of the mesh we use to deform it. In matlab I can display a 2D mesh (V,F) easily using:

axis equal;

If I save this as .eps then I get a nice vector graphics version of my mesh.

But if I don’t want the filled faces to show up in the vector output, then when I try to save the following figure as a .eps file:

axis equal;

matlab just pretends to make a postscript vector file but instead its just a .eps file containing a raster image.

Here’s what I do to get a .eps vector graphics postscript file (which can then be turned into a PDF since output to PDF is broken in matlab):

E = edges(F);
plot([V(E(:,1),1) V(E(:,2),1)]',[V(E(:,1),2) V(E(:,2),2)]','k-');
axis equal
set(gcf, 'Color', [1,1,1]);
set(gca, 'visible', 'off');

New ETHZ masters thesis project available: Server-client mesh processing

Sunday, January 1st, 2012

server client mesh processing eth masters project

Olga Sorkine and I will be hosting a master’s thesis project. The project, entitled Server-client mesh processing is now available, and we are eagerly awaiting applications.Many polygon mesh processing techniques are difficult to implement. These techniques often go unseen and unused beyond publication. Others are too computationally intensive for consumer-level laptops, so their use is limited to academic or commercial communities. Fortunately, polygon mesh processing techniques are well-suited for a server-client implementation. Since most techniques need only an input mesh and a few parameters, all computation may be done on a publicly available web server, relying on the client only to upload the input and download the output.

In this project, we will set up the necessary infrastructure on both the server and client end. The client side will be a light-weight, browser-based interface for viewing and uploading 3D surface meshes. The server side will be a clean API to many existing mesh processing implementations (e.g. smooth- ing, curvature computation, tetrahedral volume meshing). With this infrastructure, we will explore geometric data compression specific to each flavor of mesh processing in order to optimize the transfer of output data. We may also use this infrastructure to conduct large scale user-studies via Amazon’s Mechanical Turk. Such large scale user studies are rarely achieved in the mesh processing community. A server-client infrastructure should finally make human evaluation possible. Please don’t hesitate to contact me for more details.

Also, check out the full list of IGL projects.

Note: You will need to be at an ETH IP address to visit these links.

Perceived stretch from partial view of twist

Wednesday, August 31st, 2011

Here’s a twist applied to a bar.

But when the movie is cropped the viewer gets a sense that the bar is also being stretched, to the right.

Putting a texture (here just a the edges of the 3D model’s mesh), helps perceive the true deformation: twist.

It also helps reduce the perceived stretch effect on the original full-view movie.

Pseudo-color a mesh using random colors for multiple (weight) functions

Tuesday, August 23rd, 2011

Often it’s useful to visualize many functions on a mesh at the same time. Here I assign each function a random pseudocolor, then the color at each vertex of the mesh is just a blend of the function colors weighted by their respective function values. This is especially useful if your functions are already weight functions themselves (defined between 0 and 1). So, if your mesh is defined by V and F then you have a set of function values defined on the vertices W: #V by #functions, then you can use a previously posted function for randomly generating a #functions list of colors. Putting it all together with a matrix multiplication you have:


To make the rendering a little nicer you can use:

lighting gouraud
axis equal
shading interp

Here’s a visualization of bounded biharmonic weights on the Armadillo.

armadillo bbw pseudo color

And here’s an easy way to visualize segmenting the mesh based on the maximum function value:

% get indices of max value for each row

which produces something like:
armadillo bbw pseudo color max value segmentation

Bounded Biharmonic Weights: 2D MATLAB Demo

Monday, July 11th, 2011

woody matlab bbw demo I’ve finally released the matlab prototyping codebase for Bounded Biharmonic Weights. The code base runs without any prerequisite libraries on MATLAB 2011a or greater. With earlier versions of matlab, users should install MOSEK (which is available for free to academic users). This matlab package demos computing skinning weights automatically for a 2d shape. To start, add bbw_demo/ to your matlab path and issue:

>> bbw_demo

If you just want to start digging around in the weight computation code, look at:

  • boundary_conditions.m
  • biharmonic_bounded.m

Currently the bbw_demo.m file just shows basic weight computation and linear blend skinning deformation for point handles. The following features are supported in the codebase above, but are undocumented besides their file header comments:

  • Bone handles
  • Cage edges
  • Dual quaternion skinning
  • Automatic rotations at point handles via pseudo-edges
  • 3D weight computation over tetrahedral meshes
  • 2D/3D remeshing along bones and cage edges
  • …and more

The following things are not supported and probably never will be in the MATLAB demo. If I release the C++ demo then of course they will be available:

  • 3D interactive deformation
  • 2D texture mapping

contact me if you find any bugs or have any trouble using the code.

Update: Here’s a new C++ demo

Common mesh energies, sum notation and matrix notation

Monday, June 20th, 2011

I’m always re-deriving the matrix notation of common triangle mesh energies like the Dirichlet energy, Laplacian energy, or local rigidity energy. Here’s my usual strategy. I’ll assume that if your energy came from a continuous integral you can at least get it to a sum over mesh elements. Here I write my sums over vertices and their neighbors, but surely you could also sum over triangles.

mesh energies sum to matrix form
Pretty-printed pdf document

Uglier plain text, unicode version:

∑i∑j∈N(i) vi  = ∑i * #N(i) * v = VT * D * 1, 
  where D is diagonal with Dii = #N(i)
∑i∑j∈N(i) vj  = VT * A * 1, where A is the adjacency matrix
∑i∑j∈N(i) vi - vj = VT * D * 1 - VT * A * 1 = VT * (D-A) * 1

∑i∑j∈N(i) wij  = 1 * Aw * 1, where Aw adjacency with Awij = wij, if j∈N(i)
∑i∑j∈N(i) wij*vi  = ∑i vi * ∑j∈N(i) wij  = VT * Dw * 1, 
  where Dw diagonal with Dwii = ∑j∈N(i) wij
∑i∑j∈N(i) wij*vj  = VT * Aw * 1

∑i mi * ∑j∈N(i) wij vi = ∑i mi * vi * ∑j∈N(i) wij = VT * M * Dw * 1,
  where M is diagonal with Mii = wi
∑i mi * ∑j∈N(i) wij vj = VT * Aw * M * 1

∑i∑j∈N(i) vi² = VT * D * V
∑i∑j∈N(i) vivj = VT * A * V
∑i∑j∈N(i) vj² = VT * RD * V, where RD is diagonal, RDjj = ∑i 1 if j∈N(i)
  if j∈N(i) implies i∈N(j) then
  D = RD
  so ∑i∑j∈N(i) vj² = ∑i∑j∈N(i) vi² = VT * D * V
∑i∑j∈N(i) (vi - vj)² = ∑i∑j∈N(i) vi² - 2 * ∑i∑j∈N(i) vivj + ∑i∑j∈N(i) vj²
                     = 2 * ∑i∑j∈N(i) vi² - 2 * ∑i∑j∈N(i) vivj
                     = 2 * VT * D * V - 2 * VT * A * V
                     = 2 * VT * (D - A) * V
                     = 2 * VT * L * V

∑i∑j∈N(i) wij*vi²  = ∑i vi² * ∑j∈N(i) wij = VT * Dw * V
∑i∑j∈N(i) wij*vivj = VT * Aw * V
∑i∑j∈N(i) wij*vj² = VT * RDw * V, where DW is diagonal, RDwjj = ∑i wij if j∈N(i)
  if j∈N(i) implies i∈N(j) and wij = wji then 
  Dw = RDw
  so ∑i∑j∈N(i) wij*vj² = ∑i∑j∈N(i) wij*vi² = VT * Dw * V
∑i∑j∈N(i) wij(vi - vj)² = 
  = ∑i∑j∈N(i) wij*vi² - 2 * ∑i∑j∈N(i) wij*vivj + ∑i∑j∈N(i) wij*vj²
  = 2 * ∑i∑j∈N(i) wij*vi² - 2 * ∑i∑j∈N(i) wij*vivj
  = 2 * VT * Dw * V - 2 * VT * Aw * V
  = 2 * VT * (Dw - Aw) * V
  = 2 * VT * Lw * V

∑i mi * ∑j∈N(i) wij vi² = ∑i mi * vi² * ∑j∈N(i) wij = VT * M * Dw * V
∑i mi * ∑j∈N(i) wij vj² = VT * RDw * M * V
∑i mi * ∑j∈N(i) wij vivj = VT * Aw * M * V
∑i mi * ∑j∈N(i) wij(vi - vj)² = 
  = ∑i mi * ∑j∈N(i) wij*vi² 
    - 2 * ∑i mi * ∑j∈N(i) wij*vivj 
    + ∑i mi *∑j∈N(i) wij*vj²
  = 2 * ∑i mi * ∑j∈N(i) wij*vi² - 2 * ∑i mi * ∑j∈N(i) wij*vivj
  = 2 * VT * M * Dw * V - 2 * VT * Aw * M * V
  = 2 * VT * M * (Dw - Aw) * V
  = 2 * VT * M * Lw * V
  Note: let wij = cwij = cot(αij) + cot(βij) / mi and let mi = voronoi area at 
  vi, then the familar cotangent form of the discrete Dirichlet energy is
  EDirichlet = 2 * VT * M * Lcw * V

Split triangular prism into three tetrahedra

Wednesday, June 15th, 2011

triangular prism split into three tetrahedra small

Recently I began looking into tet-meshes again. Some colleagues and I were trying to wrap our heads around how many tetrahedra result from splitting up a triangular prism. we were trying to draw out the cases on the white board. Occasionally some one would say, “I’m sure it’s 3.” Then somebody else would draw another set of pictures and’d say “No, I’m sure it’s 4″.

Well, to have a definitive answer, a triangular prism may be split into three tetrahedra. If you don’t believe me, see the nice description (which actually lists all possible splits) in “How to Subdivide Pyramids, Prisms and Hexahedra into Tetrahedra” by J. Dompierre et al. in 1999. There is another interesting discussion in “Space-filling Tetrahedra in Euclidean Space” by D.M.Y. Sommerville in 1923, in which they consider the conditions under which these tetrahedra are congruent.