# Rajarshi Guha <rajarshi@presidency.com>
# 6/5/2005
#
# Tessellate the unit sphere using recursive
# subdivision. Based on C code by Steve Hollasch 
# at http://stevehollasch.com/cgindex/geometry/sphere.c
# and the vertex/face specifications taken from
# http://eeg.sourceforge.net/eegdoc/eeg_toolbox/sphere_tri.html


import math, getopt, sys

class point:
    def __init__(self,tupl):
        self.x = tupl[0]
        self.y = tupl[1]
        self.z = tupl[2]
    def __str__(self):
        return ( "(%2.4f, %2.4f, %2.4f)" % (self.x, self.y, self.z) )

class triangle:
    def __init__(self,p1,p2,p3):
        self.p1 = p1
        self.p2 = p2
        self.p3 = p3


def midpoint(p1,p2):
    x = 0.5 * (p1.x + p2.x)
    y = 0.5 * (p1.y + p2.y)
    z = 0.5 * (p1.z + p2.z)
    return( point( (x,y,z) ) )

def rep_oct():
    v = [ 
        (1,0,0),
        (-1,0,0),
        (0,1,0),
        (0,-1,0),
        (0,0,1),
        (0,0,-1) ]

    oct_tri = [ 
        triangle( point(v[0]), point(v[4]), point(v[2]) ),
        triangle( point(v[2]), point(v[4]), point(v[1]) ),
        triangle( point(v[1]), point(v[4]), point(v[3]) ),
        triangle( point(v[3]), point(v[4]), point(v[0]) ),
        triangle( point(v[0]), point(v[2]), point(v[5]) ),
        triangle( point(v[2]), point(v[1]), point(v[5]) ),
        triangle( point(v[1]), point(v[3]), point(v[5]) ),
        triangle( point(v[3]), point(v[0]), point(v[5]) )]
    return oct_tri        

def rep_tet():
    sqrt3 = 0.5773502692
    v = [ 
          (sqrt3, sqrt3, sqrt3),
          (-sqrt3, -sqrt3, sqrt3),
          (-sqrt3,sqrt3,-sqrt3),
          (sqrt3,-sqrt3,-sqrt3)
        ]
    tet_tri = [
        triangle( point(v[0]),point(v[1]),point(v[2]) ),
        triangle( point(v[0]),point(v[3]),point(v[1]) ),
        triangle( point(v[2]),point(v[1]),point(v[3]) ),
        triangle( point(v[3]),point(v[0]),point(v[2]) )]
    return tet_tri

def rep_ico():
    tau = 0.8506508084
    one = 0.5257311121
    v = [
      (  tau,  one,    0 ),  
      ( -tau,  one,    0 ),  
      ( -tau, -one,    0 ),  
      (  tau, -one,    0 ),  
      (  one,   0 ,  tau ),  
      (  one,   0 , -tau ),  
      ( -one,   0 , -tau ),  
      ( -one,   0 ,  tau ),  
      (   0 ,  tau,  one ),  
      (   0 , -tau,  one ),  
      (   0 , -tau, -one ),  
      (   0 ,  tau, -one ) ]
    ico_tri = [
	triangle( point(v[4]), point(v[8]), point(v[7]) ),
	triangle( point(v[4]), point(v[7]), point(v[9]) ),
	triangle( point(v[5]), point(v[6]), point(v[11]) ),
	triangle( point(v[5]), point(v[10]), point(v[6]) ),
	triangle( point(v[0]), point(v[4]), point(v[3]) ),
	triangle( point(v[0]), point(v[3]), point(v[5]) ),
	triangle( point(v[2]), point(v[7]), point(v[1]) ),
	triangle( point(v[2]), point(v[1]), point(v[6]) ),
	triangle( point(v[8]), point(v[0]), point(v[11]) ),
	triangle( point(v[8]), point(v[11]), point(v[1]) ),
	triangle( point(v[9]), point(v[10]), point(v[3]) ),
	triangle( point(v[9]), point(v[2]), point(v[10]) ),
	triangle( point(v[8]), point(v[4]), point(v[0]) ),
	triangle( point(v[11]), point(v[0]), point(v[5]) ),
	triangle( point(v[4]), point(v[9]), point(v[3]) ),
	triangle( point(v[5]), point(v[3]), point(v[10]) ),
	triangle( point(v[7]), point(v[8]), point(v[1]) ),
	triangle( point(v[6]), point(v[1]), point(v[11]) ),
	triangle( point(v[7]), point(v[2]), point(v[9]) ),
	triangle( point(v[6]), point(v[10]), point(v[2]) )]
    return ico_tri

