Sierpinsky-rotating-cubes.py

import turtle
import math
import random
import sys

# point is (x,y,z) float cooords
# polygon is ([points], color) is oriented for vector product to know if facing
# object is [polygons] is convex (cf. orientation)
# rotating_object is (object, (a,(u,v,w)) angle and axis of rotation (center of rotation is barycenter of object)
# figure is [rotating_objects]

def dot_product(v0, v1):
  (x0, y0, z0)= v0
  (x1, y1, z1)= v1
  return x0*x1 + y0*y1 + z0*z1

def rotate_face(polygon, angle, direction, center=(0,0,0)):
  res=[]
  #http://inside.mines.edu/fs_home/gmurray/ArbitraryAxisRotation/
  (cx, cy, cz)= center
  (u, v, w)= direction
  cos= math.cos(angle)
  sin= math.sin(angle)
  x1= (1-cos)*(cx*(v*v + w*w) - u*(cy*v + cz*w)) + sin*(cy*w - cz*v)
  x2= (1-cos)*u
  y1= (1-cos)*(cy*(w*w + u*u) - v*(cz*w + cx*u)) + sin*(cz*u - cx*w)
  y2= (1-cos)*v
  z1= (1-cos)*(cz*(u*u + v*v) - w*(cx*u + cy*v)) + sin*(cx*v - cy*u)
  z2= (1-cos)*w
  for p in polygon:
    (x,y,z)= p
    dotp= dot_product(p, direction)
    res.append((x1 + x2*dotp + x*cos + sin*(v*z - w*y)
               ,y1 + y2*dotp + y*cos + sin*(w*x - u*z)
               ,z1 + z2*dotp + z*cos + sin*(u*y - v*x)))
  return res

def rotate_object(obj, angle, direction, center=(0,0,0)):
  res= []
  for (polygon, color) in obj:
    res.append((rotate_face(polygon, angle, direction, center), color))
  return res

def rotate_figure(figure, angle, axis):
  res= []
  for (obj, rotation_data) in figure:
    res.append((rotate_object(obj, angle, axis), rotation_data))
  return res

def rotate_in_figure(figure):
  res= []
  for (obj, rotation_data) in figure:
    (angle, axis)= rotation_data
    poly_bary=[]
    for (polygon, color) in obj:
      poly_bary.append(barycenter(polygon))
    res.append((rotate_object(obj, angle, axis, barycenter(poly_bary)), rotation_data))
  return res

def substract_vectors(v0, v1):
  (x0, y0, z0)= v0
  (x1, y1, z1)= v1
  return (x0-x1, y0-y1, z0-z1)

def apply_perspective(points, camera, viewer):
  res=[]
  #http://en.wikipedia.org/wiki/3D_projection#Perspective_projection
  (vx, vy, vz)= viewer
  for p in points:
    (x, y, z) = substract_vectors(p, camera)
    res.append((vy-(vx*y)/x, vz-(vx*z)/x))
  return res

def make_rectangular_cuboid_vertices(center, dimensions):
  (cx, cy, cz)= center
  (dx, dy, dz)= dimensions
  res=[]
  for deltaX in (-dx/2, dx/2):
    for deltaY in (-dy/2, dy/2):
      for deltaZ in (-dz/2, dz/2):
        res.append((cx+deltaX, cy+deltaY, cz+deltaZ))
  return res;

def make_rectangular_cuboid_faces(center, dimensions, colors):
  points= make_rectangular_cuboid_vertices(center, dimensions)
  res=[]
  for (face_points, c) in zip(((0,2,3,1),(2,6,7,3),(0,1,5,4),(0,4,6,2),(1,3,7,5),(7,6,4,5)), colors):
    face=[] # list comprehension would be better
    for idx in face_points:
      face.append(points[idx])
    res.append((face,c))
  return res

def repeat_call(f,n):
  res=[]
  for i in range(n):
    res.append(f())
  return res

def add_vectors(v0, v1):
  (x0, y0, z0)= v0
  (x1, y1, z1)= v1
  return (x0+x1, y0+y1, z0+z1)

def custom_rotation(_, size, detail, indice):
  axis=[(1,0,0), (0,1,0), (0,0,1)]
  sign= 1
  if (indice == 0) or (indice == 5) or (indice == 6) or (indice == 3):
    sign= -1
  return (sign * math.pi/size, axis[detail%len(axis)])

def make_fractal_cuboids(center, size, color_function, rotation_function, detail, indice=0):
  size/= 3
  res= [(make_rectangular_cuboid_faces(center, (size, size, size), repeat_call(color_function, 6))
         ,rotation_function(center, size, detail, indice))]
  if detail > 0:
    indice= 0
    d= (-size, size)
    for dx in d:
      for dy in d:
        for dz in d:
          res+= make_fractal_cuboids(add_vectors(center,(dx, dy, dz)), size
                                     , color_function, rotation_function, detail-1, indice)
          indice+=1
  return res

