Compare the accuracy of various techniques for computing a triangle's area

triangle_area.py

from __future__ import print_function, division

from math import sqrt

'''
The Orient2D, Area2D, and ShewchukArea3D functions are
from "Lecture Notes on Geometric Robustness" by Jonathan Richard Shewchuk (2013).
URL: http://www.cs.berkeley.edu/~jrs/meshpapers/robnotes.pdf
'''

def Orient2D( a, b, c ):
    '''
    The determinant | ax ay 1; bx by 1; cx cy 1 | computed stably as per
    page 10 of Shewchuk.
    '''
    return ( a[0] - c[0] )*( b[1] - c[1] ) - ( b[0] - c[0] )*( a[1] - c[1] )

def ShewchukArea2D( a, b, c ):
    '''
    Given the 2D vertices of a triangle, return its area.
    '''
    
    assert len( a ) == 2
    assert len( a ) == len( b )
    assert len( a ) == len( c )
    
    return 0.5*Orient2D( a,b,c )

def ShewchukArea3D( a, b, c ):
    '''
    Given the 3D vertices of a triangle, return its area.
    '''
    
    assert len( a ) == 3
    assert len( a ) == len( b )
    assert len( a ) == len( c )
    
    return 0.5*sqrt(
        Orient2D( a[1:], b[1:], c[1:] )**2
        +
        Orient2D( ( a[2],a[0] ), ( b[2],b[0] ), ( c[2],c[0] ) )**2
        +
        Orient2D( a[:2], b[:2], c[:2] )**2
        )

from numpy import dot, asfarray, cross

def mag( v ): return sqrt( dot( v,v ) )
def mag2( v ): return dot( v,v )

def sqrt( x ):
    import math
    return float('nan') if x < 0 else math.sqrt(x)

def AreaHeronStableGivenVertices( v1, v2, v3 ):
    v1 = asfarray( v1 )
    v2 = asfarray( v2 )
    v3 = asfarray( v3 )
    
    a = mag( v2 - v1 )
    b = mag( v3 - v2 )
    c = mag( v1 - v3 )
    
    ## For some reason this makes things worse.
    # c,b,a = list(sorted([ a,b,c ]))
    
    return 0.25 * sqrt( (a+(b+c))*(c-(a-b))*(c+(a-b))*(a+(b-c)) )

def AreaHeronBasicGivenVertices( v1, v2, v3 ):
    v1 = asfarray( v1 )
    v2 = asfarray( v2 )
    v3 = asfarray( v3 )
    
    a = mag( v2 - v1 )
    b = mag( v3 - v2 )
    c = mag( v1 - v3 )
    
    s = 0.5*( a + b + c )
    
    return sqrt( s*( s - a )*( s - b )*( s - c ) )

def AreaHeronAlt2GivenVertices( v1, v2, v3 ):
    v1 = asfarray( v1 )
    v2 = asfarray( v2 )
    v3 = asfarray( v3 )
    
    a2 = mag2( v2 - v1 )
    b2 = mag2( v3 - v2 )
    c2 = mag2( v1 - v3 )
    
    return 0.25*sqrt( 2*(a2*b2 + a2*c2 + b2*c2 ) - ( a2**2 + b2**2 + c2**2 ) )

def AreaHeronAlt3GivenVertices( v1, v2, v3 ):
    v1 = asfarray( v1 )
    v2 = asfarray( v2 )
    v3 = asfarray( v3 )
    
    a2 = mag2( v2 - v1 )
    b2 = mag2( v3 - v2 )
    c2 = mag2( v1 - v3 )
    
    return 0.25*sqrt( (a2 + b2 + c2)**2 - 2*(a2**2 + b2**2 + c2**2 ) )

def AreaHeronAlt4GivenVertices( v1, v2, v3 ):
    v1 = asfarray( v1 )
    v2 = asfarray( v2 )
    v3 = asfarray( v3 )
    
    a2 = mag2( v2 - v1 )
    b2 = mag2( v3 - v2 )
    c2 = mag2( v1 - v3 )
    
    return 0.25*sqrt( 4*a2 * b2 - (a2 + b2 - c2)**2 )

def AreaCrossProduct0( v0, v1, v2 ):
    v0 = asfarray( v0 )
    v1 = asfarray( v1 )
    v2 = asfarray( v2 )
    
    return abs( .5*mag( cross( v1 - v0, v2 - v0 ) ) )

