How do I calculate the normals of a point group?

Question

I hope to come up with an interesting and tricky task this time, and would be happy to receive any expert solutions.

First of all: I'm looking for a solution that uses only simple vector math!

The exercise is as follows:

The easy part

I have three groups of points for which I want to find the common normals. The direction, whether upwards or downwards, doesn't matter. The main thing is that these normals point in the (positive or negative) direction in which the face connecting the points should point. The face from which these points originate has been rotated in different directions and no longer exists.

The points of a triangle, a square and a pentagon serve as a simple example here (here connected to a face and numbered for better illustration):

Now I want to find the common normal of this non-existent face for each point or point group (illustrated here with an arrow):

The hard part

The funny thing about this task, however, and this is where it gets tricky: the only information I have are the points themselves and the index of the point group. Nothing else.

Which enthusiastic Geometry Nodes expert has a helpful answer to this? All ideas are welcome, but as mentioned before, I am explicitly looking for a solution that is purely mathematical or based on vector math.

Here is a starting point in blend format:

_{(Blender 4.0.2+)}

---

Update/Bonustask (?)

Thanks to Nathan, the comments show that this exercise turns out to be unsolvable if the points are not exactly on the same plane (which may be the case). ...?

In this case, this exercise can be supplemented with information about which face center the respective points were connected to.

The following picture illustrates this:

Perhaps this task can be solved with this additional information.
And I'll add that in this case a detour via curves is also an option, if necessary, but I actually wanted to avoid that here.

_{(Blender 4.0.2+)}

Answer 1

(Using Blender 3.6.5)

Documentation:

Preliminary remarks:
As very well explained by Mr A in his proposal, Principal Component Analysis is a well adapted framework for this class of problems. One part of the burden of the general theory he recapped can be alleviate in the present case thanks to the coplanar nature of the points in a same group. Furthermore because only one vector, the normal, is sought, the problem can even be reduced to solving a homogeneous linear equation.
However in the following, classical formulations from this mathematical field are not used because Blender toolbox lack of linear algebra components, mostly matrices. Blender vectors are only geometrical, restricted to 3D space. The assumed choice is to match the algebra available in Geometry Nodes, even if the formulation might seem out of phase for experts. This way graph implementation is following closely mathematical developments.

