Posts Tagged ‘circle’

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

Re: Suggester shapes

Sunday, May 8th, 2011

Ken recently posted about shapes that suggest math problems. At the end of his post he leaves a riddle:

Find a shape whose area is 1/2 the area of a square and whose perimeter is the same as the perimeter of the square.

Thinking about this riddle a little bit got me thinking about other similar riddles. Like what if I ask the same question but instead of a square use a circle.

Find a shape whose area is 1/2 the area of a circle and whose perimeter is the same as the perimeter of the circle.

My idea for a solution for this new riddle was to intersect to equal-size circles. Notice as long as the circles are equal size if I take the shape that’s one of the circles minus the other circle, the perimeter will be the same as a full circle.

So then I just need to find how far apart the circles should be. I know that if the distance between the circles’ centers is zero then the left over shape will have area zero. And if the distance is twice their radii then then the left over area will be equal to a full circle. So if I want the left over area to be half the area of a circle I need the distance to be somewhere between zero and twice the radii.

I wrote out the algebra blindly hoping that the equations would collapse leaving me with a beautiful, succinct expression for the necessary distance. The problem turns out to be fairly intangible. Leaving me with a gnarly expression for the distance that is not easily solved. Here’s a picture and the expression:

shape with half area but full perimeter of a circle

I can approximate d numerically, which reveals d ≈ 0.8079455066. I’m not the first person to try this, this problem also makes an appearance on the Wolfram circle-circle intersection page.

This leaves me wondering if there is a simpler answer to my riddle. Maybe the solution doesn’t involve circle-circle intersection at all!

Barycentric versus voronoi regional area in mass matrices

Friday, June 25th, 2010

In a recent project, we were minimizing an energy over a surface approximated by a triangle mesh. Our finite element method approach reduced our energy to solving a system of linear equations with coefficients from a stiffness matrix and mass matrix. In the lumped version of a mesh’s mass matrix, diagonal entries correspond to areas on the surface around each vertex. You want all the areas to sum to the total area of the surface and you would like the area for each vertex to some how reflect the influence of its incident triangles.


One natural choice is a barycentric area. Take the barycenter (centroid, center of mass), which lies at the intersection of medians. Using those medians divide the triangle into three quadrilaterals each corresponding its incident point:

generic barycentric quads

These quadrilaterals are by definition of equal area. The barycentric mass matrix entry for vertex i is the sum of its corresponding quadrilaterals from all incident triangles, or more simply: one-third of the total area of all neighboring triangles. Notice that the contribution of a triangle to its vertices is the same for all three, regardless of there placement around the triangle.

generic barycentric region on mesh


Another, also natural, choice is voronoi area. Take the circumcenter, which lies at the center of the circumscribing circle of the triangle. It follows quickly that the circumcenter is equidistant from all vertices on the triangle. Connect this point to the midpoints along the edges of the triangle, again forming three quadrilaterals.

acute voronoi quads

Note that this time the quadrilaterals are not equal. They are not even guaranteed to be positive. For obtuse triangles, the acute angled vertices will receive negative areas (assuming you demand their sums equal the area of the triangle). Here is a plot of voronoi areas as the angle of an isosceles triangle grows from 0 (acute) to π (obtuse). The x axis is the angle of the isosceles triangle and the y axis is area. The green line shows the area of corresponding to the vertex at the angle that’s growing, the blue corresponds to the other two areas.

voronoi isoceles plot

Notice that at π/2 we might expect the other two areas to cross into negative but actually this happens later. So voronoi areas won’t produce negative coefficients until fairly obtuse.

Again, the voronoi mass matrix entry for vertex i is the sum of its corresponding quadrilaterals from all incident triangles. Since areas can be negative, this means negative entries could end up in the voronoi mass matrix.

generic voronoi on mesh

Connectivity dependence

The Voronoi area of a vertex is called so because the area represents its corresponding voronoi region: all the points on the surface for which this vertex is the closest vertex (assuming the mesh has “nice” triangles: non-obtuse). In this way the voronoi area depends only on the vertices positioning on the surface not the connectivity of the mesh.

The barycentric area, while having the advantage of always being positive, is heavily connectivity dependent as a vertex gets one-third of the area of any connected triangle.

For a simple example, here’re two meshes. One with regular positions and regular connectivity:
regular positions, regular connectivity

And one with regular positions and the same regular connectivity but a single edge at each highlighted vertex has been flipped:
regular positions, one edge flipped

