# How to calculate the average direction of 3D elements using normal vectors of their contour surface?

I have a geometry which consists of a number of slender cylindrical elements. Each cylinder is described by a lot of little triangular planes on the surface of the cylinder and I know the normal vector of each triangular plane.

I am able to calculate the orientation of the (axis of the) cylinder if there is only 1 cylinder in the file. (solved in How to calculate the direction in which a set of normal vectors (3D) are least oriented?) Now I am looking for the best way to calculate the average orientation of the cylinders if there are many cylinders in the file and when I don’t know which triangle belongs to which cylinder. If there are many cylinders with a preferential orientation, I can calculate the vector in which the cylinders are more or less aligned, but how do I calculate the average angle between this vector and the cylinders?

Is there a way to calculate the average angle of the cylinders with x-, y- and z-axis of a Cartesian coordinate system? As all cylinders have the same diameter and length, I think/hope this must be possible. I am looking for the average orientation of the cylinders (not of the normal vectors) with another vector.

Simple example: If I want to calculate the angles in the particular case of a slender cylinder along the z-axis, the average angle with the z-axis will be 0. But in the incorrect way I calculate it now, the average angle with x- and y-axis is pi/4 instead of pi/2 because the normal vectors are perpendicular to the z-axis but this can be in any direction in the x-y-plane. So, again, I am looking for a better algorithm.

-
joriki paraphrased my question very well! I want to go from a relative sense (my old version) to an absolute sense (absolute angles instead of just an indicator), so I probably should have used ‘direction’ in the title instead of ‘orientation’. – Joren Oct 24 '12 at 10:32
There's an edit link under the question; you can update the question any time you like with clarifications or additional information. In fact I think you should if you think you can now express the question more clearly. – joriki Oct 24 '12 at 11:11

I find your use of the term "orientation" slightly confusing, since it seems to oscillate between an absolute sense (roughly synonymous with "direction") and a relative sense (roughly synonymous with "alignment"). Let me paraphrase how I understand your question:

You have triangles for more than one cylinder in a file, and you don't know which triangle belongs to which cylinder. The normals of the triangles are all roughly orthogonal to the axis of their cylinder. For a given direction, you want to know how well the axes of the cylinders are aligned with this direction. In an extreme case, the axes of all cylinders might have the same direction; in that case, the test with that direction should yield "fully aligned" and a test with any direction orthogonal to it should yield "fully non-aligned". The problem is that though the cylinders may be fully aligned, the normals of the triangles form different angles with any given direction.

It seems to me the "incorrect way" that you're using isn't quite so bad. I don't fully understand how you arrive at $0$ and $\pi/4$, respectively, but an indicator that ranges between $0$ for aligned and $\pi/4$ for non-aligned seems useful. I gather that what you don't like about it is that you want it to be some sort of average angle.

