"""
A simple waterman maker for Jython
Finds all IVM balls within a radius such that volume is maximized, convexity preserved
(i.e. generates classic Waterman polyhedra). Scales/colors for POV-Ray rendering.
Convex Hull finding added courtesy of John E. Lloyd, who implemented Qhull-style
algorithms in his QuickHull3D Java package -- and dependency for this program
(try same directory, as a folder named quickhull3d).
http://www.cs.ubc.ca/spider/lloyd/java/quickhull3d.html
Thanks also to Jim & Jython team, to Guido, Tim etc.
by Kirby Urner, June 16, 2005
Version 2.1
Note: From my point of view, Steve Waterman is the originator of this family of
IVM-based polyhedra -- convex hulls with a maximum number of CCP ball vertices
of a given distance from the origin or less. That's why I (without prompting
from Steve) named them Waterman polyhedra -- and the name appears to have stuck.
Usage:
Launch Jython
>>> import wpolys
>>> wpolys.waterman(10)
[ will write w10.pov to disk, ready for rendering with POV-Ray ]
I may not be able to provide any tech support. For the hardy.
Use at your own risk.
Note: no set object in this version of Jython, so I use dictionaries
with None for values -- a way to preserve the uniqueness of keys.
"""
import quickhull3d as qh
def waterman(n):
"""
collect all ccp balls within a given sweepout radius
"""
filename = 'w'+str(n)+'.pov'
fileout = open(filename,'w')
fileout.write(header % (10 * pow(n,0.5)))
allballs = []
oldverts = []
# find balls from the outside in, 6 layers at a time,
# stop when convex hull registers no change in vertices
print 'Finding ivm lattice points'
layer = 0
for r in range(2*n,0,-2):
print '.',
balls = finder(r)
for b in balls:
newballs = permuter(*b)
allballs.extend( newballs )
layer += 1
if layer == 6 or r<=2:
hull = qh.QuickHull3D()
points = [qh.Point3d(*ball) for ball in allballs]
hull.build(points)
verts = [(p.x, p.y, p.z) for p in hull.getVertices()]
newverts = verts[:]
newverts.sort()
print '|',
if newverts == oldverts or r<=2:
break # convex hull is stable and/or all points found
else:
oldverts = newverts
layer = 0
print "\nConvex hull found\n"
faces = [list(f) for f in hull.getFaces()]
edges = getedges(faces)
povarray(n, verts, fileout)
for face in faces:
povface(face, fileout)
for edge in edges:
povedge(edge, fileout)
fileout.close()
def finder(r):
"""
Explore first octant for unique balls
"""
# main loop
okballs = []
maxx = int(pow(r,0.5))
for x in range(maxx,0,-1): # x drops from r to 1
for y in range(0,x+1): # y up to x
if x**2 + y**2 > r: # quit if too big
break
for z in range(0,y+1): # z up to y
if not (x+y+z) % 2: # even sum only
if x**2 + y**2 + z**2 == r:
okballs.append((x,y,z)) # got one
return okballs
def permuter(x,y,z):
xyz = (x,y,z)
unique = {} # import sets, use sets.Set() if < Python 2.4
# revert to dictionary for Jython (approx Python 2.1 as of now)
unique[(xyz[0],xyz[1],xyz[2])] = None
unique[(xyz[0],xyz[2],xyz[1])] = None
unique[(xyz[1],xyz[0],xyz[2])] = None
unique[(xyz[1],xyz[2],xyz[0])] = None
unique[(xyz[2],xyz[0],xyz[1])] = None
unique[(xyz[2],xyz[1],xyz[0])] = None
signs = [('+','+','-'),
('+','-','+'),
('-','+','+'),
('+','-','-'),
('-','+','-'),
('-','-','+'),
('-','-','-')]
for combo in unique.keys():
for pattern in signs:
newcombo = [eval(i + str(j))
for i,j in zip(pattern, combo)]
unique[tuple(newcombo)] = None
return unique.keys()
def getedges(faces):
alledges = {}
for f in faces:
for a,b in zip(f, f[1:]+[f[0]]):
candidate = [a,b]
candidate.sort()
alledges[tuple(candidate)] = None
return alledges.keys()
def povarray(n, verts, fileout, tex='texture6'):
"""
cut and paste into .pov file with standard (or custom) header
"""
scale = 1 # pow(2,0.5)*n
fileout.write(arrayskel % len(verts))
for b in verts[:-1]:
b = tuple([x/scale for x in b])
theline = " <%d,%d,%d>,\n" % b
fileout.write(theline)
b = verts[-1]
b = tuple([x/scale for x in b])
fileout.write(" <%d,%d,%d>\n" % b)
fileout.write("}\n")
fileout.write( arrayloop % (0.04, 'texture1') )
def povface(face, fileout):
fileout.write("polygon{%s,\n" % (len(face)+1))
for p in face:
fileout.write("Pts[%s] " % p)
fileout.write("Pts[%s] " % face[0])
fileout.write("\ntexture {texture0}\n }\n")
def povedge(edge, fileout):
fileout.write("cylinder{ Pts[%s], Pts[%s], 0.04 pigment {texture1}}\n" % (edge[0], edge[1]))
header= """//POV-Ray script
//version 3.6
global_settings { assumed_gamma 2.2 }
#include "colors.inc"
#include "metals.inc"
#include "glass.inc"
#declare Cam_factor = %s;
#declare Camera_X = 0.5 * Cam_factor;
#declare Camera_Y = 0.5 * Cam_factor;
#declare Camera_Z = 1 * Cam_factor;
#declare texture0 = T_Copper_2A;
#declare texture1 = Col_Dark_Green;
camera { location
up <0, 1.0, 0> right <-4/3, 0, 0>
direction <0, 0, 3> look_at <0, 0, 0>
rotate <0,0,0>}
light_source { color White shadowless}
light_source { color White shadowless}
light_source { color White shadowless}
//Background:
background {color Black}
"""
arrayskel = """
# declare NumPoints = %s;
# declare Pts = array [NumPoints] {
"""
arrayloop = """
#declare ptno=0;
#while(ptno