def normalize(p):
    mag = p.x*p.x + p.y*p.y + p.z*p.z
    if (mag != 0):
        mag = 1.0 / math.sqrt(mag)
        x = p.x * mag
        y = p.y * mag
        z = p.z * mag
        return( point( (x,y,z) ) )
    else: 
        return(p)
    
def usage():
    print """
    Usage: python tess.py [OPTIONS]

    Generates the points on a unit sphere obtained by recursive
    subdivision of the triangles of an initial representation which
    may be tetrahedral, octagonal or icosahedral.

    By default the levels of subdivision is set to 4 and the points
    are dumped to stdout. The resultant vertices are clockwise ordered

    -h,--help       This message
    -l,--level      Number of levels of subdivision
    -o,--output     A file to dump the coordinates to
    -s,--scale      A factor by which the final coordinates are scaled (default is 1)
    -t,--type       Type of initial representation: tet, oct or ico
    """

if __name__ == '__main__':

    maxlevel = 4
    outfile = ''
    reptype = 'oct'
    scalefactor = 1
    
    try:
        opt,args = getopt.getopt(sys.argv[1:], 't:o:l:hs:',\
        ['level=','help','output=', 'type=', 'scale='])
    except getopt.GetoptError:
        usage()
        sys.exit(0)
        
    for o,a in opt:
        if o in ('-h','--help'):
            usage()
            sys.exit(0)
        if o in ('-l','--level'):
            maxlevel = int(a)
        if o in ('-o','--output'):
            outfile = a
        if o in ('-s','--scale'):
            scalefactor = float(a)
        if o in ('-t','--type'):
            if a not in ['tet','oct','ico']:
                usage()
                sys.exit()
            else:
                reptype = a

    oldtess = None
    if reptype == 'tet':
        oldtess = rep_tet()
    elif reptype == 'oct':
        oldtess = rep_oct()
    elif reptype == 'ico':
        oldtess = rep_ico()

    # Start the subdivision
    for j in range(1, maxlevel):

        oldN = len(oldtess)
        newN = len(oldtess)*4
        newtess = [ triangle( point((0,0,0)), point((0,0,0)), point((0,0,0)) ) ] * newN

        # Subdivide each triangle 
        for i  in range(0, oldN):
            old = oldtess[i]

            p1 = normalize(midpoint(old.p1, old.p3))
            p2 = normalize(midpoint(old.p1, old.p2))
            p3 = normalize(midpoint(old.p2, old.p3))

            newtess[i*4] = triangle( old.p1, p2, p1 )
            newtess[i*4 + 1] = triangle( p2, old.p2, p3 )
            newtess[i*4 + 2] = triangle( p1,p2,p3 )
            newtess[i*4 + 3] = triangle( p1,p3,old.p3 )

        oldtess = newtess            
    
    print 'Number of triangles = %d' % (len(oldtess))


    if scalefactor != 1:
        for i in range(0,len(oldtess)):
            oldtess[i].p1.x = oldtess[i].p1.x * scalefactor
            oldtess[i].p1.y = oldtess[i].p1.y * scalefactor
            oldtess[i].p1.z = oldtess[i].p1.z * scalefactor

            oldtess[i].p2.x = oldtess[i].p2.x * scalefactor
            oldtess[i].p2.y = oldtess[i].p2.y * scalefactor
            oldtess[i].p2.z = oldtess[i].p2.z * scalefactor

            oldtess[i].p3.x = oldtess[i].p3.x * scalefactor
            oldtess[i].p3.y = oldtess[i].p3.y * scalefactor
            oldtess[i].p3.z = oldtess[i].p3.z * scalefactor

    if (outfile != ''):
        f = open(outfile,'w')
        for tri in oldtess:
            f.write("%f %f %f\n" % (tri.p1.x, tri.p1.y, tri.p1.z))
            f.write("%f %f %f\n" % (tri.p2.x, tri.p2.y, tri.p2.z))
            f.write("%f %f %f\n" % (tri.p3.x, tri.p3.y, tri.p3.z))
        f.close()            
    else:            
        for tri in oldtess:
            print tri.p1
            print tri.p2
            print tri.p3
            print