Posts Tagged ‘triangle’

How does matlab’s hypot function work?

Saturday, January 7th, 2017

The hypot function in matlab purports to be more numerically stable at computing the hypotenuse of a (right-)triangle in 2D. The example the help hypot gives is:

a = 3*[1e300 1e-300];
b = 4*[1e300 1e-300];
c1 = sqrt(a.^2 + b.^2)
c2 = hypot(a,b)

where you see as output:

c1 =

   Inf     0


c2 =

       5e+300       5e-300

this is a compiled built-in function so you can’t just open hypot to find out what its doing. It might just be pre-dividing by the maximum absolute value. There’s probably a better reference for this, but I found it in: “Vector length and normalization difficulties” by Mike Day.

Continuing this example:

m = max(abs([a;b]))
c3 = m.*sqrt((a./m).^2 + (b./m).^2)


c3 =

       5e+300       5e-300

While matlab’s hypot only accepts two inputs, pre-dividing by the maximum obviously extends to any dimension.

Triangle area in nD

Saturday, January 7th, 2017

My typical go-to formula for the area of a triangle with vertices in nD is Heron’s formula. Actually this is a formula for the area of a triangle given its side lengths. Presumably it’s easy to compute side lengths in any dimension: e.g., a = sqrt((B-C).(B-C)).

Kahan showed that the vanilla Heron’s formula is numerically poor if your (valid) triangle side lengths are given as floating point numbers. He improved this formula by sorting the side-lengths and carefully rearranging the terms in the formula to reduce roundoff.

However, if you’re vertices don’t actually live in R^n but rather F^n where F is the fixed precision floating-point grid, then computing side lengths will in general incur some roundoff error. And it may so happen (and it does so happen) that this roundoff error invalidates the triangle side lengths. Here, invalid means that the side lengths don’t obey the triangle inequality. This might happen for nearly degenerate (zero area) triangles: one (floating point) side length ends up longer than the sum of the others. Kahan provides an assertion to check for this, but there’s no proposed remedy that I could find. One could just declare that these edge-lengths must of come from a zero-area triangle. But returning zero might let the user happily work with very nonsensical triangle side lengths in other contexts not coming from an embedded triangle. You could have two versions: one for “known embeddings” (with assertions off and returning 0) and one for “intrinsic/metric only” (with assertions on and returning NaN). Or you could try to fudge in an epsilon a nd hope you choose it above the roundoff error threshold but below the tolerance for perceived “nonsense”. These are messy solutions.

An open question (open in the personal sense of “I don’t know the answer”) is what is the most robust way to compute triangle area in nD with floating point vertex positions. A possible answer might be: compute the darn floating point edge lengths, use Kahan’s algorithm and replace NaNs with zeros. That seems unlikely to me because (as far as I can tell) there are 4 sqrt calls, major sources of floating point error.

Re-reading Shewchuk’s robustness notes, I saw his formula for triangle area for points A,B,C in 3D:

Area3D(A,B,C) = sqrt( Area2D(Axy,Bxy,Cxy)^2 + Area2D(Ayz,Byz,Cyz)^2 + Area2D(Azx,Bzx,Czx)^2 ) / 2

where Area2D(A,B,C) is computed as the determinant of [[A;B;C] 1]. This formula reduces the problem of computing 3D triangle area into computing 2D triangle area on the “shadows” (projections) of the triangle on each axis-aligned plane. This lead me to think of a natural generalization to nD:

