Test the numerical stability of different techniques for computing normals.

stable_normals.py

from __future__ import print_function, division

from math import *
from numpy import *

from pprint import pprint

def get_rotation_matrix( axis, radians ):
    '''
    Given the axis, angle rotation defined by 'axis' and 'radians', returns a
    3x3 numpy.matrix M suitable for multiplying a 3-vector v: M*v.
    
    NOTE: the angle is in radians, not degrees as glRotate*() takes.
    
    EXAMPLE:
        ## A CCW rotation of pi/2 radians about the z axis in the xy plane. Maps <1,0,0> to <0,1,0>.
        dot( asarray( get_rotation_matrix( (0,0,1), pi/2 ) ), (1,0,0) )
        
        ## The matrix transforming a plane with arbitrary normal to the xy plane:
        asarray( get_rotation_matrix( *get_axis_angle_rotation( plane_normal, (0,0,1) ) ) )
    
    tested
    
    This function comes from vecutil.py: https://github.com/yig/vecutil/blob/master/vecutil.py
    '''
    
    axis = asfarray( axis )
    
    ## This code comes from hmatrix.h / adt/mat3t.h
    
    u = normalize( axis )
    
    u1 = u[0]
    u2 = u[1]
    u3 = u[2]
    
    U = asfarray([
        [u1*u1, u1*u2, u1*u3],
        [u2*u1, u2*u2, u2*u3],
        [u3*u1, u3*u2, u3*u3]
        ])
    
    S = asfarray([
        [0, -u3, u2],
        [u3,  0, -u1],
        [-u2,  u1, 0]
        ])
    
    result = U + (identity(3) - U) * cos(radians)  + S * sin(radians)
    return result

def magnitude2( n ):
    return dot( n, n )

def magnitude( n ):
    return sqrt( magnitude2( n ) )

def normalize( n ):
    n = n.copy()
    
    mag = magnitude( n )
    
    if mag != 0.:
        n /= mag
    
    return n

def random_unit_vector():
    ## Normalize any vector.
    ## What about all zeros?
    # return normalize( random.random(3) )
    
    while True:
        v = random.random(3)
        ## Skip the zero vector.
        if ( v == 0. ).all(): continue
        ## Reject points outside the unit sphere.
        if magnitude2( v ) > 1.: continue
        ## If we're still here, we are good. Return the normalized vector.
        return normalize( v )

def normal_naive0( f ):
    n = cross( f[1] - f[0], f[2] - f[0] )
    return n

def normal_naive1( f ):
    n = cross( f[2] - f[1], f[0] - f[1] )
    return n

def normal_naive2( f ):
    n = cross( f[0] - f[2], f[1] - f[2] )
    return n

def normal_all_average_before_normalizing( f ):
    ## Average without normalizing
    n1 = cross( f[1] - f[0], f[2] - f[0] )
    n2 = cross( f[2] - f[1], f[0] - f[1] )
    n3 = cross( f[0] - f[2], f[1] - f[2] )
    n = n1 + n2 + n3
    return n

def normal_all_average_after_normalizing( f ):
    n1 = normalize( cross( f[1] - f[0], f[2] - f[0] ) )
    n2 = normalize( cross( f[2] - f[1], f[0] - f[1] ) )
    n3 = normalize( cross( f[0] - f[2], f[1] - f[2] ) )
    n = n1 + n2 + n3
    return n

def normal_svd( f ):
    ## Subtract the centroid
    centroid = average( f, axis = 0 )
    centered = f - centroid
    U, sigma, V = linalg.svd( centered )
    ## SVD will return the unit axes if all singular values are 0.
    ## It would be cheating if the target normal were (0,0,1).
    if ( sigma < 1e-8 ).all(): print( 'SVD would be cheating if the target normal were (0,0,1).' )
    #print( 'U:', U )
    #print( 'sigma:', sigma )
    #print( 'V:', V )
    n = V[2]
    ## ALERT! We don't know whether it should be + or - n.
    return n