You could probably find a function that maps your current indicator to something like an average angle, but generally speaking, averaging trigonometric functions tends to be easier and to have nicer properties than averaging angles (or absolute values of angles, which seems to be what you're doing), so I'll do something similar to what you did, but using trigonometric functions whose averages we can easily calculate exactly, so we'll get an expression in closed form relating the values of the indicator to something like an average angle.

For this to work, the triangles for each cylinder need to be roughly equidistributed about its axis. If the normals of the triangles tend to preferentially point in one of the directions orthogonal to the axis, that will make the cylinder appear more aligned with the third direction (orthogonal to the axis and the preferred normal direction) than it actually is.

So consider a cylinder with unit vector $\vec a$ along its axis, a unit vector $\vec u$ along which we want to measure the cylinder's alignment, and unit vectors $\vec n_1$ and $\vec n_2$ orthogonal to $\vec a$, which we choose such that $\vec n_1$ is also orthogonal to $\vec u$ (which we can always do). We can write the normal of a triangle of the cylinder as $\vec n=\cos\phi\,\vec n_1+\sin\phi\,\vec n_2$, and equidistribution of the normals around the cylinder's axis amounts to equidistribution of $\phi$ over $[0,2\pi]$. We have $\vec u\cdot\vec n_1=0$ (since we chose $\vec n_1$ orthogonal to $\vec u$) and $\vec u\cdot\vec n_2=\pm\sin\theta$, were $\theta$ is the angle between $\vec a$ and $\vec u$ (you may want to draw a sketch to see why that's the case). Thus $\vec u\cdot\vec n=\pm\sin\theta\sin\phi$. Now if we average this over all the equidistributed normals, we'll just get $0$ from the average over $\phi$, independent of $\theta$; that's not what we want. However, if we square before averaging, we get the average of $\sin^2\theta\sin^2\phi$ over $\phi$, which is $\frac12\sin^2\theta$. That's good, because we shouldn't be distinguishing between negative and positive values of $\theta$ or between $\theta$ and $\pi-\theta$ anyway.

Now to turn this into something like an average angle, we just have to apply the inverse operations. That is, if $\langle(\vec u\cdot\vec n)^2\rangle$ denotes the average of the square of the scalar product of the alignment direction $\vec u$ with the triangle normals $\vec n$, the "average angle" $\hat\theta$ would be derived from $\langle(\vec u\cdot\vec n)^2\rangle=\langle\frac12\sin^2\theta\rangle$ as

$$\hat\theta=\arcsin\sqrt{2\langle(\vec u\cdot\vec n)^2\rangle}\;.$$

You'll need to clamp any arguments above $1$ for the arcsine to $1$; this might occur if the normals aren't perfectly equidistributed.

If all the cylinders are perfectly aligned along $\vec u$, then all the normals will be orthogonal to $\vec u$, so $\langle(\vec u\cdot\vec n)^2\rangle$ will be zero, and so will $\hat\theta$. On the other hand, if all the cylinders are perfectly aligned orthogonal to $\vec u$ and their normals are equidistributed about their axes, then $\langle(\vec u\cdot\vec n)^2\rangle$ will be $1/2$, so $\hat\theta$ will be $\pi/2$ as desired. More generally, if all the cylinder axes form an angle $\theta$ with $\vec u$, then $\hat\theta$ will be $\theta$. If the cylinders form different angles with $\vec u$, then $\hat\theta$ will be an average angle in the sense that $\sin^2\theta$ (or equivalently $\cos^2\theta$) is averaged.

-
The normals of the triangles are more or less equidistributed around the axis of the cylinder, but not every cylinder in the file has the same direction. Looks like the angle between u and n should be larger than 45° to calculate the arcsin(). This is not always the case (some cylinders have a different direction). Why do I have to clamp ^theta to the value $pi/2$? I think this method may work for directions u which are more or less equal to the main direction of the cylinders, but for other directions u the arcsin will lead too often to $pi/2$ if I calculate it this way. What is your opinion? – Joren Oct 24 '12 at 10:37
@Joren: I think that's a fundamental misunderstanding of the approach. I'm aware that not all cylinders have the same direction. The arcsine is calculated not from a single dot product $\vec u\cdot\vec n$ but from the average of the squares of all the dot products. The answer explains why this averages to at most $1/2$. Please read the answer again carefully and perhaps try to implement it to see how it works. – joriki Oct 24 '12 at 10:57
@Joren: To answer your question: If the cylinders are perfectly aligned orthogonal to the test direction and the normals aren't perfectly equidistributed around their axes but tend to point in the test direction, then $\sin^2\phi$ may average to slightly more than $1/2$ and the argument of the arcsine may be greater than $1$. Since this is an artefact of the failure of equidistribution, as long as approximate equidistribution can be assumed it still shows that the cylinders are almost perfectly orthogonal to the test direction, so it makes sense to clamp the argument to $1$. – joriki Oct 24 '12 at 11:05
@Joren: By the way, note that your current indicator also relies on the equidistribution of the normals. As I wrote, you could take your average value of $\pi/4$ to indicate "non-aligned", but you only get that average if the normals are equidistributed. – joriki Oct 24 '12 at 11:14
Thank you for the clarification. Taking the average of the squares of the scalar product before taking the arcsin seems to work! I tested the code for some small test cases and the result is in accordance with the real direction. Thanks again! – Joren Oct 24 '12 at 11:49