Posts Tagged ‘area’

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

and

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)

produces

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).

Signed polygon area in matlab

Sunday, October 30th, 2016

Suppose you have a polygon’s 2D corners stored in P, then the signed area is given by:

signed_polyarea = @(P) 0.5*(P(1:end,1)'*P([2:end 1],2)-P(1:end,2)'*P([2:end 1],1));

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)), ...
  (edge_norms(:,1).^2+edge_norms(:,2).^2-edge_norms(:,3).^2)./(2*edge_norms(:,1).*edge_norms(:,2))];
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:


mycos:=(a,b,c)->(b^2+c^2-a^2)/(2*b*c);
barycentric:=(a,b,c)->a*mycos(a,b,c);
normalized_barycentric:=(a,b,c)->barycentric(a,b,c)/(barycentric(a,b,c)+barycentric(b,c,a)+barycentric(c,a,b));
area:=(a,b,c)->sqrt((a+b-c)*(a-b+c)*(-a+b+c)*(a+b+c))/4;                                            
inner_triangle:=(a,b,c)->normalized_barycentric(a,b,c)*area(a,b,c);
quad:=(a,b,c)->(inner_triangle(b,c,a)+inner_triangle(c,b,a))/2;

Which then you can run:


simplify(quad(a,b,c));

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).