March Madness: Spherify!

I can’t just sit on the sidelines glowering at all the happy coders all month. Here’s my first code snippet for the month: spherify!

It’s a simple little python script that uses an image file as a heightmap, maps it on to a sphere, and spits out an STL file.  It can be charitably described as “crude”.  Usage:

Usage: spherify.py [options] source

Options:
  -h, --help            show this help message and exit
  -o OUTPATH, --output=OUTPATH
                        output to given path
  -r RADIUS, --radius=RADIUS
                        radius of lowest point in mm
  --height=HEIGHT       additional radius of highest point in mm

You’ll need to have the python imaging library installed to use it. Source after the break.

#!/usr/bin/env python

# Spherify is a simple script that takes in an
# image and uses it as a heightmap for the surface
# of a sphere.  The output is an STL file.

from PIL import Image
from optparse import OptionParser
import sys
from math import pi, sin, cos, sqrt

def f2s(f):
    "Convert a float to a string with given precision."
    # Really just to strip out sci. notation, which makes
    # STL cry
    return "%8.8f" % f

def normal(v1,v2,v3):
    "Compute an approximate normal for a tri"
    norm = [v1[i] + v2[i] + v3[i] for i in range(3)]
    l = sqrt(sum(map(lambda x:x**2,norm)))
    return tuple(map(lambda x:x/l,list(norm)))

def spherify(image,out,radius,height):
    "Write out the given image as a sphere in STL format"
    out.write("solid spherified\n")  # start object
    data = image.getdata()
    (min,max) = image.getextrema()
    (w,h) = image.size
    y_rpp = pi/(h+2)    # Y radians per pixel (adapted for poles)
    x_rpp = 2*pi/w      # X radians per pixel
    def to_vertex(x,y,v=None):
        "convert from (X,Y,value) to a 3-tuple coordinate"
        if v == None:
            v = data[(x%w) + (y-1)*w]
        phi = y_rpp*y
        theta = x_rpp*x
        r = radius + (height*v)/max
        c3 = r * cos(phi)
        a = r * sin(phi)
        c2 = a * -cos(theta)
        c1 = a * sin(theta)
        return (c1,c2,c3)
    def print_facet(v1,v2,v3):
        "print out a facet"
        n = normal(v1,v2,v3)
        out.write("  facet normal "+" ".join(map(f2s,n))+"\n")
        out.write("    outer loop\n")
        out.write("      vertex "+" ".join(map(f2s,v1))+"\n")
        out.write("      vertex "+" ".join(map(f2s,v2))+"\n")
        out.write("      vertex "+" ".join(map(f2s,v3))+"\n")
        out.write("    endloop\n")
        out.write("  endfacet\n")
    # top triangles
    top=to_vertex(0,0,data[0])
    for x in range(w):
        v1=top
        v2=to_vertex(x,1)
        v3=to_vertex(x+1,1)
        print_facet(v1,v2,v3)
    # middle triangles
    for y in range(h-1):
        for x in range(w):
            v1=to_vertex(x,y+1)
            v2=to_vertex(x+1,y+1)
            v3=to_vertex(x,y+2)
            v4=to_vertex(x+1,y+2)
            print_facet(v1,v2,v3)
            print_facet(v2,v3,v4)
    # bottom triangles
    bot=to_vertex(0,h+2,data[w*h-1])
    for x in range(w):
        v1=to_vertex(x,h+1,data[x+(h-1)*w])
        v2=to_vertex(x+1,h+1,data[((x+1)%w)+(h-1)*w])
        v3=bot
        print_facet(v1,v2,v3)    
    out.write("endsolid\n") # end object


def main():
    parser = OptionParser(usage="usage: %prog [options] source")
    parser.add_option("-o","--output",dest="outPath",
                      help="output to given path")
    parser.add_option("-r","--radius",dest="radius",
                      type="int", default="35",
                      help="radius of lowest point in mm")
    parser.add_option("--height",dest="height",
                      type="int", default="10",
                      help="additional radius of highest point in mm")
    (options,args) = parser.parse_args()
    if len(args) != 1:
        parser.error("Please provide a single input file.")
    image = Image.open(args[0])
    image = image.convert("L")
    outFile = sys.stdout
    if (options.outPath):
        outFile = open(options.outPath)
    spherify(image,outFile,options.radius,options.height)


if __name__ == "__main__":
    main()