AreaND(A,B,C) = sqrt( ∑_{i<j} ( Area2D(Aij,Bij,Cij)^2 ) / 2

This formula computes the area of the shadow on all (N choose 2) axis-aligned planes. Since the sqrt receives a sum of squares as in argument there’s no risk of getting a NaN. There’s a clear drawback that this formula is O(N^2) vs Heron’s/Kahan’s O(N).

Stable triangle normals

Wednesday, February 11th, 2015

Here’s one way to compute a triangle’s normal using floating point arithmetic that is stable with respect to the order the triangle’s vertices are given in. To demonstrate the issue, consider:

% Vertex positions
A = [1/8 1 1/6];
B = [1/3 1/5 1/7];
C = [1/4 1/9 1/2];
% normals via cross product of different edge vectors
NB = cross(A-B,C-B);
NC = cross(B-C,A-C);
NA = cross(C-A,B-A);


ans =
            0            0  -2.7756e-17
ans =
   5.5511e-17   1.3878e-17            0

Computing the normals three different ways, I get three different floating point errors. The right thing to do I guess is analyze the vertex positions and design a normal computation that minimizes floating point error (the solution hopefully being stable via uniqueness). I took a lazier route and take a careful average of the three solutions. The might also have the effect of cutting down floating point error a bit. But at least it’s stable.

Unfortunately it’s not enough to naively sum N1, N2, and N3 entry-wise and divide by 3. The summation could incur roundoff error depending on the order of the three arguments. So first we must sort the arguments. If going to the trouble of sorting we might as well sort in a way that reduces roundoff error: descending order of absolute value. The matlab syntax for this is a little clunky:

NABC = [NA(:) NB(:) NC(:)];
[~,I] = sort([NABC],2,'descend');
sNABC = reshape( ...
N = 1/3*reshape(sum(sNABC,2),shape);

Solid sweeps along line segments

Saturday, February 7th, 2015

I made a nice realization today that the boolean code I’ve been working on can also be used to compute “solid sweeps” of meshes along line segments: the Minkowski sum of a polyhedron and a line segment.

Check out linear_sweep.m in gptoolbox. It works in 2D (using Triangle and the winding number to classify the boundary of the swept region):

octopus swept area

and in 3D:

knight swept volume

the 3D version is currently only robust if the sweep vector goes outside the bounding box of the input shape. In those cases, I’m using a technique analogous to my robust booleans and it works quite well. But if the sweep vector is small, I’m resorting to using TetGen and computing winding numbers: TetGen will often fail. Hopefully this limitation is not for long, since I have a clear idea how to handle small vectors. Just need to code it up.

Update: The 3D version is now robust, too.

A Simple Method for Correcting Facet Orientations in Polygon Meshes Based on Ray Casting

Sunday, November 23rd, 2014

correct facet orientations on a bike, cell phone, chair, and truck

We’ve finally published our paper A Simple Method for Correcting Facet Orientations in Polygon Meshes Based on Ray Casting in the Journal of Computer Graphics Tools. The paper was written by Kenshi Takayama, Alec Jacobson, Ladislav Kavan, and Olga Sorkine-Hornung.

Abstract: We present a method for fixing incorrect orientations of facets in an input polygon mesh, a problem often seen in popular 3D model repositories, such that the front side of facets is visible from viewpoints outside of a solid shape represented or implied by the mesh. As opposed to previously proposed methods which are rather complex and hard to reproduce, our method is very simple, only requiring sampling visibilities by shooting many rays. We also propose a simple heuristic to handle interior facets that are invisible from exterior viewpoints. Our method is evaluated extensively with the SHREC Generic 3D Warehouse dataset containing 3168 manually designed meshes, and is demonstrated to be very effective.

You can find an implementation of this method in the embree extension of libigl: igl/embree/reorient_facets_raycast.h. If you store your mesh by its vertices in rows of a #V by 3 matrix V and triangle indices in the rows of a #F by 3 matrix F. Then you can quickly reorient your triangles to all point outward consistently with a single function call:

#include <igl/embree/reorient_facets_raycast.h>
Eigen::MatrixXI FF; // #F by 3 list of output triangle indices, some rows potentially reversed
Eigen::VectorXi I; // #F list of booleans revealing whether facet was reversed

As a preprocessor to our generalized winding numbers, this forms a powerful tool for determining the inside from outside for arbitrary meshes. This is especially important for creating volumetric tetrahedral meshes.

Triangle mesh for image as surface

Sunday, October 26th, 2014

Here’s a little chuck of matlab code I use to make a surface mesh out of an image:

% load image as grayscale
im = rgb2gray(imresize(im2double(imread('hans-hass.jpg')),0.75));
% create triangle mesh grid
[V,F] = create_regular_grid(size(im,2),size(im,1),0,0);
% Scale (X,Y) to fit image
V(:,1) = V(:,1)*size(im,2);
V(:,2) = (1-V(:,2))*size(im,1);
V(:,3) = im(:);

Classic bump deformation example mesh

Saturday, June 28th, 2014

This is code I end up re-writing every once in a while to demonstrate k-harmonic deformation on a simple 2d domain. The domain is just a square. This codes meshes the domain so that two inner circles show up in the edge set. Then I can select the vertices on and outside the outer circle and on and inside the inner circle to create a bump.

R = 1;
r = 0.15;
outer_sam = 200;
inner_sam = ceil(outer_sam/R*r);
inner_theta = linspace(0,2*pi,inner_sam+1)';
outer_theta = linspace(0,2*pi,outer_sam+1)';
P_inner = r*[cos(inner_theta(2:end)) sin(inner_theta(2:end))];
P_outer = R*[cos(outer_theta(2:end)) sin(outer_theta(2:end))];
E_inner = [1:size(P_inner,1);2:size(P_inner,1) 1]';
E_outer = [1:size(P_outer,1);2:size(P_outer,1) 1]';
P_bb = 1.1*[-R -R;-R R;R R;R -R];
E_bb = [1 2;2 3;3 4;4 1];
[P_bb,E_bb] = resample_polygon(P_bb,E_bb,2*pi*R/outer_sam);
P = [P_inner;P_outer;P_bb];
E = [E_inner;size(P_inner,1)+[E_outer;size(P_outer,1)+E_bb]];
[V,F] = triangle(P,E,[],'Flags',sprintf('-Yq32a%0.17f',(2*pi*R/outer_sam)^2),'Quiet'); 

bump domain matlab

Compute and visualize self-intersections and intersections between two meshes

Saturday, April 26th, 2014

I add the self-intersection computation code from our winding numbers project to libigl. I made a small example (libigl/examples/intersections) that will compute and visualize all triangles participating in self-intersections in red. Here are the self-intersections of the Klein bottle:


I also wrote up a function to compute intersections between two meshes. Here’s that same example run on two meshes:

intersections between two meshes

Linearly interpolate colors in a triangle using SVG (Gouraud shading)

Monday, August 12th, 2013

After a long battle with the SVG documentation here is a triangle colored by linearly interpolated colors specified at each vertex.

<svg xmlns="" version="1.200000" width="100%" height="100%" viewBox="0 0 100.000000 86.600000" xmlns:xlink="">
  <g transform="matrix(1 0 0 -1 0 86.600000)">

      <linearGradient id="fadeA-1" gradientUnits="userSpaceOnUse" x1="50.000000" y1="0.000000" x2="50.000000" y2="86.600000">
        <stop offset="0%" stop-color="#FF0000"/>
        <stop offset="100%" stop-color="#000000" />
      <linearGradient id="fadeB-1" gradientUnits="userSpaceOnUse" x1="0.000000" y1="86.60000" x2="75.000000" y2="43.300000">
        <stop offset="0%" stop-color="#00FF00"/>
        <stop offset="100%" stop-color="#000000" />
      <linearGradient id="fadeC-1" gradientUnits="userSpaceOnUse" x1="100.000000" y1="86.60000" x2="25.000000" y2="43.300000">
        <stop offset="0%" stop-color="#0000FF"/>
        <stop offset="100%" stop-color="#000000" />

      <path id="pathA-1" d="M 50.000000,0.000000 L 0.000000,86.600000 100.000000,86.600000 Z" fill="url(#fadeA-1)"/>
      <path id="pathB-1" d="M 50.000000,0.000000 L 0.000000,86.600000 100.000000,86.600000 Z" fill="url(#fadeB-1)"/>
      <filter id="Default">
        <feImage xlink:href="#pathA-1" result="layerA" x="0" y="0" />
        <feImage xlink:href="#pathB-1" result="layerB" x="0" y="0" />
        <feComposite in="layerA" in2="layerB" operator="arithmetic" k1="0" k2="1.0" k3="1.0" k4="0" result="temp"/>
        <feComposite in="temp" in2="SourceGraphic"   operator="arithmetic" k1="0" k2="1.0" k3="1.0" k4="0"/>
    <g stroke="none" stroke-width="0" shape-rendering="crispEdges" >
      <path d="M 50.000000,0.000000 L 0.000000,86.600000 100.000000,86.600000 Z" fill="url(#fadeC-1)" filter="url(#Default)" />

SVG render

JPG rasterization

What you should see above if your browser supports SVG resterized linear interpolation of colors in triangle svg

Update: This still isn’t quite linear interpolation. For some reason I get a tendency toward black in the middle. I tried using fancier arithmetic (multiplication) with masks but this made things worse (quantization artifacts). There is also a nasty relative coordinates issue with the filters I used above.

Update: For reference, here’s what I’d like to see:

linear interpolation

Consistent (tet) mesh neighbor information

Thursday, January 10th, 2013

I’ve been using tet mesh from a matlab interface I wrote. It produces a list of tet indices, I call TT in matlab, where each row is a four-tuple of tet indices. It can also output a .neigh file which contains a list of each tet’s neighbor by index. I call this TN in matlab. This is also a list of rows of four, when now a positive index of j on row i means that tet i is a neighbor with tet j. Neighboring means sharing a face. Sharing a face means tet i and tet j have 3 indices in common.

One would expect that these come in a consistent ordering. For example, if j appears as the kth element on row i, then I might expect tet j to be the “kth neighbor” of tet i. So if k=3 I expect that tet i shares the face opposite its kth vertex (kth column in TT).

Doch leider ist es nicht wahr. (actually it is, see below).

Instead it seemed that tetgen was using a different ordering. To find out which consistent ordering scheme it was using I wrote a little script:

TE = [repmat(1:size(TT,1),1,size(TT,2))' TN(:)];
% verify edges are symmetric
A = sparse(TE(TE(:,2)>0,1),TE(TE(:,2)>0,2),1);

all_perms = perms(1:4);
% loop over all perms
for p = 1:size(all_perms,1)
  TI = [T(:,[2 3 4]); T(:,[3 4 1]); T(:,[4 1 2]); T(:,[1 2 3])];
  TI = [ ...
    T(:,all_perms(p,[2 3 4])); ...
    T(:,all_perms(p,[3 4 1])); ...
    T(:,all_perms(p,[4 1 2])); ...
    T(:,all_perms(p,[1 2 3]))];
  A = sparse(TE(TE(:,2)>0,1),TE(TE(:,2)>0,2),sum(TI(TE(:,2)>0),2));
  fprintf('%d: %d\n',p,full(max(max(abs(A-A')))));
  if max(max(abs(A-A'))) == 0
    % found it

First I do a little sanity check to make sure that tet’s agree that they are neighbors. Then I consider every (4! = 24) possible consistent ordering of tets by checking whether their shared faces agree.

But again! Foiled. Tetgen seems to output neighbors in an inconsistent manner. To resolve this I wrote another small script:

% Get rid of boundary edges
TE = TE(TE(:,2)>0,:);
TNI = zeros(size(TE,1),1);
for c = 1:size(TT,2)
  [TNIc] = all(~bsxfun(@eq,TT(TE(:,1),c),TT(TE(:,2),:)),2);
  % Assert face-manifoldness
  assert(all(TNI(TNIc) == 0));
  TNI(TNIc) = c;
% Assert that edges correspond to face sharers
% Rebuild TN
new_TN = full(sparse(TE(:,1),TNI,TE(:,2),size(TN,1),size(TN,2)));
new_TN(new_TN==0) = -1;

As a consequence the first script can be used to check for a consistent ordering and the second script can be used to check for face-manifoldness and correct neighbor relations.

Update: Turns out tetgen was fine and I had a bug. I was reordering the columns in my TT, without reordering the columns in TN. But these scripts did help me find and eliminate my bug!