def test_many_normals( fs, ns_should_be, perps_should_be = None ):
    names = 'naive0', 'naive1', 'naive2', 'all_average_before_normalizing', 'all_average_after_normalizing', 'svd'
    diff_errors = { name: 0. for name in names }
    perp_errors = { name: 0. for name in names }
    ground_truth_perp_error = 0.
    
    if perps_should_be is None: perps_should_be = zeros( ( len( fs ), 0, 3 ) )
    
    assert len( fs ) == len( ns_should_be )
    assert len( fs ) == len( perps_should_be )
    N = len( fs )
    
    for f, n_should_be, perp_should_be in zip( fs, ns_should_be, perps_should_be ):
        for name in names:
            func = globals()[ 'normal_' + name ]
            prenormalized_n = func( f )
            ## A patch for svd, which doesn't know the sign:
            if dot( prenormalized_n, n_should_be ) < 0.: prenormalized_n = -prenormalized_n
            
            n = normalize( prenormalized_n )
            
            diff_err = abs( n - n_should_be ).sum()
            diff_errors[ name ] += diff_err
            
            for p in perp_should_be:
                assert prenormalized_n.shape == p.shape
                perp_err = abs( dot( prenormalized_n, p ) )
                mag2 = magnitude2( prenormalized_n )
                ## We do have to appropriately normalize this error.
                if mag2 > 0.:
                    perp_err *= 1./sqrt( mag2 )
                perp_errors[ name ] += perp_err
        
        ## Let's also evaluate n_should_be with perp_should_be
        for p in perp_should_be:
            perp_err = abs( dot( n_should_be, p ) )
            ground_truth_perp_error += perp_err
    
    print( '# Total L1 error |n - what n should be|:' )
    # pprint( errors )
    pprint( sorted( diff_errors.iteritems(), key = lambda e: e[1] ) )
    
    ## Average
    average_diff_errors = {}
    for name, error in diff_errors.iteritems():
        ## Divide by 3, because there are three components for the normals.
        assert 3 == ns_should_be.shape[1]
        average_diff_errors[ name ] = error / (3*N)
    
    print( '# Average L1 error |n - what n should be|:' )
    # pprint( errors )
    pprint( sorted( average_diff_errors.iteritems(), key = lambda e: e[1] ) )
    
    ###########
    ## If we have no perpendicular vectors, return.
    if perps_should_be.shape[1] == 0: return
    
    print( '# Total prenormalized perpendicular dot error |dot( prenormalized n, perpendicular vector )|/|prenormalized n|:' )
    # pprint( errors )
    pprint( sorted( perp_errors.iteritems(), key = lambda e: e[1] ) )
    
    ## Average
    average_perp_errors = {}
    for name, error in perp_errors.iteritems():
        average_perp_errors[ name ] = error / ( N * perps_should_be.shape[1] )
    
    print( '# Average prenormalized perpendicular dot error |dot( prenormalized n, perpendicular vector )|/|prenormalized n|:' )
    # pprint( errors )
    pprint( sorted( average_perp_errors.iteritems(), key = lambda e: e[1] ) )
    
    print( '# Ground truth normal perpendicular dot error |dot( true n, perpendicular vector )|:' )
    print( 'total:', ground_truth_perp_error )
    print( 'average:', ground_truth_perp_error / ( N * perps_should_be.shape[1] ) )

def test_one_normal( f ):
    names = 'naive0', 'naive1', 'naive2', 'all_average_before_normalizing', 'all_average_after_normalizing', 'svd'
    for name in names:
        func = globals()[ 'normal_' + name ]
        prenormalized_n = func( f )
        n = normalize( prenormalized_n )
        print( name, n, '(prenormalized:', prenormalized_n, ')' )

def test_thin( s = 100 ):
    print( '### test_thin( s = %r )' % ( s, ) )
    
    ## The following line allows us to take string input.
    s = float( s )
    
    f = asfarray([
        ( 0,0,0 ),
        ( 1,0,0 ),
        ( 1./s,1./s,0 )
        ])
    
    test_one_normal( f )