And here are the barycentric (left) and voronoi (right) areas for the highlighted vertices on the regular connectivity mesh.
regular with areas

Now, notice how when the edge is flipped the barycentric area changes drastically, and actually changes for all vertices incident on the flipped edge. The voronoi however stays the same for this vertex and all vertices.
one-edge flipped with areas

Hybrid approach to vonoroi

To handle the problem of negative voronoi areas for some obtuse triangles, there exists a hybrid approach. For acute triangles, use voronoi as is. For obtuse triangles snap the circumcenter to the midpoint along the edge opposite the large angle and use the corresponding quads (I described this in an earlier post). Notice for a right angle triangle the circumcenter corresponds with this midpoint, so as a triangle becomes obtuse there isn’t a jump in area distribution. A keen eye will notice that I am doing this in the voronoi area picture above for the obtuse triangle.

Areas of quadrilaterals within a triangle determined by its circumcenter (or perpendicular bisectors), using matlab

Thursday, February 11th, 2010

First, I need to locate the circumcenter of the triangle. The circumcenter is the center of the circle circumscribing the three points of the triangle. Thus, it is equidistance from all three points.
circumscribed triangle circumcenter
I do not need the actual world coordinates of the circumcenter, rather it is okay for me to just know the relative distance to each corner. The law of cosines tells me that the relative distance of the circumcenter from the line opposite point A (called a in the picture) is cos A (that is cosine of the angle at A). Likewise for B and C. Knowing these I can convert to barycentric coordinates by multiplying by the opposite edge lengths.

