Re: covering a sphere with patches

from visual import *

import math
from math import *

#
# general vector and quaternion functions
#
def add(v0, v1):
	return [ v1[0] + v0[0], v1[1] + v0[1], v1[2] + v0[2] ]

def sub(v0, v1):
	return [ v1[0] - v0[0], v1[1] - v0[1], v1[2] - v0[2] ]

def mag(v):
	return math.sqrt(v[0] * v[0] + v[1] * v[1] + v[2] * v[2] )

def mul(scalar, v):
	return [ scalar * v[0], scalar * v[1], scalar * v[2] ]

def normalize(v):
	return mul(1/mag(v), v)

def cross(v0, v1):
	return [ v0[1]*v1[2] - v0[2]*v1[1], v0[2]*v1[0] - v0[0]*v1[2], v0[0]*v1[1] - v0[1]*v1[0] ]

def dot(v0, v1):
	return v0[0] * v1[0] + v0[1] * v1[1] + v0[2] * v1[2]

def createQuaternionFromAngleAndAxis(angle, axis):
	angle = angle / 2
	axis = mul(math.sin(angle), normalize(axis))
	return axis + [ math.cos(angle) ]
	
# return a quarternion, which represents the rotation which is needed
# to get the target vector, if the source vector is rotated
def angleBetween(source, target):
	# the following code is based on an article by Stan Melax in Game Programming Gems,
	# see this book: http://books.google.de/books?id=hiBFUv_FT0wC&pg=PA214
	# first calculate a quarternion, which represents the rotation between v0 and v1
	
	v0 = normalize(source)
	v1 = normalize(target)
	d = dot(v0, v1)
	
	# special case: if the dot product is near -1, v0 and v1 are pointing
	# in opposite directions and we can use any perpendicular axis and 180 degree
	if d < 1e-5 - 1:
		axis = cross([0, 1, 0], v0)
		if mag(axis) < 1e-5:
			axis = cross([0, 0, 1], v0)
		return createQuaternionFromAngleAndAxis(pi, axis)
	else:
		# normal case
		s = sqrt((1+d) * 2)
		c = cross(v0, v1)
		qx = c[0] / s
		qy = c[1] / s
		qz = c[2] / s
		qw = s / 2
		return [ qx, qy, qz, qw ]

# return a new vector, which is the vector v rotated by the quaternion q
def rotate(q, v):
	rw = -q[0]*v[0]-q[1]*v[1]-q[2]*v[2]
	rx = q[3]*v[0]+q[1]*v[2]-q[2]*v[1]
	ry = q[3]*v[1]+q[2]*v[0]-q[0]*v[2]
	rz = q[3]*v[2]+q[0]*v[1]-q[1]*v[0]
	return [
		-rw*q[0]+rx*q[3]-ry*q[2]+rz*q[1],
		-rw*q[1]+ry*q[3]-rz*q[0]+rx*q[2],
		-rw*q[2]+rz*q[3]-rx*q[1]+ry*q[0] ]

# scale a vector
# v: vector
# length: destination length
def scaleToLength(v, length):
	return mul(length / mag(v), v)


#
# general mesh functions
#

# translate all points of the specified array
def translateArray(coords, v):
	coords2 = []
	for i in range(len(coords)):
		coords2 = coords2 + [ add(v, coords[i]) ]
	return coords2

# scale all points of the specified array to the specified length
def scaleArrayToLength(coords, length):
	coords2 = []
	for i in range(len(coords)):
		coords2 = coords2 + [ scaleToLength(coords[i], length) ]
	return coords2
	
# triangulate a mesh, if a polygon has more than 4 vertices (convex polygons, only)
def triangulate(coords, faces):
	faces2 = []
	for i in range(len(faces)):
		face = faces[i]
		if len(face) > 4:
			center = [0, 0, 0]
			for j in range(len(face)):
				coord = coords[face[j]]
				center[0] = center[0] + coord[0]
				center[1] = center[1] + coord[1]
				center[2] = center[2] + coord[2]
			center[0] = center[0] / len(face)
			center[1] = center[1] / len(face)
			center[2] = center[2] / len(face)
			centerIndex = len(coords)
			coords = coords + [center]
			for j in range(len(face)):
				oldIndex1 = face[j]
				if j == len(face) - 1:
					oldIndex2 = face[0]
				else:
					oldIndex2 = face[j + 1]
				faces2 = faces2 + [[oldIndex1, oldIndex2, centerIndex]]
		else:
			faces2 = faces2 + [face]
	return coords, faces2

# rotate all points of the specified array
def rotateArray(coords, q):
	coords2 = []
	for i in range(len(coords)):
		coords2 = coords2 + [ rotate(q, coords[i]) ]
	return coords2