def AreaCrossProduct1( v0, v1, v2 ):
    v0 = asfarray( v0 )
    v1 = asfarray( v1 )
    v2 = asfarray( v2 )
    
    return abs( .5*mag( cross( v2 - v1, v0 - v1 ) ) )

def AreaCrossProduct2( v0, v1, v2 ):
    v0 = asfarray( v0 )
    v1 = asfarray( v1 )
    v2 = asfarray( v2 )
    
    return abs( .5*mag( cross( v0 - v2, v1 - v2 ) ) )

def test():
    ## Data from: Miscalculating Area and Angles of a Needle-like Triangle (from Lecture Notes for Introductory Numerical Analysis Classes) (W. Kahan 2014)
    ## URL: https://www.cs.berkeley.edu/~wkahan/Triangle.pdf
    
    lines = '''10 10 10 43.30127019 43.30127020 60 60 60
-3 4 2 2.905 -31337 151.045 151.045 -31337
100000 99999.99979 0.00029 17.6 9.999999990 0 2.02E-7 1.145915591E-7
100000 100000 1.00005 50010.0 50002.50003 0 5.73072E-4 5.72986443E-4
99999.99996 99999.99994 0.00003 -31337 1.118033988 0 -31337 1.281172578E-8
99999.99996 0.00003 99999.99994 -31337 1.118033988 48.18968509 -31337 48.18968510
10000 5000.000001 15000 0 612.3724358 180.000 180.000 179.9985965
99999.99999 99999.99999 200000 0 -31337 180 180 -31337
5278.64055 94721.35941 99999.99996 -31337 0 -31337 -31337 180
100002 100002 200004 0 0 -31337 180 180
31622.77662 0.000023 31622.77661 0.447 0.327490458 90 70.5 64.22853824
31622.77662 0.0155555 31622.77661 246.18 245.9540000 90.00 89.963187 89.96315276'''
    
    triangles = []
    for line in lines.split('\n'):
        a,b,c,HeronBasicArea,HeronStableArea,Cpp,Cp,Cdegrees = list(map( float, line.split() ))
        triangles.append({ 'a': a, 'b': b, 'c': c, 'HeronBasicArea': HeronBasicArea, 'HeronStableArea': HeronStableArea, 'Cdegrees': Cdegrees })
    
    from math import pi, cos, sin
    
    def triangle_vertices( a,b,c, Cdegrees ):
        Crad = Cdegrees*pi/180.
        v0 = asfarray( (0,0,0) )
        v1 = asfarray( (a,0,0) )
        v2 = asfarray( (b*cos(Crad),b*sin(Crad),0) )
        return v0, v1, v2
    
    names = 'AreaHeronBasicGivenVertices', 'AreaHeronStableGivenVertices', 'AreaHeronAlt2GivenVertices', 'AreaHeronAlt3GivenVertices', 'AreaHeronAlt4GivenVertices', 'AreaCrossProduct0', 'AreaCrossProduct1', 'AreaCrossProduct2', 'ShewchukArea3D'
    errors = { name: 0. for name in names }
    
    for i, tri in enumerate( triangles ):
        print( '--- Triangle', i, '---' )
        
        vs = triangle_vertices( tri['a'], tri['b'], tri['c'], tri['Cdegrees'] )
        sides_original = tri['a'], tri['b'], tri['c']
        sides_recomputed = mag( vs[1] - vs[0] ), mag( vs[0] - vs[2] ), mag( vs[2] - vs[1] )
        count_error = abs( asfarray( sides_original ) - sides_recomputed ).sum() < 1e-1
        
        print( 'Sides (original  ):', *sides_original )
        print( "Sides (recomputed):", *sides_recomputed )
        print( "[Good match]" if count_error else "[Bad match]" )
        print( 'Vertices:\n\t%r\n\t%r\n\t%r' % ( tuple( vs[0] ), tuple( vs[1] ), tuple( vs[2] ) ) )
        
        Aaccurate = tri['HeronStableArea']
        print( '%- 30s:' % 'Accurate area', tri['HeronStableArea'] )
        print( '%- 30s:' % 'Predicted Unstable Heron area', tri['HeronBasicArea'] )
        
        print()
        
        for name in names:
            A = globals()[ name ]( *vs )
            print( '%- 30s:' % name, '%r' % A )
            
            if count_error: errors[ name ] += abs( A - Aaccurate )
        
        print()
    
    print( '--- Total error ---' )
    from pprint import pprint
    pprint( errors )

if __name__ == '__main__':
    test()