def test_random_rotation( N = 10000, s = 1e20 ):
    print( '### test_random_rotation( N = %r, s = %r )' % ( N, s ) )
    
    ## The following line allows us to take string input.
    N = int( N )
    s = float( s )
    
    ## Pick three random points on the z = 0 plane.
    fs = random.random( (N,3,3) )
    fs[:,:,0] /= s
    fs[:,:,2] = 0.
    
    ## Normals.
    ns_should_be = zeros( ( N, 3 ) )
    ns_should_be[:,2] = 1.
    
    ## Perpendicular vectors.
    perps_should_be = zeros( ( N, 5, 3 ) )
    perps_should_be[:,0,:] = (1,0,0)
    perps_should_be[:,1,:] = (0,1,0)
    
    for i in range( N ):
        ## Generate a random rotation matrix.
        axis = random_unit_vector()
        angle = pi*random.random()
        R = get_rotation_matrix( axis, angle )
        
        ## Rotate the face vertices, the normal, and the first two perpendicular vectors.
        fs[i] = [ dot( R, v ) for v in fs[i] ]
        ns_should_be[i] = dot( R, ns_should_be[i] )
        for j in range(2):
            perps_should_be[i][j] = dot( R, perps_should_be[i][j] )
    
    ## Three more perpendicular vectors are the triangles' edges.
    perps_should_be[:,2,:] = fs[:,2,:] - fs[:,1,:]
    perps_should_be[:,3,:] = fs[:,1,:] - fs[:,0,:]
    perps_should_be[:,4,:] = fs[:,0,:] - fs[:,2,:]
    
    ## Use only the triangles' edge vectors.
    perps_should_be = perps_should_be[:,2:,:]
    
    test_many_normals( fs, ns_should_be, perps_should_be )

def test_random_plane( N = 10000, s = 1e20 ):
    print( '### test_random_plane( N = %r, s = %r )' % ( N, s ) )
    
    ## The following line allows us to take string input.
    N = int( N )
    s = float( s )
    
    ## Generate random points on a random plane.
    ## If we choose a plane normal with z > 0, we can lift a point in the z=0 plane
    ## appropriately.
    
    ## Pick three random points on the z = 0 plane.
    fs = random.random( (N,3,3) )
    fs[:,:,0] /= s
    fs[:,:,2] = 0.
    
    ## Space for normals.
    ns_should_be = zeros( ( N, 3 ) )
    
    ## Space for perpendicular vectors.
    perps_should_be = zeros( ( N, 5, 3 ) )
    perps_should_be[:,0,:] = (1,0,0)
    perps_should_be[:,1,:] = (0,1,0)
    
    ## Generate a random unit normal with z>0, generate an offset, and lift the xy points.
    for i in range( N ):
        ## Get a normal vector with positive z.
        while True:
            normal = random_unit_vector()
            if normal[2] > .1: break
        
        ## Generate a random offset in [0,1).
        offset = random.random()
        
        ## The function returning the z value for a given x,y.
        def plane_z( xy ):
            return ( dot( xy, normal[:2] ) + offset )/-normal[2]
        
        ## Generate z values for the xy points.
        fs[i,:,2] = [ plane_z( v[:2] ) for v in fs[i] ]
        ## The normal is the plane's normal.
        ns_should_be[i] = normal
        
        ## The first two perpendicular vectors are the unit +x, +y directions on the plane.
        perps_should_be[i,0,2] = plane_z( perps_should_be[i,0,:2] )
        perps_should_be[i,1,2] = plane_z( perps_should_be[i,1,:2] )
    
    ## Three more perpendicular vectors are the triangles' edges.
    perps_should_be[:,2,:] = fs[:,2,:] - fs[:,1,:]
    perps_should_be[:,3,:] = fs[:,1,:] - fs[:,0,:]
    perps_should_be[:,4,:] = fs[:,0,:] - fs[:,2,:]
    
    ## Use only the triangles' edge vectors.
    perps_should_be = perps_should_be[:,2:,:]
    
    test_many_normals( fs, ns_should_be, perps_should_be )


def usage():
    import sys
    print( 'Usage:', sys.argv[0], 'thin [stretch-factor (default = 100)]', file = sys.stderr )
    print( 'Usage:', sys.argv[0], 'random_rotation [number-of-triangles (default = 2000) [stretch-factor (default = 1e20)]]', file = sys.stderr )
    print( 'Usage:', sys.argv[0], 'random_plane [number-of-triangles (default = 2000) [stretch-factor (default = 1e20)]]', file = sys.stderr )
    sys.exit(-1)

if __name__ == '__main__':
    import sys
    
    if len( sys.argv ) == 1:
        usage()
    
    which = sys.argv[1]
    func = globals()[ 'test_' + which ]
    func( *sys.argv[2:] )