# subdivide a polyhedron
# coords: coordinates
# faces: faces of the polyhedron (must be triangles)
# radius: radius of the sphere in which the polyhedron is inscribed
def subdivide(coords, faces, radius):
	lineHash = {}
	coords2 = coords
	faces2 = []
	for i in range(len(faces)):
		# get face coordinate indices
		face = faces[i]
		iv0 = face[0]
		iv1 = face[1]
		iv2 = face[2]
		
		# get face coordinates
		v0 = coords[iv0]
		v1 = coords[iv1]
		v2 = coords[iv2]
		
		# build line keys
		if iv0 < iv1:
			l1 = iv0 * len(coords) + iv1
		else:
			l1 = iv1 * len(coords) + iv0
		if iv1 < iv2:
			l2 = iv1 * len(coords) + iv2
		else:
			l2 = iv2 * len(coords) + iv1
		if iv0 < iv2:
			l3 = iv0 * len(coords) + iv2
		else:
			l3 = iv2 * len(coords) + iv0
			
		# try to find mid point index of a line in line hashtable
		if l1 in lineHash:
			iv0v1 = lineHash[l1]
		else:
			v0v1 = scaleToLength(mul(0.5, add(v0, v1)), radius)
			iv0v1 = len(coords2)
			lineHash[l1] = iv0v1
			coords2 = coords2 + [ v0v1 ]
		if l2 in lineHash:
			iv1v2 = lineHash[l2]
		else:
			v1v2 = scaleToLength(mul(0.5, add(v1, v2)), radius)
			iv1v2 = len(coords2)
			lineHash[l2] = iv1v2
			coords2 = coords2 + [ v1v2 ]
		if l3 in lineHash:
			iv2v0 = lineHash[l3]
		else:
			v2v0 = scaleToLength(mul(0.5, add(v2, v0)), radius)
			iv2v0 = len(coords2)
			lineHash[l3] = iv2v0
			coords2 = coords2 + [ v2v0 ]
			
		# create new faces
		faces2 = faces2 + [ [ iv0, iv0v1, iv2v0 ] ]
		faces2 = faces2 + [ [ iv1, iv1v2, iv0v1 ] ]
		faces2 = faces2 + [ [ iv2, iv2v0, iv1v2] ]
		faces2 = faces2 + [ [ iv0v1, iv1v2, iv2v0] ]
	return coords2, faces2

# keep all polygons with z coordinate less than the limit parameter
def half(coords, faces, limit):
	coords2 = []
	faces2 = []
	coordsHash = {}
	for i in range(len(faces)):
		face = faces[i]
		iv0 = face[0]
		iv1 = face[1]
		iv2 = face[2]
		v0 = coords[iv0]
		v1 = coords[iv1]
		v2 = coords[iv2]
		if v0[2] < limit or v1[2] < limit or v2[2] < limit:
			# map old coordinates to new coordinates to avoid duplicates
			if iv0 in coordsHash:
				iv0 = coordsHash[iv0]
			else:
				coordsHash[iv0] = len(coords2)
				iv0 = len(coords2)
				coords2 = coords2 + [ v0 ]
			if iv1 in coordsHash:
				iv1 = coordsHash[iv1]
			else:
				coordsHash[iv1] = len(coords2)
				iv1 = len(coords2)
				coords2 = coords2 + [ v1 ]
			if iv2 in coordsHash:
				iv2 = coordsHash[iv2]
			else:
				coordsHash[iv2] = len(coords2)
				iv2 = len(coords2)
				coords2 = coords2 + [ v2 ]
			faces2 = faces2 + [ [ iv0, iv1, iv2 ] ]
	return coords2, faces2


#
# platonic solids
# coordinates and faces from http://people.sc.fsu.edu/~burkardt/c_src/plato_ply/plato_ply.c
#

def tetrahedronData():
	coords = [ [-1, -1, -1], [1,  1, -1], [1, -1,  1], [-1,  1,  1] ]
	faces = [ [3, 2, 1], [2, 0, 1], [0, 2, 3], [3, 1, 0] ]
	return coords, faces

def octahedronData():
	coords = [ [1,  0,  0], [0, -1,  0], [-1,  0,  0], [0,  1,  0], [0,  0,  1], [0,  0, -1] ]
	faces = [ [1, 0, 4], [2, 1, 4], [3, 2, 4], [0, 3, 4], [0, 1, 5], [1, 2, 5], [2, 3, 5], [3, 0, 5] ]
	return coords, faces

def cubeData():
	coords = [ [-1, -1, -1], [1, -1, -1], [1,  1, -1], [-1,  1, -1], [-1, -1,  1], [1, -1,  1], [1,  1,  1], [-1,  1,  1] ]
	faces = [ [3, 2, 1, 0], [6, 7, 4, 5], [5, 1, 2, 6], [0, 4, 7, 3], [6, 2, 3, 7], [4, 0, 1, 5] ]
	return coords, faces

