triangulate_pareByArea(XYZ, v, maxTriangleArea)

% v = triangulate_pareByArea(XYZ, v, maxTriangleArea)
%
% pare results of triangulate by area of triangles
% generated. this is useful because if the original XYZ were
% on a regular grid, erroneous triangles completing the
% convex shape enclosing all the XYZ will be larger than the
% typical area, and can be removed.
%
% note, only the area of the projection to the XY plane is
% used.

function v = triangulate_pareByArea(XYZ, v, maxTriangleArea)

N = size(v,2);

tris = cat(3, XYZ(:,v(1,:)), XYZ(:,v(2,:)), XYZ(:,v(3,:)));
tris = permute(tris, [1 3 2]);

accepted = zeros(N,1);

for n=1:N
	tri = tris(:,:,n);
	A = tri(:,1) - tri(:,2);
	B = tri(:,2) - tri(:,3);
	C = tri(:,3) - tri(:,1);
	a = sqrt(sum(A.^2));
	b = sqrt(sum(B.^2));
	c = sqrt(sum(C.^2));
	area = 0.25 * sqrt((a+b+c)*(b+c-a)*(c+a-b)*(a+b-c))
	if area <= maxTriangleArea accepted(n) = 1; end
end

v = v(:, find(accepted));