How to Estimate the Normal of a Point Cloud

14 Feb 2026

This is something I spent multiple days on recently trying to figure out, so I thought I'd write it down. Hopefully it helps someone out.

# 1. Covariance matrix

First calculate the covariance matrix of the given set of points. This is fairly easy with the following pseudo-code:

usize n = points.len;

f32 mean = getMean(points);

f32 xx, xy, xz;
f32 yx, yy, yz;
f32 zx, zy, zz;
for i in n {
	p = points[i] - mean;
	xx += p.x * p.x;
	xy += p.x * p.y;
	xz += p.x * p.z;
	// And so on ...
}

nm1 = n - 1;
xx /= nm1;
xy /= nm1;
xz /= nm1;
// And so on ...

The covariance matrix is the covariance of x with x, x with y, x with z, and so on, over all of the points. Covariance can be defined as:

sum((x - avg(x))(y - avg(y)))
-----------------------------
           n - 1

# 2. Minimum eigenvalue

To get the normal, we need the minimum eigenvector. To get the minimum eigenvector, we can use the minimum eigenvalue.

Since the covariance matrix is guaranteed to be a symmetric 3x3 matrix, we can use the algorithm described here (we only need eig3, which is the minimum). I've copied the algorithm here:

% Given a real symmetric 3x3 matrix A, compute the eigenvalues
% Note that acos and cos operate on angles in radians

p1 = A(1,2)^2 + A(1,3)^2 + A(2,3)^2
if (p1 == 0)
   % A is diagonal.
   eig1 = A(1,1)
   eig2 = A(2,2)
   eig3 = A(3,3)
else
   q = trace(A)/3               % trace(A) is the sum of all diagonal values
   p2 = (A(1,1) - q)^2 + (A(2,2) - q)^2 + (A(3,3) - q)^2 + 2 * p1
   p = sqrt(p2 / 6)
   B = (1 / p) * (A - q * I)    % I is the identity matrix
   r = det(B) / 2

   % In exact arithmetic for a symmetric matrix  -1 <= r <= 1
   % but computation error can leave it slightly outside this range.
   if (r <= -1)
      phi = pi / 3
   elseif (r >= 1)
      phi = 0
   else
      phi = acos(r) / 3
   end

   % the eigenvalues satisfy eig3 <= eig2 <= eig1
   eig1 = q + 2 * p * cos(phi)
   eig3 = q + 2 * p * cos(phi + (2*pi/3))
   eig2 = 3 * q - eig1 - eig3     % since trace(A) = eig1 + eig2 + eig3
end

# 3. Minimum eigenvector

To get the eigenvector from the eigenvalue, we can use the cross product trick. The code would look something like this:

Mat3x3 v:  = ...; // Covariance matrix
f32 eig3 = ...; // Minimum eigen value

Mat3x3 m = v - eig3 * Matri3x3.identity;

// Return the cross product of two linearly independant rows,
// which will have the largest cross product
Vec3 c1 = cross(m.row(0), m.row(1));
Vec3 c2 = cross(m.row(1), m.row(2));
Vec3 c3 = cross(m.row(2), m.row(0));
f32 lc1 = c1.lengthSquared();
f32 lc2 = c2.lengthSquared();
f32 lc3 = c3.lengthSquared();

// Normalize the vector because it's a normal
if (lc1 > lc2 && lc1 > lc3) return lc1.normalized();
if (lc2 > lc1 && lc2 > lc3) return lc2.normalized();
else return lc3.normalized();

# Conclusion

In summary, all we're doing is getting the normalized minimum eigenvalue of the covariance matrix of the point cloud. There are useful NumPy functions for this, which most people use in their posts on this topic, but it took me ages to find the approriate algorithms to implement it myself. The Python code would look like this, where A is the list of points:

V = np.cov(A)
eigvals, eigvecs = np.linalg.eigh(V)
normal = eigvals[:, np.argmin(eigvecs)]
normal /= np.norm(normal)

← Back