Mathematical formulation:
Only the points in a same group, sharing the same Group Index and the same normal, are considered. Managing more than one group is detailed in the implementation.
Let $n$ be the number of points, $M_{i, 1\leq i \leq n}$ the points, $C$ the face centre the points were connected to, $\vec{n}$ the sought normal. If $C$ were not available, it would have be replaced by the points barycentre as in lemon proposal.
$C$ and $\vec{n}$ are defining the plane the points lay on, because $\vec{CM_i}$ is perpendicular to $\vec{n}$. Using the dot product "$\cdot$", it yields for every point $i$ : $$\vec{CM_i} \cdot \vec{n} = 0 \label{CMdorN} \tag{1}$$ All these scalar equations are collected in a matrix form, with matrix $Q$ assembled from each vector $\vec{CM_i}$ written as a line (i.e. the transpose of a column vector) and vector $N$ assembled from the coordinates of $\vec{n}$.
$Q$ is defined as: $$ Q = \left[ {\begin{array}{c} \vec{CM_1} \\ \vec{CM_2} \\ \vdots \\ \vec{CM_n} \end{array}} \right] = \left[ {\begin{array}{ccc} x_1 & y_1 & z_1 \\ x_2 & y_2 & z_2 \\ \vdots & \vdots & \vdots \\ x_n & y_n & z_n \end{array}} \right] \tag{2}$$ where $x_i, y_i, z_i$ are the coordinates of $\vec{CM_i}$.
$N$ is defined as: $$N = \left[ {\begin{array}{c} n_x \\ n_y \\ n_z \end{array}} \right] \tag{3}$$ Equation ($\ref{CMdorN}$) is then rewritten as: $$ Q N = 0_n \label{QNeq0} \tag{4}$$ where $0_n$ is the column vector of size $n$ full of zeros.
As $Q$ is rectangular, and not square, the "best" solution $N$ that could be achieved is the one minimizing the norm of $Q N$. It is to notice that if some points $M_i$ are out of the plane, the right-hand side of Equation ($\ref{QNeq0}$) is not equal to zero. But if it is small, the minimization of the norm of $Q N$ returns a satisfactory "mean" normal.
Let $Q^T$ be the transpose of $Q$, such that: $$ Q^T = \left[ {\begin{array}{cccc} x_1 & x_2 & \ldots & x_n \\ y_1 & y_2 & \ldots & y_n \\ z_1 & z_2 & \ldots & z_n \end{array}} \right] \tag{5}$$ $N$ is then sought as solution of the homogeneous linear equation: $$ Q^{T}Q N = 0_3 \label{QtQNeq0} \tag{6}$$ with: $$ Q^{T}Q = \left[ {\begin{array}{ccc} \sum_{i=1}^n x_i^2 & \sum_{i=1}^n x_i y_i & \sum_{i=1}^n x_i z_i \\ \sum_{i=1}^n y_i x_i & \sum_{i=1}^n y_i^2 & \sum_{i=1}^n y_i z_i \\ \sum_{i=1}^n z_i x_i & \sum_{i=1}^n z_i y_i & \sum_{i=1}^n z_i^2 \end{array}} \right] \tag{7}$$ $Q^{T}Q$ is called the covariance matrix.
Let $\vec{V_x}, \vec{V_y}, \vec{V_z}$ be vectors defined as: $$ \vec{V_x} = \left( {\begin{array}{c} \sum_{i=1}^n x_i^2 \\ \sum_{i=1}^n x_i y_i \\ \sum_{i=1}^n x_i z_i \end{array}} \right), \vec{V_y} = \left( {\begin{array}{c} \sum_{i=1}^n y_i x_i \\ \sum_{i=1}^n y_i^2 \\ \sum_{i=1}^n y_i z_i \end{array}} \right), \vec{V_z} = \left( {\begin{array}{c} \sum_{i=1}^n z_i x_i \\ \sum_{i=1}^n z_i y_i \\ \sum_{i=1}^n z_i^2 \end{array}} \right) \tag{8}$$ Equation ($\ref{QtQNeq0}$) is then rewritten as: $$ \left[ {\begin{array}{c} \vec{V_x} \\ \vec{V_y} \\ \vec{V_z} \end{array}} \right] N = 0 \label{VNeq0} \tag{9}$$ In vectorial formalism, Equation ($\ref{VNeq0}$) reads: $$ \left\{ {\begin{array}{c} \vec{V_x} \cdot \vec{n} = 0 \\ \vec{V_y} \cdot \vec{n} = 0 \\ \vec{V_z} \cdot \vec{n} = 0 \end{array}} \right. \label{VxyzdotNeq0} \tag{10}$$ It states that $\vec{n}$ should be perpendicular to the space spanned by the three vectors $(\vec{V_x}, \vec{V_y}, \vec{V_z})$.
As all the points $M_i$ and $C$ are coplanar, it is known that this space is degenerated of dimension 2, even of dimension 1 if the points are aligned. It means that without loosing generality, one might assume that $\vec{V_z}$ is a linear combination of the two others: $$ \vec{V_z} = \alpha \vec{V_x} + \beta \vec{V_y}\tag{11}$$ So: $$ \vec{V_z} \cdot \vec{n} = \alpha \vec{V_x} \cdot \vec{n} + \beta \vec{V_y} \cdot \vec{n} \tag{12}$$ The last line of Equation ($\ref{VxyzdotNeq0}$) is thus just a consequence of the first two lines. Remains only: $$ \left\{ {\begin{array}{c} \vec{V_x} \cdot \vec{n} = 0 \\ \vec{V_y} \cdot \vec{n} = 0 \end{array}} \right. \label{VxydotNeq0} \tag{13}$$ It states that $\vec{n}$ should be perpendicular to $\vec{V_x}$ and $\vec{V_y}$. As a consequence $\vec{n}$ should be collinear to the cross product "$\times$" of both: $$ \vec{n} \ \| \ \vec{V_x} \times \vec{V_y} \label{VxCrossVy} \tag{14}$$ If the points are aligned, only one line of Equation ($\ref{VxyzdotNeq0}$) remains, and $\vec{n}$ is undetermined.
To not rely on an arbitrary choice of the two vectors in Equation ($\ref{VxCrossVy}$), Gram-Schmidt process is use to orthonormalize the set of vectors $(\vec{V_x}, \vec{V_y}, \vec{V_z})$:

• The first vector labelled $\vec{e_1}$ is chosen among $(\vec{V_x}, \vec{V_y}, \vec{V_z})$ as the one maximizing its norm ; it is normalized.
• Then all three vectors are orthogonalized to $\vec{e_1}$ by: $$ \vec{U} = \vec{V} - (\vec{V} \cdot \vec{e_1}) \ \vec{e_1} \tag{15}$$
• The second vector labelled $\vec{e_2}$ is chosen among $(\vec{U_x}, \vec{U_y}, \vec{U_z})$ as the one maximizing its norm ; it is normalized.

Finally, the normal $\vec{n}$ is calculated from: $$\vec{n} = \vec{e_1} \times \vec{e_2} \tag{16}$$

Implementation:

Step 1: Assembly of the covariance matrix
Objective: Compute the set of vectors $(\vec{V_x}, \vec{V_y}, \vec{V_z})$. 1.1. Computation of $$\vec{CM_i} = \left( {\begin{array}{c} x_i \\ y_i \\ z_i \end{array}} \right)$$ 1.2. Computation of the vector $$ x_i \left( {\begin{array}{c} x_i \\ y_i \\ z_i \end{array}} \right) = \left( {\begin{array}{c} x_i^2 \\ x_i y_i \\ x_i z_i \end{array}} \right) $$ 1.3. Summation of the output of step 1.2 using the Group ID to accumulate the values in different bins, one per Group ID. So each sum is restricted to the points sharing the same Group ID, yielding: $$\vec{V_x} = \left( {\begin{array}{c} \sum_{i=1}^n x_i^2 \\ \sum_{i=1}^n x_i y_i \\ \sum_{i=1}^n x_i z_i \end{array}} \right)$$ 1.4. Computation of the vector $$ y_i \left( {\begin{array}{c} x_i \\ y_i \\ z_i \end{array}} \right) = \left( {\begin{array}{c} y_i x_i \\ y_i^2 \\ y_i z_i \end{array}} \right) $$ 1.5. Summation of the output of step 1.4 yielding: $$\vec{V_y} = \left( {\begin{array}{c} \sum_{i=1}^n y_i x_i \\ \sum_{i=1}^n y_i^2 \\ \sum_{i=1}^n y_i z_i \end{array}} \right)$$ 1.6. Computation of the vector $$ z_i \left( {\begin{array}{c} x_i \\ y_i \\ z_i \end{array}} \right) = \left( {\begin{array}{c} z_i x_i \\ z_i y_i \\ z_i^2 \end{array}} \right) $$ 1.7. Summation of the output of step 1.6 yielding: $$\vec{V_z} = \left( {\begin{array}{c} \sum_{i=1}^n z_i x_i \\ \sum_{i=1}^n z_i y_i \\ \sum_{i=1}^n z_i^2 \end{array}} \right)$$

Step 2: Gram-Schmidt normalization
Objective: Find among $(\vec{V_x}, \vec{V_y}, \vec{V_z})$ the vector with maximal norm, and normalize it. 2.1. Computation of the three norms.
2.2. Sorting of cases following the pseudo-code:

if |Vx| > |Vy|:
  if |Vx| > |Vz|:
    return Vx
  else
    return Vz
else
  if |Vy| > |Vz|:
    return Vy
  else
    return Vz

2.3. Normalization: $$ \vec{e} = {\vec{V}} \ / \ {\|\vec{V}\|} $$

Step 3: Gram-Schmidt orthogonalization
Objective: Keep only the part of vector $\vec{V}$ perpendicular to vector $\vec{e}$. 3.1. Computation of the projection of $\vec{V}$ on $\vec{e}$ as $\vec{V} \cdot \vec{e}$.
3.2. Computation of the part of $\vec{V}$ aligned with $\vec{e}$ as $(\vec{V} \cdot \vec{e}) \ \vec{e}$.
3.3. Computation of the part of $\vec{V}$ perpendicular to $\vec{e}$ by removing its part aligned with $\vec{e}$.

Step 4: Computation of the normal
Objective: Find vectors $\vec{e_1}$ and $\vec{e_2}$ by an incomplete orthonormalization of $(\vec{V_x}, \vec{V_y}, \vec{V_z})$, then compute the normal. 4.1. Computation of $\vec{e_1}$.
4.2. Subtraction of direction $\vec{e_1}$ from $(\vec{V_x}, \vec{V_y}, \vec{V_z})$.
4.3. Computation of $\vec{e_2}$.
4.4. Computation of the normal as $\vec{n} = \vec{e_1} \times \vec{e_2}$.

Result:

Resources:

Answer 2

I don't believe this is a solvable problem in its fully general form, where we can feed an algorithm a set of any points, in any order, in any position, and get the same face normal Blender provided for the face from which those points were derived. There's missing information, which is how those points are connected by edges.

We can see how, connecting two identical sets of points with different edges, Blender calculates different face normals.

From comments:

I would also be interested in a solution where the points per group are on one level.

If all points are on the same plane, then we can choose any three points (that are not collinear) to determine the plane, and from that, determine the normal of that plane.

Note that the origin of the displayed vector is the object origin, rather than the center of the three points.   Every "edge" is orthogonal to the normal, so if we have two edges that do not lie in the same line, then we can simply find a vector that is orthogonal to both of those edges via the cross-product. This can't tell us whether the normal points up or down-- I think you realize that and are okay with it.

There's also the problem of, what if the points we pick are collinear? You might already know that none of your edges are degenerate and that your faces are convex, but if not, we may end up with parallel vectors and a meaningless cross product. In that case, we can abandon any vectors that are parallel with the first one we choose:

I shift the whole thing to the space of our first point (arbitrarily chosen as index = 0) then check the normalized position of each point; I discard any points for which this vector is parallel (or antiparallel) to the vector to index = 1 (our arbitrarily chosen baseline vector.) I then discard our first point (its vector to itself in this space is the zero vector, useless, and its position is baked into the coordinate space anyways.) With what's left, I can take a cross-product to find the normal. I'm showing lines representing the normal as well as the two vectors used to generate it on the right; on the left are our source points, plenty of which are collinear.

In the case where there all of our points are collinear, this isn't a solvable problem-- the normal is undefined.

Answer 3

There is a statistical technique called Principal Component Analysis (PCA) that fits this problem. It helps in data analysis by expressing the data in a new coordinate system that captures the variability of data better than the original coordinate system.

In this problem, the data points are the vertices that make up the face; the original coordinate system is the Global or Local coordinate system; the new coordinate system is a rotated coordinate system aligned to the plane of the vertices. The first two axes of it capture the biggest variability of the data which means the third axis is the best candidate to point along (or in the opposite direction of) the normal.

As Nathan said, even if a direction is found, we will always have two options. We'll need another method to decide where the normal should be pointing (inside or outside). Normally, it will be possible to decide based on the volume of the shape as a whole or the neighboring faces.

Also as Nathan said, there is no assurance that the data points (vertices in this case) are distributed to have principal components. However, in scientific data analysis, there is usually a sound assumption that needs data to prove, so there should be a pattern to be captured with this technique. Similarly in this 3D problem, we have to assume that the data comes from a reasonable source (broken model, 3D scanning, etc.) that there is indeed an axis that has the least variability, and that the vertices are not collinear.

Additionally, if we set the vertices to point along one vector, and the face is actually not flat, it will look like it is flat in shading as light calculation depends on the normals.

Assuming that the data is valid, PCA can be successfully performed to obtain the normal direction. The steps are quite involved, but the main steps are listed below.

1. Normalize the data which requires calculating the mean and standard deviation of each variable (x, y, and z components). Normalization is done by subtracting the mean and dividing by the standard deviation.
2. Calculate the covariance matrix. It captures the correlation between each pair of variables.
3. Calculate the eigenvectors and eigenvalues. The eigenvector of the variable that has the smallest eigenvalue is the best option to use for the normal. The other two capture most of the variability of the data which mean that they form the approximate plane of the face.

This method is assured to give deterministic results and will give the best candidate vector for the normal vector.

I think this YouTube video explains the concept of this technique very well, and this one explains how to perform the calculation steps well too.

Sadly, the complexity of the calculation is prohibiting considering that it has to be done to each face. The most difficult step is the last one which involves solving 3rd degree equations. It is not the solution that OP wants, but I share it here in case someone finds it interesting or useful.

Answer 4

(Using Blender 3.6.5)

This proposal is based on the Convex Hull node applied on a subset of the input geometry extracted with a Separate Geometry node. The sought normal per subset is the normal to the first (i.e. Index 0) face. As illustration, it is draw from the subset barycentre.
Because I do not know how to iterate over the Group Index, it is controlled through the Group Input. Perhaps Blender 4.0 Repeat Zone could save this approach... Resources:

Answer 5

We know that the general case can not be solved.

But, we can try to do something if:

• all points of a group are not colinear
• the center point of a group is not at the same position as one point of this group (ok for any convex shapes). But we can test that, eliminate this point and continue without this point (this aspect is not implemented below).
• original face is planar
• no dup vertices
• other?

From that we need two points per group id, that can be the center C of this group and the nearest point N from this center.

Then we can calculate for each point P the cross product of PN and CN which is :

• either invalid if PN and CN are aligned (in that case we set temporary the normal to 0,0,0)
• the normal if not. But this normal can be the original normal or its opposite

Accumulate the normal and get leading / trailing difference, and we can compare that to the previous found normal.

This comparison is partially implemented here. But doing completely is comparing X, Y and Z component in sequence and return true for the first that meet the condition (so not only testing X, as in the image above).

We now have homogenuous normals, but some normal set to 0, 0, 0 may remain.

To eliminate them, we can accumulate again so that we have the sum of the normals, the sum of the non zero and get the mean.