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);
NA-NB
NB-NC
producing
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( ...
NABC(sub2ind(size(NABC),repmat(1:size(NABC,1),size(NABC,2),1)',I)),size(NABC));
N = 1/3*reshape(sum(sNABC,2),shape);