If you do not care about speed too much you can sort the neighboring points by ang = atan2(a,b,n)
angle:
n - view direction
ang - CW? angle
sin(ang) = ((n x a).b)/(|a|*|b|)
cos(ang) = (a.b) /(|a|*|b|)
if (sin(ang)>=0) ang= acos(cos(ang))
if (sin(ang)< 0) ang=2*Pi - acos(cos(ang))
ang = <0,2*Pi)
So if you got point p0
and its neighbors p1,p2,p3,...
then for example set view vector as:
n = (p1-p0) x (p2-p1)
n = n / |n|
and then compute angle for each point
a1 = 0
a2 = atan2(p1-p0,p2-p0,n)
a3 = atan2(p1-p0,p3-p0,n)
...
then simply sort the points p1,p2,... by their angles and then you can simply triangulate (as triangle fan using p0 as base point of the fan)...
triangle(q0,q1,q2)
triangle(q0,q2,q3)
triangle(q0,q3,q4)
...
where q0,q1,q2,q3...
are the points p0,p1,p2,p3,...
in sorted order
From this you simply can apply your normal computation ... Also check this:
In case you want to weight the normals by triangle size then simply do not normalize the partial normals and simply normalize the resulting sum normal instead of dividing it by cnt
...
btw this can be done also without any goniometrics ... you can use bubble sort and simply sort they points so the comply selected winding rule so for example:
if (dot(n,cross(p(i)-p0,p(j)-p0)) > 0) swap(p(i),p(j))
However this will prevent the use of faster sorting algorithms so unless some spatial subdivision of points is used it will be slow for too many points... with spatial subdivision or small amount of point this should be fast. Also this approach does not require any additional memory space as no angle is required and sorting can be done in-place.