Calculating surface normal in Python using Newell's Method

Question

I'm trying to implement Newell's Method to calculate the surface normal vector in Python, based on the following pseudocode from here.

Begin Function CalculateSurfaceNormal (Input Polygon) Returns Vector

   Set Vertex Normal to (0, 0, 0)

   Begin Cycle for Index in [0, Polygon.vertexNumber)

      Set Vertex Current to Polygon.verts[Index]
      Set Vertex Next    to Polygon.verts[(Index plus 1) mod Polygon.vertexNumber]

      Set Normal.x to Sum of Normal.x and (multiply (Current.y minus Next.y) by (Current.z plus Next.z))
      Set Normal.y to Sum of Normal.y and (multiply (Current.z minus Next.z) by (Current.x plus Next.x))
      Set Normal.z to Sum of Normal.z and (multiply (Current.x minus Next.x) by (Current.y plus Next.y))

   End Cycle

   Returning Normalize(Normal)

End Function

Here's my code:

Point3D = collections.namedtuple('Point3D', 'x y z')

def surface_normal(poly):
    n = [0.0, 0.0, 0.0]

    for i, v_curr in enumerate(poly):
        v_next = poly[(i+1) % len(poly)]
        n[0] += (v_curr.y - v_next.y) * (v_curr.z - v_next.z)
        n[1] += (v_curr.z - v_next.z) * (v_curr.x - v_next.x)
        n[2] += (v_curr.x - v_next.x) * (v_curr.y - v_next.y)

    normalised = [i/sum(n) for i in n]

    return normalised

def test_surface_normal():
    poly = [Point3D(0.0, 0.0, 0.0),
            Point3D(0.0, 1.0, 0.0),
            Point3D(1.0, 1.0, 0.0),
            Point3D(1.0, 0.0, 0.0)]

    assert surface_normal(poly) == [0.0, 0.0, 1.0]

This fails at the normalisation step since the n at that point is [0.0, 0.0, 0.0]. If I'm understanding correctly, it should be [0.0, 0.0, 1.0] (confirmed by Wolfram Alpha).

What am I doing wrong here? And is there a better way of calculating surface normals in python? My polygons will always be planar so Newell's Method isn't absolutely necessary if there's another way.

Answer 1

I thought it worth clarifying a few points about the informative implementation of Newell's method above. On first reading https://stackoverflow.com/users/1706564/jamie-bull update I thought the change of sign referred to the first equation only. But in fact it is all of them. Second, to provide a comparison to cross-product method. Finally, deal with divide by zero issue - but could just divide by epsilon. I also used numpy instead.

import numpy as np

def surface_normal_newell(poly):

    n = np.array([0.0, 0.0, 0.0])

    for i, v_curr in enumerate(poly):
        v_next = poly[(i+1) % len(poly),:]
        n[0] += (v_curr[1] - v_next[1]) * (v_curr[2] + v_next[2]) 
        n[1] += (v_curr[2] - v_next[2]) * (v_curr[0] + v_next[0])
        n[2] += (v_curr[0] - v_next[0]) * (v_curr[1] + v_next[1])

    norm = np.linalg.norm(n)
    if norm==0:
        raise ValueError('zero norm')
    else:
        normalised = n/norm

    return normalised


def surface_normal_cross(poly):

    n = np.cross(poly[1,:]-poly[0,:],poly[2,:]-poly[0,:])

    norm = np.linalg.norm(n)
    if norm==0:
        raise ValueError('zero norm')
    else:
        normalised = n/norm

    return normalised


def test_surface_normal3():
    """ should return:

        Newell:
        Traceback (most recent call last):
          File "demo_newells_surface_normals.py", line 96, in <module>
            test_surface_normal3()
          File "demo_newells_surface_normals.py", line 58, in test_surface_normal3
            print "Newell:", surface_normal_newell(poly) 
          File "demo_newells_surface_normals.py", line 24, in surface_normal_newell
            raise ValueError('zero norm')
        ValueError: zero norm
    """
    poly = np.array([[1.0,0.0,0.0],
                     [1.0,0.0,0.0],
                     [1.0,0.0,0.0]])
    print "Newell:", surface_normal_newell(poly) 


def test_surface_normal2():
    """ should return:

        Newell: [ 0.08466675 -0.97366764 -0.21166688]
        Cross : [ 0.08466675 -0.97366764 -0.21166688]
    """
    poly = np.array([[6.0,1.0,4.0],
                     [7.0,0.0,9.0],
                     [1.0,1.0,2.0]])
    print "Newell:", surface_normal_newell(poly)
    print "Cross :", surface_normal_cross(poly)


def test_surface_normal1():
    """ should return:

        Newell: [ 0.  0. -1.]
        Cross : [ 0.  0. -1.]
    """
    poly = np.array([[0.0,1.0,0.0],
                     [1.0,1.0,0.0],
                     [1.0,0.0,0.0]])
    print "Newell:", surface_normal_newell(poly) 
    print "Cross :", surface_normal_cross(poly)


print "Test 1:"
test_surface_normal1()
print "\n"

print "Test 2:"
test_surface_normal2()
print "\n"

print "Test 3:"
test_surface_normal3()

Answer 2

If you want an alternate approach to Newell's Method you can use the Cross-Product of 2 non-parallel vectors. That should work for just about any planar shape you'd provide it. I know the theory says it is for convex polygons, but examples we've looked at on Wolfram Alpha return an appropriate surface normal for even concave polygons (e.g. a bowtie polygon).

Answer 3

Ok, the problem was actually a stupid one.

The lines like:

n[0] += (v_curr.y - v_next.y) * (v_curr.z - v_next.z)

should be:

n[0] += (v_curr.y - v_next.y) * (v_curr.z + v_next.z) 

The values in the second set of brackets should be added, not subtracted.