cosines = [ ...
  (edge_norms(:,3).^2+edge_norms(:,2).^2-edge_norms(:,1).^2)./(2*edge_norms(:,2).*edge_norms(:,3)), ...
  (edge_norms(:,1).^2+edge_norms(:,3).^2-edge_norms(:,2).^2)./(2*edge_norms(:,1).*edge_norms(:,3)), ...
barycentric = cosines.*edge_norms;
normalized_barycentric = barycentric./[sum(barycentric')' sum(barycentric')' sum(barycentric')'];

A quality of the barycentric trilinear coordinates of a point in a triangle is that they specify the relative area of inner triangles form by connecting the corners to that point. In our case, the triangles made by connecting the circumcenter to each corner. Here, the aquamarine, red, and purple triangles corresponding to points A, B, and C.
circumscribed triangle partial triangles
Now I find the area of the triangle and then multiple against the normalized barycentric coordinates of the circumcenter to get the areas of each of these inner triangles:

areas = 0.25*sqrt( ...
  (edge_norms(:,1) + edge_norms(:,2) - edge_norms(:,3)).* ...
  (edge_norms(:,1) - edge_norms(:,2) + edge_norms(:,3)).* ...
  (-edge_norms(:,1) + edge_norms(:,2) + edge_norms(:,3)).* ...
  (edge_norms(:,1) + edge_norms(:,2) + edge_norms(:,3)));
partial_triangle_areas = normalized_barycentric.*[areas areas areas];

Finally the quadrilateral at each point is built from half of the neighboring inner triangles: the green, blue and yellow quads corresponding to A, B, and C.
circumscribed triangle quadrilaterals
Note, that the fact that these inner triangles are cut perfectly in half follows immediately from the fact the the circumcenter is equidistance from each corner. So the last step is just to add up the halves:

quads = [ (partial_triangle_areas(:,2)+ partial_triangle_areas(:,3))*0.5 ...
  (partial_triangle_areas(:,1)+ partial_triangle_areas(:,3))*0.5 ...
  (partial_triangle_areas(:,1)+ partial_triangle_areas(:,2))*0.5];

Note: I’ve suffered trying to use mathematic or maple to reduce these steps to a simple algebra formula (I must be mistyping something). As far as I can tell the above steps check out (of course, it only makes sense for cases where the circumcenter is within or on the triangle: i.e. the triangle is not obtuse). I have checked out my equations on maple with the following functions:


Which then you can run:


To see a beautiful long algebra formula for the quadrilateral at point A.
Which in matlab corresponds to:

quads = ...
  [-sqrt(-(edge_norms(:,1) + edge_norms(:,2) + edge_norms(:,3)).* ...
    (edge_norms(:,1) - edge_norms(:,2) - edge_norms(:,3)).* ...
    (edge_norms(:,1) - edge_norms(:,2) + edge_norms(:,3)).* ...
    (edge_norms(:,1) + edge_norms(:,2) - edge_norms(:,3))).* ...
    (2.*edge_norms(:,2).^2.*edge_norms(:,3).^2  + edge_norms(:,1).^2.* ...
    edge_norms(:,2).^2  - edge_norms(:,2).^4  + edge_norms(:,1).^2.* ...
    edge_norms(:,3).^2  - edge_norms(:,3).^4)./ ...
    (8.*(edge_norms(:,1) + edge_norms(:,2) + edge_norms(:,3)).* ...
    (edge_norms(:,1) - edge_norms(:,2) - edge_norms(:,3)).* ...
    (edge_norms(:,1) - edge_norms(:,2) + edge_norms(:,3)).* ...
    (edge_norms(:,1) + edge_norms(:,2) - edge_norms(:,3))), ...
  -sqrt(-(edge_norms(:,2) + edge_norms(:,3) + edge_norms(:,1)).* ...
    (edge_norms(:,2) - edge_norms(:,3) - edge_norms(:,1)).* ...
    (edge_norms(:,2) - edge_norms(:,3) + edge_norms(:,1)).* ...
    (edge_norms(:,2) + edge_norms(:,3) - edge_norms(:,1))).* ...
    (2.*edge_norms(:,3).^2.*edge_norms(:,1).^2  + edge_norms(:,2).^2.* ...
    edge_norms(:,3).^2  - edge_norms(:,3).^4  + edge_norms(:,2).^2.* ...
    edge_norms(:,1).^2  - edge_norms(:,1).^4)./ ...
    (8.*(edge_norms(:,2) + edge_norms(:,3) + edge_norms(:,1)).* ...
    (edge_norms(:,2) - edge_norms(:,3) - edge_norms(:,1)).* ...
    (edge_norms(:,2) - edge_norms(:,3) + edge_norms(:,1)).* ...
    (edge_norms(:,2) + edge_norms(:,3) - edge_norms(:,1))), ...
  -sqrt(-(edge_norms(:,3) + edge_norms(:,1) + edge_norms(:,2)).* ...
    (edge_norms(:,3) - edge_norms(:,1) - edge_norms(:,2)).* ...
    (edge_norms(:,3) - edge_norms(:,1) + edge_norms(:,2)).* ...
    (edge_norms(:,3) + edge_norms(:,1) - edge_norms(:,2))).* ...
    (2.*edge_norms(:,1).^2.*edge_norms(:,2).^2  + edge_norms(:,3).^2.* ...
    edge_norms(:,1).^2  - edge_norms(:,1).^4  + edge_norms(:,3).^2.* ...
    edge_norms(:,2).^2  - edge_norms(:,2).^4)./ ...
    (8.*(edge_norms(:,3) + edge_norms(:,1) + edge_norms(:,2)).* ...
    (edge_norms(:,3) - edge_norms(:,1) - edge_norms(:,2)).* ...
    (edge_norms(:,3) - edge_norms(:,1) + edge_norms(:,2)).* ...
    (edge_norms(:,3) + edge_norms(:,1) - edge_norms(:,2)))];

Note: It seems also that the wikipedia page for the circumscribed circle has a direct formula for the barycentric coordinates of the circumcenter. Though, wolfram’s mathworld has an ostensibly different formula. Upon a quick review Wolfram’s is at least easily verifiable to be the same as mine (recall these only need to be relative to each other).

Average ratio of incenters to circumcenters, using matlab

Monday, February 8th, 2010

The incenter radius of a triangle is given by the formula:

      2*A      with A = area
r = -------

the circumcenter radius is given by

R = -----

We want the ratio of the incenter to circumcenter so that’s

r         8A2
_ = ---------------
R   ((a+b+c)*(abc))

So then if you have a triangle mesh with a variable edge_norms which is an n by 3 matrix of the lengths of edges in your mesh, a way to get the incenter radii is, first applying Heron’s formula to find the areas of each triangle:

semi_perimeters = (edge_norms(:,1) + edge_norms(:,2) + edge_norms(:,3))*0.5;
areas = sqrt( ...
  semi_perimeters.* ...
  (semi_perimeters-edge_norms(:,1)).* ...
  (semi_perimeters-edge_norms(:,2)).* ...
ratios = 8*areas.*areas./(sum(edge_norms')'.*prod(edge_norms')')

Then just take the mean:


Note: I coded up the calculations of the incenters and circumcenters to check this out so I figure I might as well post that, too: