Posts Tagged ‘triangle’

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);
assert(size(A,1)==size(A,2));
assert(max(max(A'-A))==0);

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
    break;
  end
end

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;
end
% Assert that edges correspond to face sharers
assert(all(TNI>0));
% 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!

CGAL intersections.h conflicts with Boolean_set_operations_2.h

Thursday, November 15th, 2012

I was trying to transform the cgal 2d polygon intersection example into an example that computes 3d triangle intersections.

This gave some typical pages of CGAL compiler errors:


In file included from /opt/local/include/CGAL/intersection_3.h:49:0,
                 from /opt/local/include/CGAL/Kernel/function_objects.h:34,
                 from /opt/local/include/CGAL/Cartesian/function_objects.h:28,
                 from /opt/local/include/CGAL/Cartesian/Cartesian_base.h:68,
                 from /opt/local/include/CGAL/Simple_cartesian.h:28,
                 from /opt/local/include/CGAL/Exact_predicates_exact_constructions_kernel.h:28,
                 from CGAL_includes.hpp:4:
/opt/local/include/CGAL/Triangle_3_Triangle_3_intersection.h: In instantiation of 'void CGAL::internal::intersection_coplanar_triangles_cutoff(const typename Kernel::Point_3&, const typename Kernel::Point_3&, const typename Kernel::Point_3&, const Kernel&, std::list<typename Kernel::Point_3>&) [with Kernel = CGAL::Simple_cartesian<CGAL::Gmpq>; typename Kernel::Point_3 = CGAL::Point_3<CGAL::Simple_cartesian<CGAL::Gmpq> >; typename R_::Point_3 = CGAL::Point_3<CGAL::Simple_cartesian<CGAL::Gmpq> >]':
/opt/local/include/CGAL/Triangle_3_Triangle_3_intersection.h:98:3:   Error: required from 'CGAL::Object CGAL::internal::intersection_coplanar_triangles(const typename Kernel::Triangle_3&, const typename Kernel::Triangle_3&, const Kernel&) [with Kernel = CGAL::Simple_cartesian<CGAL::Gmpq>; typename Kernel::Triangle_3 = CGAL::Triangle_3<CGAL::Simple_cartesian<CGAL::Gmpq> >]'
/opt/local/include/CGAL/Triangle_3_Triangle_3_intersection.h:136:53:   Error: required from 'CGAL::Object CGAL::internal::intersection(const typename Kernel::Triangle_3&, const typename Kernel::Triangle_3&, const Kernel&) [with Kernel = CGAL::Simple_cartesian<CGAL::Gmpq>; typename Kernel::Triangle_3 = CGAL::Triangle_3<CGAL::Simple_cartesian<CGAL::Gmpq> >]'
/opt/local/include/CGAL/Kernel/function_objects.h:2535:49:   Error: required from 'CGAL::CommonKernelFunctors::Intersect_3<K>::Object_3 CGAL::CommonKernelFunctors::Intersect_3<K>::operator()(const T1&, const T2&) const [with T1 = CGAL::Triangle_3<CGAL::Simple_cartesian<CGAL::Gmpq> >; T2 = CGAL::Triangle_3<CGAL::Simple_cartesian<CGAL::Gmpq> >; K = CGAL::Simple_cartesian<CGAL::Gmpq>; CGAL::CommonKernelFunctors::Intersect_3<K>::Object_3 = CGAL::Object]'
/opt/local/include/CGAL/Lazy.h:1756:51:   Error: required from 'CGAL::Lazy_construction_object<LK, AC, EC>::result_type CGAL::Lazy_construction_object<LK, AC, EC>::operator()(const L1&, const L2&) const [with L1 = CGAL::Triangle_3<CGAL::Epeck>; L2 = CGAL::Triangle_3<CGAL::Epeck>; LK = CGAL::Epeck; AC = CGAL::CommonKernelFunctors::Intersect_3<CGAL::Simple_cartesian<CGAL::Interval_nt<false> > >; EC = CGAL::CommonKernelFunctors::Intersect_3<CGAL::Simple_cartesian<CGAL::Gmpq> >; CGAL::Lazy_construction_object<LK, AC, EC>::result_type = CGAL::Object]'
/opt/local/include/CGAL/Triangle_3_Triangle_3_intersection.h:180:42:   Error: required from 'CGAL::Object CGAL::intersection(const CGAL::Triangle_3<R>&, const CGAL::Triangle_3<R>&) [with K = CGAL::Epeck]'
main.cpp:41:49:   Error: required from here
/opt/local/include/CGAL/Triangle_3_Triangle_3_intersection.h:62:64: Error: error: call of overloaded 'intersection(CGAL::CartesianKernelFunctors::Construct_line_3<CGAL::Simple_cartesian<CGAL::Gmpq> >::Line_3, CGAL::CartesianKernelFunctors::Construct_line_3<CGAL::Simple_cartesian<CGAL::Gmpq> >::Line_3, const CGAL::Simple_cartesian<CGAL::Gmpq>&)' is ambiguous
/opt/local/include/CGAL/Triangle_3_Triangle_3_intersection.h:62:64: note: candidates are:
In file included from /opt/local/include/CGAL/intersection_3_1.h:427:0,
                 from /opt/local/include/CGAL/intersection_3.h:29,
                 from /opt/local/include/CGAL/Kernel/function_objects.h:34,
                 from /opt/local/include/CGAL/Cartesian/function_objects.h:28,
                 from /opt/local/include/CGAL/Cartesian/Cartesian_base.h:68,
                 from /opt/local/include/CGAL/Simple_cartesian.h:28,
                 from /opt/local/include/CGAL/Exact_predicates_exact_constructions_kernel.h:28,
                 from CGAL_includes.hpp:4:
/opt/local/include/CGAL/Intersections_3/intersection_3_1_impl.h:200:1: Error: note: CGAL::Object CGAL::internal::intersection(const typename K::Line_3&, const typename K::Line_3&, const K&) [with K = CGAL::Simple_cartesian<CGAL::Gmpq>; typename K::Line_3 = CGAL::Line_3<CGAL::Simple_cartesian<CGAL::Gmpq> >]
In file included from CGAL_includes.hpp:11:0:
/opt/local/include/CGAL/Boolean_set_operations_2.h:904:23: Error: note: OutputIterator CGAL::intersection(InputIterator, InputIterator, OutputIterator, unsigned int) [with InputIterator = CGAL::Line_3<CGAL::Simple_cartesian<CGAL::Gmpq> >; OutputIterator = CGAL::Simple_cartesian<CGAL::Gmpq>]
In file included from /opt/local/include/CGAL/intersection_3.h:49:0,
                 from /opt/local/include/CGAL/Kernel/function_objects.h:34,
                 from /opt/local/include/CGAL/Cartesian/function_objects.h:28,
                 from /opt/local/include/CGAL/Cartesian/Cartesian_base.h:68,
                 from /opt/local/include/CGAL/Simple_cartesian.h:28,
                 from /opt/local/include/CGAL/Exact_predicates_exact_constructions_kernel.h:28,
                 from CGAL_includes.hpp:4:
/opt/local/include/CGAL/Triangle_3_Triangle_3_intersection.h: In instantiation of 'void CGAL::internal::intersection_coplanar_triangles_cutoff(const typename Kernel::Point_3&, const typename Kernel::Point_3&, const typename Kernel::Point_3&, const Kernel&, std::list<typename Kernel::Point_3>&) [with Kernel = CGAL::Simple_cartesian<CGAL::Interval_nt<false> >; typename Kernel::Point_3 = CGAL::Point_3<CGAL::Simple_cartesian<CGAL::Interval_nt<false> > >; typename R_::Point_3 = CGAL::Point_3<CGAL::Simple_cartesian<CGAL::Interval_nt<false> > >]':
/opt/local/include/CGAL/Triangle_3_Triangle_3_intersection.h:98:3:   Error: required from 'CGAL::Object CGAL::internal::intersection_coplanar_triangles(const typename Kernel::Triangle_3&, const typename Kernel::Triangle_3&, const Kernel&) [with Kernel = CGAL::Simple_cartesian<CGAL::Interval_nt<false> >; typename Kernel::Triangle_3 = CGAL::Triangle_3<CGAL::Simple_cartesian<CGAL::Interval_nt<false> > >]'
/opt/local/include/CGAL/Triangle_3_Triangle_3_intersection.h:136:53:   Error: required from 'CGAL::Object CGAL::internal::intersection(const typename Kernel::Triangle_3&, const typename Kernel::Triangle_3&, const Kernel&) [with Kernel = CGAL::Simple_cartesian<CGAL::Interval_nt<false> >; typename Kernel::Triangle_3 = CGAL::Triangle_3<CGAL::Simple_cartesian<CGAL::Interval_nt<false> > >]'
/opt/local/include/CGAL/Kernel/function_objects.h:2535:49:   Error: required from 'CGAL::CommonKernelFunctors::Intersect_3<K>::Object_3 CGAL::CommonKernelFunctors::Intersect_3<K>::operator()(const T1&, const T2&) const [with T1 = CGAL::Triangle_3<CGAL::Simple_cartesian<CGAL::Interval_nt<false> > >; T2 = CGAL::Triangle_3<CGAL::Simple_cartesian<CGAL::Interval_nt<false> > >; K = CGAL::Simple_cartesian<CGAL::Interval_nt<false> >; CGAL::CommonKernelFunctors::Intersect_3<K>::Object_3 = CGAL::Object]'
/opt/local/include/CGAL/Lazy.h:385:22:   Error: required from 'CGAL::Lazy_rep_2<AC, EC, E2A, L1, L2>::Lazy_rep_2(const AC&, const EC&, const L1&, const L2&) [with AC = CGAL::CommonKernelFunctors::Intersect_3<CGAL::Simple_cartesian<CGAL::Interval_nt<false> > >; EC = CGAL::CommonKernelFunctors::Intersect_3<CGAL::Simple_cartesian<CGAL::Gmpq> >; E2A = CGAL::Cartesian_converter<CGAL::Simple_cartesian<CGAL::Gmpq>, CGAL::Simple_cartesian<CGAL::Interval_nt<false> > >; L1 = CGAL::Triangle_3<CGAL::Epeck>; L2 = CGAL::Triangle_3<CGAL::Epeck>]'
/opt/local/include/CGAL/Lazy.h:1715:73:   Error: required from 'CGAL::Lazy_construction_object<LK, AC, EC>::result_type CGAL::Lazy_construction_object<LK, AC, EC>::operator()(const L1&, const L2&) const [with L1 = CGAL::Triangle_3<CGAL::Epeck>; L2 = CGAL::Triangle_3<CGAL::Epeck>; LK = CGAL::Epeck; AC = CGAL::CommonKernelFunctors::Intersect_3<CGAL::Simple_cartesian<CGAL::Interval_nt<false> > >; EC = CGAL::CommonKernelFunctors::Intersect_3<CGAL::Simple_cartesian<CGAL::Gmpq> >; CGAL::Lazy_construction_object<LK, AC, EC>::result_type = CGAL::Object]'
/opt/local/include/CGAL/Triangle_3_Triangle_3_intersection.h:180:42:   Error: required from 'CGAL::Object CGAL::intersection(const CGAL::Triangle_3<R>&, const CGAL::Triangle_3<R>&) [with K = CGAL::Epeck]'
main.cpp:41:49:   Error: required from here
/opt/local/include/CGAL/Triangle_3_Triangle_3_intersection.h:62:64: Error: error: call of overloaded 'intersection(CGAL::CartesianKernelFunctors::Construct_line_3<CGAL::Simple_cartesian<CGAL::Interval_nt<false> > >::Line_3, CGAL::CartesianKernelFunctors::Construct_line_3<CGAL::Simple_cartesian<CGAL::Interval_nt<false> > >::Line_3, const CGAL::Simple_cartesian<CGAL::Interval_nt<false> >&)' is ambiguous
/opt/local/include/CGAL/Triangle_3_Triangle_3_intersection.h:62:64: note: candidates are:
In file included from /opt/local/include/CGAL/intersection_3_1.h:427:0,
                 from /opt/local/include/CGAL/intersection_3.h:29,
                 from /opt/local/include/CGAL/Kernel/function_objects.h:34,
                 from /opt/local/include/CGAL/Cartesian/function_objects.h:28,
                 from /opt/local/include/CGAL/Cartesian/Cartesian_base.h:68,
                 from /opt/local/include/CGAL/Simple_cartesian.h:28,
                 from /opt/local/include/CGAL/Exact_predicates_exact_constructions_kernel.h:28,
                 from CGAL_includes.hpp:4:
/opt/local/include/CGAL/Intersections_3/intersection_3_1_impl.h:200:1: Error: note: CGAL::Object CGAL::internal::intersection(const typename K::Line_3&, const typename K::Line_3&, const K&) [with K = CGAL::Simple_cartesian<CGAL::Interval_nt<false> >; typename K::Line_3 = CGAL::Line_3<CGAL::Simple_cartesian<CGAL::Interval_nt<false> > >]
In file included from CGAL_includes.hpp:11:0:
/opt/local/include/CGAL/Boolean_set_operations_2.h:904:23: Error: note: OutputIterator CGAL::intersection(InputIterator, InputIterator, OutputIterator, unsigned int) [with InputIterator = CGAL::Line_3<CGAL::Simple_cartesian<CGAL::Interval_nt<false> > >; OutputIterator = CGAL::Simple_cartesian<CGAL::Interval_nt<false> >]
make: *** [obj/main.o] Error 1

I fixed this by making sure in the files where I include:


#include <CGAL/intersections.h>

I don’t include:


#include <CGAL/Boolean_set_operations_2.h>

Triangle wave for indices (integers)

Friday, May 4th, 2012

Often I have N things and I want to assign them to n<<N groups/colors/positions etc. If I don’t care too much about their order or distribution I might randomly assign them. If I care that sequential items go into different groups and that groups get the same number of items I might use:


ids = mod(indices,n)+1

Which plotted as a function of the indices looks like a sawtooth wave:
index mod sawtooth

Another way to get non-repeating numbers, with near-equal distribution is to use a triangle wave, counting up to n and then back down to 1 and repeating:


ids = abs(mod(indices,n*2-2)-n+1)+1

which looks like:
index mod triangle

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;
  end
  figure(f);

  % 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
  view(2);

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

  % loop until mouse up or ESC is pressed
  done = false;
  while(~done)
    drawnow;
  end

  % 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);
    append_current_point();
  end

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

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

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

  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;
    else
      % update plot
      set(p,'Xdata',P(:,1),'Ydata',P(:,2));
    end
  end

  function finish()
    done = true;
    set(gcf,'windowbuttonmotionfcn','');
    set(gcf,'windowbuttonupfcn','');
    set(gcf,'windowbuttondownfcn','');
    set(gcf,'keypressfcn','');
  end

end

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

intrinsic cotangent formula

Monday, January 23rd, 2012

Here’s a derivation of the cotangent of an angle α of a triangle ABC given just the lengths of each side of the triangle.
acute triangle

First of all we remind ourselves what cotangent is in terms of cosine and sine:


        cos α
cot α = ------
        sin α

Now look at cosine and sine separately. By the law of cosines we have that


a2 = b2 + c2 - 2bc cos α

Rearranging things we have that


        -a2 + b2 + c2
cos α = ------------
             2bc

Now looking at the sine, we start with the familiar area of the triangle treating b as the base.


     1
A = --- bh
     2

Using SOH-CAH-TOA, we replace the height:


     1
A = --- bc sin α
     2

Rearranging this we have:


         2A
sin α = ----
         bc

Now put the cosine and sine derivations together to get:


        cos α    -a2 + b2 + c2   bc
cot α = ------ = ------------  ----
        sin α        2bc        2A

Finally arriving at:


        -a2 + b2 + c2
cot α = ------------
             4A

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.

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.

Convert quad mesh OBJ to triangle mesh OBJ using regular expressions (with bash and sed)

Wednesday, December 15th, 2010

Here’s a simple bash one-liner that converts a quad mesh stored in an OBJ file into a triangle mesh stored in a OBJ file. It splits every quad (a,b,c,d) into two triangles (a,b,c) and (a,c,d).

I issue this in the bash shell of my mac (you may have to tweak the sed syntax):


cat quads.obj | sed -E -e "s/f ([0-9\/]+) ([0-9\/]+) ([0-9\/]+) ([0-9\/]+)/f \1 \2 \3`echo -e "\r"`f \1 \3 \4/g" > triangles.obj

Or I open the file with vim and use:


:%s/f \([0-9\/]\+\) \([0-9\/]\+\) \([0-9\/]\+\) \([0-9\/]\+\)/f \1 \2 \3\rf \1 \
3 \4/g

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