def draw_face_2D(points, color):
  turtle.up()
  turtle.goto(points[0])
  turtle.fillcolor(color)
  turtle.down()
  turtle.begin_fill()
  for p in points:
    turtle.goto(p)
  turtle.end_fill()

def cross_product(v0, v1):
  (x0, y0, z0)= v0
  (x1, y1, z1)= v1
  return (y0*z1 - z0*y1, z0*x1 - x0*z1, x0*y1 - y0*x1)

def points_to_vector(p0, p1):
  (x0, y0, z0)= p0
  (x1, y1, z1)= p1
  return (x1-x0, y1-y0, z1-z0)

def draw_colored_faces(colored_faces, camera, viewer):
  for (face_points, color) in colored_faces:
    if dot_product(cross_product(points_to_vector(face_points[1], face_points[0])
                                ,points_to_vector(face_points[1], face_points[2]))
                   ,points_to_vector(camera, viewer)) >0 :
      draw_face_2D(apply_perspective(face_points, camera, viewer), color)

def barycenter(points):
  (cx, cy, cz)= (0., 0., 0.)
  for (x,y,z) in points:
    cx+= x
    cy+= y
    cz+= z
  return (cx/len(points), cy/len(points), cz/len(points))

def dist_square(p0, p1):
  (x0, y0, z0)= p0
  (x1, y1, z1)= p1
  return (x1-x0)*(x1-x0)+(y1-y0)*(y1-y0)+(z1-z0)*(z1-z0)

def make_dist_key(viewer):
  def dist_key(colored_face):
    (points, color)= colored_face
    return dist_square(barycenter(points), viewer)
  return dist_key

def get_colored_faces_from_figure(figure):
  res=[]
  for (obj, _) in figure:
    res+= obj
  return res

def draw_colored_oriented_faces_viewer_order(colored_faces, camera, viewer):
  colored_faces.sort(key= make_dist_key(viewer), reverse= True)
  for (face_points, color) in colored_faces:
    if dot_product(cross_product(points_to_vector(face_points[1], face_points[0])
                                ,points_to_vector(face_points[1], face_points[2]))
                   ,points_to_vector(camera, viewer)) >0 :
      draw_face_2D(apply_perspective(face_points, camera, viewer), color)

def draw(figure, camera, viewer):
  draw_colored_oriented_faces_viewer_order(get_colored_faces_from_figure(figure), camera, viewer)

def random_color():
  return (random.random(), random.random(), random.random())

def colorize_faces(faces, color_function):
  res=[]
  for f in faces:
    res.append((f, color_function()))
  return res

def scale_vector(v, s):
  (x, y, z)= v
  return (x*s, y*s, z*s)

def scale_figure(figure, s):
  res= []
  for (obj, rotation_data) in figure:
    scaled_obj= []
    for (points, color) in obj:
      scaled_points= []
      for point in points:
        scaled_points.append(scale_vector(point, s))
      scaled_obj.append((scaled_points, color))
    res.append((scaled_obj, rotation_data))
  return res

# for i in test_prefix-2-00*.ps; do convert $i ${i%.ps}.png; done
# convert -background white -alpha remove test_prefix-2-00*.png test_prefix.gif
def pic_to_file(filename):
  turtle.getscreen().getcanvas().postscript(file=filename)


def demo(details=2, outfile_prefix=None):
  turtle.title("3d cube demo")
  turtle.hideturtle()
  cuboids= make_fractal_cuboids((10, 10, 10), 200, random_color, custom_rotation, details)
  turtle.tracer(0)
  axis = 1, 1, 1
  camera = 500, 0, 0
  viewer = 550, 0, 0
  angle = 0.01
  scale_delta= 1.02
  scale= 1
  index=0
  while True:
    turtle.clear()
    draw(scale_figure(cuboids, scale), camera, viewer)
    turtle.update()
    if outfile_prefix:
      pic_to_file(outfile_prefix+"-"+str(details)+"-"+format(index, '04d')+".ps")
      index+= 1
    cuboids=rotate_figure(rotate_in_figure(cuboids), angle, axis)
    scale*= scale_delta
    if scale>1.5:
      scale= 1.5
      scale_delta= 0.9
    if scale < 0.5:
      scale= 0.5
      scale_delta= 1.1

if __name__ == "__main__":
  details= 2
  outfile_prefix= None
  if len(sys.argv) == 3:
    outfile_prefix= sys.argv[2]
  if len(sys.argv) >= 2:
    details= int(sys.argv[1])
  demo(details, outfile_prefix)