def icosahedronData():
	a = (sqrt(5) + 1) / 2
	b = 1
	len = sqrt(a*a + b*b)
	a = a / len
	b = b / len
	coords = [ [0, -b,  a], [a,  0,  b], [a,  0, -b], [-a,  0, -b], [-a,  0,  b], [-b,  a,  0], [b,  a,  0], [b, -a,  0], [-b, -a,  0], [0, -b, -a], [0,  b, -a], [0,  b,  a] ]
	faces = [ [1, 2, 6], [1, 7, 2], [3, 4, 5], [4, 3, 8], [6, 5, 11], [5, 6, 10], [9, 10, 2], [10, 9, 3], [7, 8, 9], [8, 7, 0], [11, 0, 1], [0, 11, 4], [6, 2, 10], [1, 6, 11], [3, 5, 10], [5, 4, 11], [2, 7, 9], [7, 1, 0], [3, 9, 8], [4, 8, 0] ]
	return coords, faces

def dodecahedronData():
	a = 1 + (1 + sqrt(5)) / 2
	b = 1
	c = (1 + sqrt(5)) / 2
	len = sqrt(a*a + b*b)
	a = a / len
	b = b / len
	c = c / len
	coords = [ [-c, -c, c], [a, b, 0], [a, -b, 0], [-a, b, 0], [-a, -b, 0], [0, a, b], [0, a, -b], [b, 0, -a], [-b, 0, -a], [0, -a, -b], [0, -a, b], [b, 0, a], [-b, 0, a], [c, c, -c], [c, c, c], [-c, c, -c], [-c, c, c], [c, -c, -c], [c, -c, c], [-c, -c, -c] ]
	faces = [ [14, 11, 18, 2, 1], [2, 17, 7, 13, 1], [15, 8, 19, 4, 3], [4, 0, 12, 16, 3], [16, 5, 6, 15, 3], [13, 6, 5, 14, 1], [18, 10, 9, 17, 2], [19, 9, 10, 0, 4], [17, 9, 19, 8, 7], [13, 7, 8, 15, 6], [16, 12, 11, 14, 5], [18, 11, 12, 0, 10] ]
	return coords, faces



# create a icosphere
# radius: radius of the sphere
# subdivisions: number of triangle subdivisions
def createIcosphere(radius, subdivisions):
	coords, faces = icosahedronData()
	q = angleBetween(coords[0], [0, 0, 1])
	coords = rotateArray(coords, q)
	coords = scaleArrayToLength(coords, radius)
	for i in range(subdivisions):
		coords, faces = subdivide(coords, faces, radius)
	return coords, faces

def VDrawMesh(coords, faces, color):
    lc = len(coords)
    carray = array(lc*[lc*[0]])
    #print carray
    for coord in coords:
        sphere(pos=coord, radius=0.1,color=(1,0,0))

    for face in faces:
        lp = len(face)
        for i in range(-1,lp):
            p1 = face[i]
            j = i+1
            if j == lp: j = -1
            p2 = face[j]
            if carray[p1][p2] == 1 or carray[p2][p1] == 1: continue
            wire = cylinder(pos=coords[p1], axis=vector(coords[p2])-vector(coords[p1]), 
                            radius=0.05, color=color)
            carray[p1][p2] = 1
            carray[p2][p1] = 1

# create platonic solids
solidRadius = 2
coords, faces = cubeData()
coords = scaleArrayToLength(coords, solidRadius)
coords = translateArray(coords, [ -5, 0, 0 ])

vcube = VDrawMesh(coords, faces, (0,0,1))    
    
coords, faces = octahedronData()
coords = scaleArrayToLength(coords, solidRadius)
coords = translateArray(coords, [ 0, 0, 0 ])

voctahedron = VDrawMesh(coords, faces, (0,1,0))    

coords, faces = dodecahedronData()
coords = scaleArrayToLength(coords, solidRadius)
coords = translateArray(coords, [ 5, 0, 0 ])

vdodecahedron = VDrawMesh(coords, faces, (1,1,1))    

coords, faces = icosahedronData()
coords = scaleArrayToLength(coords, solidRadius)
coords = translateArray(coords, [ 0, 5, 0 ])

vicosahedron = VDrawMesh(coords, faces, (0,1,1))    

coords, faces = tetrahedronData()
coords = scaleArrayToLength(coords, solidRadius)
coords = translateArray(coords, [ 5, 5, 0 ])

vtetrahedron = VDrawMesh(coords, faces, (1,0,1))    


# create geodesic dome
domeRadius = 7
domeSubdivisions = 3

coords, faces = createIcosphere(domeRadius, domeSubdivisions)
vfull = VDrawMesh(coords, faces, (1,1,0))