#!/usr/bin/env python
#----------
#
# Placed into the public domain in November 2008 by Bob Harris.
#
# solution_to_pov_pieces--
#    Convert a packed polyomino solution to pieces that can be pasted into a
#    a povray file.
#
# The input comes in two forms.  In both formats, a "-" indicates an empty cell.
#
# Note that both files can contain 'directives' which are indicated by the first
# character being a ">".  Directives are not currently documented-- for their
# meaning you need to look at do_directive.
#
# The default format is like this (y,z labels are not part of the input file):
#
#    ka k4 p9 rq rq rq m1   (y = 6)  (z = 0)
#    ka k4 p9 rq rq m1 pf
#    ka ka ka k2 k2 nz pf
#    ni ni k2 k2 k2 nz nz
#    ls p2 p2 p2 mb k5 nz
#    ls mo mo mo lz k5 mb
#    ls mo mo lz lz k5 k5   (y = 0)
#
#    ka k4 k4 k4 m1 m1 m1            (z = 1)
#    ka k4 p9 p9 rq m1 m1
#    ni p9 p9 p9 pf pf pf
#    ni ni ni k2 pf pf nz
#    ls ls p2 p2 mb mb nz
#    ls ls p2 lz mb mb mb
#    mo mo lz lz lz k5 k5
#
#     ...
#
#    lk lk lm ly ly ja ja            (z = 6)
#    lk lk lm ly ly lu lu
#    lp jw jw jw ly lu js
#    lp jw jw m9 m9 m9 js
#    lt lt nu m9 jz js js
#    lt nu nu lo jz js oq
#    lt nu nu lo jz jz oq
#
# The other input format (--z-across) has a single letter per piece, like
# this 49 piece packing of a 7x7x7:
#
#    gdfkXFF  ddffXFF  deefXXX  deeeSSL  BMMtSLL  BBMtSJJ  BBMtJJJ   (y = 6)
#    ggkkXFF  dgkfXbb  dgkfbbb  eellZLL  BHMVSSL  BMMtSIL  oottJIJ
#    mgqqKFK  mgkfAAb  llkiAAA  llliZnn  HHVVZZZ  HHsVDII  ootVIII
#    mvqqKKK  mvqiKKb  CCiiAAc  CCCiZnn  HHVVZRR  oossDRR  ojjsDDR
#    mmpqTTT  vvqpQQT  vvihYcc  vCChEnU  uuusEnn  juusDRG  jjjsDDR
#    wmppQcT  wrhpQcT  wrhhYcc  OONhEEU  OaNNEEU  jauWPEG  aauWPPG
#    rrrpQQT  wrrpQYY  wwNhYYY  OwNNUUU  OONWPPU  OWWWPGG  aaaWPGG   (y = 0)
#    (z=0)                                                 (z=6)
#
# The output is a seies of piece definitions, like this one:
#
#    // piece q
#    #declare piece43 = array[7] { <1,0,0>, <0,1,0>, <1,1,0>, <0,2,0>, ... }
#    #declare loc43 = <2,2,0>;             (position of bounding box)
#    #declare com43 = <13,21,11>/14.0;     (center of mass)
#    #declare col43 = <.02,.89,.40>;       (random color)
#
#----------

import sys,re,random

debug = []

cameraDist  = None
cameraTrans = None
imageAspect = None
drawHidden  = None
useClock    = None
softColors  = None
showAxes    = None
showBox     = None

xMax = yMax = zMax = numPieces = 0
pieceNames = None

orderHints  = {}
givenColors = []

#----------
#
# solution_to_pov_pieces--
#    main program
#
#----------

def usage(s=None):
    message = """
solution_to_pov_pieces
  <name>.<ext>            input file (a polyomino packing file)
                          output will be in <name>.ini and <name>.pov
  --cameradist=<num>      set the distance from the camera to the puzzle
  --cameratrans=<x,y,z>   set camera translation vector
  --aspect=<w/h>          set image aspect ratio
  --z-across              z-axis goes left to right in input file
                          (default is top to bottom)
  --drawhidden            "draw" hidden pieces (slows down rendering)
  --colors=<num>          allow at least <num> colors
  --colors=all            every piece gets a different color
  --softcolors            choose colors that are generally softer
  --neighborhood=[18,26]  specify neighborhood size for different piece colors
                          (default is 6)
  --noclock               create a single-frame pov with all pieces placed
  --axes                  show x,y,z axes
  --box                   show bounding box
  --seed=<seed>           specify random number seed
"""
    if (s == None): sys.exit (message)
    else:           sys.exit ("%s\n%s" % (s,message))


def main():
    global debug
    global cameraDist,cameraTrans,imageAspect
    global drawHidden,useClock,softColors,showAxes,showBox
    global xMax,yMax,zMax,numPieces
    global pieceNames

    ##########
    # parse the command line
    ##########

    zAcross      = False
    cameraDist   = None
    cameraTrans  = None
    imageAspect  = None
    seed         = None
    neighborhood = 6
    numColors    = None
    useClock     = True
    useClock     = True
    drawHidden   = False
    showAxes     = False
    showBox      = False
    inputName    = None

    # pick off options

    args = sys.argv[1:]
    while (len(args) > 0):
        arg = args.pop(0)
        val = None
        fields = arg.split("=",1)
        if (len(fields) == 2):
            arg = fields[0]
            val = fields[1]
            if (val == ""):
                usage("missing a value in %s=" % arg)

        if (arg == "--z-across") and (val == None):
            zAcross = True
        elif (arg == "--cameradist") and (val != None):
            cameraDist = val
        elif (arg == "--cameratrans") and (val != None):
            cameraTrans = val.split(",")
        elif (arg == "--aspect") and (val != None):
            imageAspect = val
        elif (arg == "--drawhidden") and (val == None):
            drawHidden = True
        elif (arg == "--neighborhood") and (val == "6"):
            neighborhood = 6
        elif (arg == "--neighborhood") and (val == "18"):
            neighborhood = 18
        elif (arg == "--neighborhood") and (val == "26"):
            neighborhood = 26
        elif (arg == "--colors") and (val == "all"):
            numColors = "all"
        elif (arg == "--colors") and (val != None):
            numColors = int(val)
        elif (arg == "--noclock") and (val == None):
            useClock = False
        elif (arg == "--softcolors") and (val == None):
            softColors = True
        elif (arg == "--axes") and (val == None):
            showAxes = True
        elif (arg == "--box") and (val == None):
            showBox = True
        elif (arg == "--seed") and (val != None):
            try:    seed = int(val)
            except: seed = val
        elif (arg == "--debug") and (val != None):
            debug.append(val)
        elif (arg == "--debug") and (val == None):
            debug.append("debug")
        elif (arg.startswith("--")) and (val == None):
            usage("unknown argument: %s" % arg)
        elif (inputName == None) and (val == None):
            inputName = arg
        else:
            usage("unknown argument: %s" % arg)

    if (inputName == None):
        usage("you must provide an input file name")

    dot = inputName.rfind(".")
    if (dot == -1): outputName = inputName
    else:           outputName = inputName[:dot]

    if (seed != None):
        random.seed(seed)

    ##########
    # read the solution
    ##########

    f = file(inputName,"rt")
    if (zAcross): cells = read_across_packing(f)
    else:         cells = read_packing(f)
    f.close()

    xMax = len(cells)
    yMax = len(cells[0])
    zMax = len(cells[0][0])

    components = connected_components(cells)

    pieceNames = components.keys()
    pieceNames.sort()
    numPieces = len(pieceNames)

    # write the ini file

    f = file(outputName+".ini","wt")
    print_ini(f,outputName+".pov")
    f.close()

    # start creating the pov file

    f = file(outputName+".pov","wt")
    print_pov_main(f)

    print >>f,"// packing:"
    print >>f,"//"
    if (zAcross): print_cells_across(f,"//  ",cells)
    else:         print_cells       (f,"//  ",cells)

    ##########
    # assign colors
    ##########

    # assign colors

    if (numColors == "all"):
        coloring = {}
        for (ix,piece) in enumerate(pieceNames):
            coloring[piece] = ix+1
        numColors = len(pieceNames)
    else:
        g = neighbor_graph(cells,components,neighborhood=26)
        (numColors,coloring) = color_graph(g,numColors)
        print >>f,"// coloring graph:"
        print >>f,"//\t(%d colors)" % numColors
        print_graph_and_colors(f,"//   ", g,coloring)
        print >>f

    # create colors

    colors = generate_color_set(numColors)

    # print the color set

    for color in range(numColors):
        (r,g,b) = colors[color]
        print >>f,"#declare col%d = %s;" % (color+1,cubie_string(r,g,b))
    print >>f

    ##########
    # write the piece descriptions
    ##########

    nameToNum = {}
    numToName = {}

    for (num,name) in enumerate(pieceNames):
        nameToNum[name] = num+1
        numToName[nameToNum[name]] = name
        cubies = [(x,y,z) for (x,y,z) in components[name]]
        n = len(cubies)
        xx = min([x for (x,y,z) in cubies])    # (position of bounding box)
        yy = min([y for (x,y,z) in cubies])
        zz = min([z for (x,y,z) in cubies])
        (cx,cy,cz) = (0,0,0)                # (will be center of mass)
        for ix in range(n):
            (x,y,z) = cubies[ix]
            (x,y,z) = (x-xx,y-yy,z-zz)
            cubies[ix] = (x,y,z)
            (cx,cy,cz) = (cx+x,cy+y,cz+z)

        (cx,cy,cz) = (2*cx+n,2*cy+n,2*cz+n)    # (add 1/2 in each direction)

        print >>f,"// piece %s" % name
        print >>f,"#declare piece%d = array[%d] { %s }" % (nameToNum[name],n,cubies_string(cubies))
        print >>f,"#declare loc%d = %s;"                % (nameToNum[name],cubie_string(xx,yy,zz))
        print >>f,"#declare com%d = %s/%d.0;"           % (nameToNum[name],cubie_string(cx,cy,cz),2*n)
        print >>f

    ##########
    # write the piece placements
    ##########

    # determine the placement order

    order = piece_order(cells,components)

    # determine when pieces become completely hidden by subsequent placements

    if (not drawHidden):
        cellOccupied = zero_array(xMax,yMax,zMax)
        for (ix,name) in enumerate(order):
            for (x,y,z) in components[name]:
                cellOccupied[x][y][z] = ix+1

        hiddenTime = {}
        for name in pieceNames:
            hiders = hiding_cells(cells,components[name])
            if (hiders == None): continue
            hiddenTime[name] = max([cellOccupied[x][y][z] for (x,y,z) in hiders])

    for (ix,name) in enumerate(order):
        num = nameToNum[name]

        if (useClock):
            if (drawHidden) or (name not in hiddenTime):
                print >>f,"#if (clock >= %d)" % (ix+1)
            else:
                print >>f,"#if (clock >= %d & clock < %d)" \
                        % (ix+1,hiddenTime[name])

        print >>f,"object { make_piece(piece%d,col%d) translate loc%d }" \
            % (num,coloring[name],num)
        if (useClock):
            print >>f,"#end"
        print >>f

    f.close()

##########
# read a puzzle description (see this file's header for formats)
##########

def read_packing(f):
    planes = []
    plane  = []
    for line in f:
        line = line.rstrip()
        if (line.startswith(">")):
            do_directive(line)
        elif (line != ""):
            plane.insert(0,line.split())
        elif (plane != []):
            planes.append(plane)
            plane = []
            continue

    if (plane != []):
        planes.append(plane)

    return [[[planes[z][y][x] for z in range(len(planes))]
                              for y in range(len(planes[0]))]
                              for x in range(len(planes[0][0]))]


def read_across_packing(f):
    cells = []
    for line in f:
        if (line.startswith(">")):
            do_directive(line)
        elif (line != ""):
            line = line.rstrip().split()
            cells.insert(0,line)

    return [[[cells[y][z][x] for z in range(len(cells[0]))]
                             for y in range(len(cells))]
                             for x in range(len(cells[0][0]))]


def do_directive(line):
    line = line[1:].strip()
    if (do_color_directive(line)): return
    if (do_order_directive(line)): return
    assert False, "unknown directive: \"%s\"" % line


# handle a color directive

colorPatt = \
      '^color' \
    + ' +(?P<red>[0-9]*\.[0-9]*)' \
    + ' *, *(?P<green>[0-9]*\.[0-9]*)' \
    + ' *, *(?P<blue>[0-9]*\.[0-9]*)' \
    + ' *' \
    + '$'
colorRe = re.compile(colorPatt)

def do_color_directive(s):
    global givenColors

    m = colorRe.match(s)
    if (m == None): return False

    red   = float(m.group('red'))
    green = float(m.group('green'))
    blue  = float(m.group('blue'))

    assert (0 <= red   <= 1)
    assert (0 <= green <= 1)
    assert (0 <= blue  <= 1)

    givenColors.append((red,green,blue))
    return True


# handle an order directive

orderPatt = \
      '^order' \
    + ' *(?P<order>[0-9]+)' \
    + ' *\(' \
    + ' *(?P<pieces>.+)' \
    + ' *\)' \
    + '$'
orderRe = re.compile(orderPatt)

def do_order_directive(s):
    global orderHints

    m = orderRe.match(s)
    if (m == None): return False

    order  = m.group('order')
    pieces = m.group('pieces')

    order  = int(m.group('order'))
    pieces = pieces.split()

    for name in pieces:
        orderHints[name] = order

    return True

##########
# deal with puzzle "pieces"
##########

# returns a hash from piece name to a list of (x,y,z) cubies

def connected_components(cells):
    components = {}

    for z in range(zMax):
        for y in range(yMax):
            for x in range(xMax):
                name = cells[x][y][z]
                if (name == "-"): continue
                if (name not in components): components[name] = []
                components[name].append((x,y,z))

    for name in components.keys():
        cubies = [(z,y,x) for (x,y,z) in components[name]]
        cubies.sort()
        components[name] = [(x,y,z) for (z,y,x) in cubies]

    return components

# returns a hash from piece name to a list of neighboring piece names
# neighborhood is one of 6,18,26 depending on the size of neighborhood desired

def neighbor_graph(cells,components,neighborhood=6):
    g = {}
    for piece in components.keys():
        g[piece] = []

    for piece in components.keys():
        ngbrs = g[piece]
        for (x,y,z) in components[piece]:
            if (x > 0):
                ngbr = cells[x-1][y][z]
                if (ngbr not in [piece,"-"]) and (ngbr not in ngbrs):
                    ngbrs.append(ngbr)
                    g[ngbr].append(piece)
            if (y > 0):
                ngbr = cells[x][y-1][z]
                if (ngbr not in [piece,"-"]) and (ngbr not in ngbrs):
                    ngbrs.append(ngbr)
                    g[ngbr].append(piece)
            if (z > 0):
                ngbr = cells[x][y][z-1]
                if (ngbr not in [piece,"-"]) and (ngbr not in ngbrs):
                    ngbrs.append(ngbr)
                    g[ngbr].append(piece)

            if (neighborhood < 18): continue
            if (x > 0) and (y > 0):
                ngbr = cells[x-1][y-1][z]
                if (ngbr not in [piece,"-"]) and (ngbr not in ngbrs):
                    ngbrs.append(ngbr)
                    g[ngbr].append(piece)
            if (x < xMax-1) and (y > 0):
                ngbr = cells[x+1][y-1][z]
                if (ngbr not in [piece,"-"]) and (ngbr not in ngbrs):
                    ngbrs.append(ngbr)
                    g[ngbr].append(piece)
            if (x > 0) and (z > 0):
                ngbr = cells[x-1][y][z-1]
                if (ngbr not in [piece,"-"]) and (ngbr not in ngbrs):
                    ngbrs.append(ngbr)
                    g[ngbr].append(piece)
            if (x < xMax-1) and (z > 0):
                ngbr = cells[x+1][y][z-1]
                if (ngbr not in [piece,"-"]) and (ngbr not in ngbrs):
                    ngbrs.append(ngbr)
                    g[ngbr].append(piece)
            if (y > 0) and (z > 0):
                ngbr = cells[x][y-1][z-1]
                if (ngbr not in [piece,"-"]) and (ngbr not in ngbrs):
                    ngbrs.append(ngbr)
                    g[ngbr].append(piece)
            if (y < yMax-1) and (z > 0):
                ngbr = cells[x][y+1][z-1]
                if (ngbr not in [piece,"-"]) and (ngbr not in ngbrs):
                    ngbrs.append(ngbr)
                    g[ngbr].append(piece)

            if (neighborhood < 26): continue
            if (x > 0) and (y > 0) and (z > 0):
                ngbr = cells[x-1][y-1][z-1]
                if (ngbr not in [piece,"-"]) and (ngbr not in ngbrs):
                    ngbrs.append(ngbr)
                    g[ngbr].append(piece)
            if (x < xMax-1) and (y > 0) and (z > 0):
                ngbr = cells[x+1][y-1][z-1]
                if (ngbr not in [piece,"-"]) and (ngbr not in ngbrs):
                    ngbrs.append(ngbr)
                    g[ngbr].append(piece)
            if (x > 0) and (y < yMax-1) and (z > 0):
                ngbr = cells[x-1][y+1][z-1]
                if (ngbr not in [piece,"-"]) and (ngbr not in ngbrs):
                    ngbrs.append(ngbr)
                    g[ngbr].append(piece)
            if (x < xMax-1) and (y < yMax-1) and (z > 0):
                ngbr = cells[x+1][y+1][z-1]
                if (ngbr not in [piece,"-"]) and (ngbr not in ngbrs):
                    ngbrs.append(ngbr)
                    g[ngbr].append(piece)

    return g

# returns the number of colors and a hash from piece name to a color index;
# this may not return a minimal coloring, but it is usually fine in practice

def color_graph(g,numColors=0):
    if (type(numColors) != type(0)) or (numColors < 0): numColors = 0

    coloring = {}

    # make a list of pieces, sorted by decreasing number of neighbors

    uncolored = [(-len(g[piece]),piece) for piece in g]
    uncolored.sort()

    # assign colors to pieces in order of most neighborly

    for (ngbrCount,piece) in uncolored:
        bad  = [coloring[ngbr] for ngbr in g[piece] if (ngbr in coloring)]
        good = [color for color in range(1,numColors+1) if (color not in bad)]
        if (good == []):
            numColors += 1
            color = numColors
        else:
            color = random.choice(good)
        coloring[piece] = color

    # permute the colors so that color index 1 is the most used, etc.

    colorUsed = {}
    for color in range(1,numColors+1):
        colorUsed[color] = 0

    for piece in g:
        color = coloring[piece]
        colorUsed[color] += 1

    usage = [(colorUsed[color],color) for color in range(1,numColors+1)]
    usage.sort()
    usage.reverse()

    oldToNew = {}
    for ix in range(len(usage)):
        newColor = ix+1
        (count,oldColor) = usage[ix]
        oldToNew[oldColor] = newColor
        # print "color %d used %d times" % (newColor,count)

    for piece in g:
        coloring[piece] = oldToNew[coloring[piece]]

    return (numColors,coloring)

# print a neighbor graph

def print_graph(f,prefix,g):
    if (prefix == None): prefix = ""

    names = g.keys()
    names.sort()

    for piece in names:
        ngbrs = [n for n in g[piece]]
        ngbrs.sort()
        print >>f,"%s%s: %s" % (prefix,piece," ".join(ngbrs))

# print a neighbor graph with color info

def print_graph_and_colors(f,prefix,g,coloring):
    if (prefix == None): prefix = ""

    names = g.keys()
    names.sort()

    for piece in names:
        color = coloring[piece]
        ngbrs = [n for n in g[piece]]
        ngbrs.sort()
        ngbrs = ["%s(%s)" % (n,coloring[n]) for n in ngbrs]
        print >>f,"%s%s(%s): %s" % (prefix,piece,color," ".join(ngbrs))


##########
# generate a set of colors more-or-less evenly spaced through the color space
#
# this could certainly be improved (e.g. compute distances in a better color
# space than rgb);  but it works fine in practice
##########

resol    = 1000
distance = None

def generate_color_set(numColors):
    global distance

    numLocked = len(givenColors)
    if (numLocked >= numColors):
        bestColors = []
        for ix in range(numColors):
            (r,g,b) = givenColors[ix]
            bestColors.append((r,g,b))
        return bestColors

    numLocked =  len(givenColors)
    colors    =  [(int(resol*r),int(resol*g),int(resol*b)) for (r,g,b) in givenColors]
    colors    += [random_color() for ix in range(numLocked,numColors)]
    numColors =  len(colors)

    # compute all distances

    distance    = zero_array(numColors,numColors)
    distanceRow = zero_array(numColors)

    for ix in range(numColors-1):
        (r,g,b) = colors[ix]
        for iy in range(ix+1,numColors):
            d = color_dist((r,g,b),colors[iy])
            distance[ix][iy] = d

    # try to improve the set

    bestD = 0
    bestColors = [(r,g,b) for (r,g,b) in colors]

    newColors = 0
    while (newColors < 100):
        # select the closest pair of colors

        (d,ix,iy) = closest_pair(numLocked)
        if (ix >= numLocked):
            if (total_distance(iy) > total_distance(ix)): (ix,iy) = (iy,ix)

        # choose a replacement color for colors[iy]

        (r,g,b) = random_color()

        # compute distances to it

        for ix in range(numColors):
            d = color_dist((r,g,b),colors[ix])
            distanceRow[ix] = d

        minD = min([distanceRow[ix] for ix in range(numColors) if (ix!=iy)])

        # if it's an improvement, keep it;  otherwise, randomly keep it

        if (minD < d):
            keepProb = float(minD)/d
            if (random.random >= keepProb): continue

        newColors += 1
        colors[iy] = (r,g,b)
        for ix in range(iy):             distance[ix][iy] = distanceRow[ix]
        for ix in range(iy+1,numColors): distance[iy][ix] = distanceRow[ix]

        if (d > bestD):
            bestD = 0
            bestColors = [(r,g,b) for (r,g,b) in colors]

    for ix in range(numColors):
        (r,g,b) = bestColors[ix]
        bestColors[ix] = (float(r)/resol,float(g)/resol,float(b)/resol)

    return bestColors

# find the closest pair of colors in the set, with the restiction that only the
# first member of the pair can be a 'locked' color

def closest_pair(numLocked):
    n = len(distance)
    dists = []
    for ix in range(n-1):
        rest = max(ix+1,numLocked)
        dists.extend([(distance[ix][iy],ix,iy) for iy in range(rest,n)])
    dists.sort()
    return dists[0]

# compute the total distance from a color to the other colors in the set

def total_distance(ix):
    n = len(distance)
    d = 0
    d += sum([distance[iy][ix] for iy in range(ix)])
    d += sum([distance[ix][iy] for iy in range(ix+1,n)])
    return d

# find the euclidean distance between two colors, in rgb space (actually, the
# square of the distance)

def color_dist((r1,g1,b1),(r2,g2,b2)):
    return (r1-r2)*(r1-r2) + (g1-g2)*(g1-g2) + (b1-b2)*(b1-b2)

# choose a random color from rgb space
# for 'hard' colors, we avoid the dark octant
# for 'soft' colors, we avoid the four darkest octants

def random_color():
    r = random.random()
    g = random.random()
    b = random.random()

    primary = random.randint(1,3)

    if (softColors):
        if (primary != 1): r = .5 + r/2
        if (primary != 2): g = .5 + g/2
        if (primary != 3): b = .5 + b/2
    else:
        if   (primary == 1): r = .5 + r/2
        elif (primary == 2): g = .5 + g/2
        else:                b = .5 + b/2

    return (int((resol+1)*r),int((resol+1)*g),int((resol+1)*b))

##########
# find a good order in which to place the pieces, trying to minimize the number
# of obscured empty cells
#
# We make a very crude approximation here to what obscures what.  We assume
# that
#    (1) the viewing angle is along the line (x,y,z) = (t,t,k-t)
#    (2) all viewing lines are parallel to that
#    (3) the cell (x,y,z) obscures the cells along four such lines:
#          (x+t,y+t,z-t),(x+1+t,y+t,z-t),(x+t,y+1+t,z-t),(x+t,y+t,z-1-t)
##########

# returns a list of piece names

def piece_order(cells,components):
    # compute the partial cells each piece obscures

    pieceObscures = {}
    for piece in components.keys():
        cubies = components[piece]
        obs = []
        ignore = [piece,"-"]
        for (x,y,z) in cubies:
            (xx,yy,zz) = (x+1,y+1,z-1)
            while (xx < xMax) and (yy < yMax) and (zz >= 0):
                if ((xx,yy,zz) not in obs) and (cells[xx][yy][zz] not in ignore):
                    obs.append((xx,yy,zz))
                (xx,yy,zz) = (xx+1,yy+1,zz-1)

            (xx,yy,zz) = (x+1,y,z)
            while (xx < xMax) and (yy < yMax) and (zz >= 0):
                if ((xx,yy,zz) not in obs) and (cells[xx][yy][zz] not in ignore):
                    obs.append((xx,yy,zz))
                (xx,yy,zz) = (xx+1,yy+1,zz-1)

            (xx,yy,zz) = (x,y+1,z)
            while (xx < xMax) and (yy < yMax) and (zz >= 0):
                if ((xx,yy,zz) not in obs) and (cells[xx][yy][zz] not in ignore):
                    obs.append((xx,yy,zz))
                (xx,yy,zz) = (xx+1,yy+1,zz-1)

            (xx,yy,zz) = (x,y,z-1)
            while (xx < xMax) and (yy < yMax) and (zz >= 0):
                if ((xx,yy,zz) not in obs) and (cells[xx][yy][zz] not in ignore):
                    obs.append((xx,yy,zz))
                (xx,yy,zz) = (xx+1,yy+1,zz-1)

        pieceObscures[piece] = obs

        if ("order" in debug):
            print "// %s = %s" % (piece,cubies)
            print "//    obscures %s" % (obs)
            print "//"

    # mark each cell as empty and unobscured, or not-to-be-used

    notUsed  = 0
    empty    = 1
    obscured = 2
    full     = 3

    cellObscured = zero_array(xMax,yMax,zMax)
    for piece in components.keys():
        for (x,y,z) in pieceObscures[piece]:
            cellObscured[x][y][z] = empty

    # determine the piece placement order

    order = []
    piecesLeft = [piece for piece in pieceNames]

    while (piecesLeft != []):

        if ("order" in debug):
            print
            print_cells_across(sys.stdout,"",cells)
            print
            print_obscured_across(sys.stdout,"",cellObscured)
            print

        # for each unused piece, compute how much obscurity would change if we
        # placed it;  this is the number of unoccupied cells that it would
        # obscure minus the number of obscured cells it would occupy

        minPiece = None
        minOrd   = None
        minObs   = None

        for piece in piecesLeft:
            mustChange = False

            # if this piece is not in the minimum order-group, skip it

            if (piece not in orderHints): orderHints[piece] = 0

            mustChange = (minOrd == None) or (orderHints[piece] < minOrd)
            if (not mustChange) and (orderHints[piece] > minOrd): continue

            # compute the change in obscurity from this piece

            cubies = components[piece]
            obs = sum([1 for (x,y,z) in pieceObscures[piece] if (cellObscured[x][y][z] == empty)]) \
                - sum([1 for (x,y,z) in cubies               if (cellObscured[x][y][z] == obscured)])
            if ("order" in debug):
                print "obs(%s) = %d (order %d)" % (piece,obs,orderHints[piece])
            mustChange = mustChange or (minObs == None) or (obs < minObs)
            if (not mustChange) and (obs > minObs): continue

            # compute the tie breakers

            xx = min([x for (x,y,z) in cubies])    # (position of bounding box)
            yy = min([y for (x,y,z) in cubies])
            zz = min([z for (x,y,z) in cubies])

            if (not mustChange):
                if   (zz >  minZ): continue
                elif (zz == minZ):
                    if   (yy <  minY): continue
                    elif (yy == minY):
                        if   (xx <  minX): continue

            minPiece = piece
            minOrd   = orderHints[piece]
            minObs   = obs
            minX     = xx
            minY     = yy
            minZ     = zz

        assert (minPiece != None), "internal error in piece_order()"

        # 'place' this piece

        if ("order" in debug):
            print "place %s (%d) %s" % (minPiece,minObs,cubies)

        order.append(minPiece)
        piecesLeft.remove(minPiece)

        for (x,y,z) in components[minPiece]:
            cellObscured[x][y][z] = full

        for (x,y,z) in pieceObscures[minPiece]:
            if (cellObscured[x][y][z] == empty): cellObscured[x][y][z] = obscured

    return order

##########
# figure out what cells will hide a piece;  this is the union of the seven
# neighbors of each cell
#    (x-1,y,z),(x,y-1,z),(x,y,z+1),(x-1,y-1,z),(x-1,y,z+1),(x,y-1,z+1),
#    (x-1,y+1,z+1)
#
# if any of them are out of bounds, we return None as the result, because the
# piece can never be completely hidden.
##########

def hiding_cells(cells,cubies):
    cubies = [(x,y,z) for (x,y,z) in cubies]
    allHiders = []

    while (True):
        # find how many of the cubies are hidden

        hidden = {}
        for (x,y,z) in cubies:
            hidden[(x-1,y  ,z  )] = True
            hidden[(x-1,y-1,z  )] = True
            hidden[(x  ,y  ,z+1)] = True
            hidden[(x-1,y-1,z  )] = True
            hidden[(x-1,y  ,z+1)] = True
            hidden[(x  ,y-1,z+1)] = True
            hidden[(x-1,y-1,z+1)] = True

        hiders = [(x,y,z) for (x,y,z) in hidden if ((x,y,z) not in cubies)
                                               and ((x,y,z) not in allHiders)]

        allHiders += hiders

        # if any of these is out of bounds, give up

        if (min([x for (x,y,z) in hiders]) <  0):    return None
        if (min([y for (x,y,z) in hiders]) <  0):    return None
        if (max([z for (x,y,z) in hiders]) >= zMax): return None

        # are any of these cells unoccupied?  if so, remove them from the
        # hider list, and go back and see what hides them

        cubies = [(x,y,z) for (x,y,z) in hiders if (cells[x][y][z] == "-")]
        if (cubies == []): break

        for (x,y,z) in cubies:
            allHiders.remove((x,y,z))

    return allHiders

##########
# create a 3D "array" out of nested lists (this is in lieu of relying on other
# packages like Numeric, numpy, etc.)
##########

def zero_array(xMax,yMax=None,zMax=None):
    if (yMax == None):
        return [0 for x in range(xMax)]
    elif (zMax == None):
        return [[0 for y in range(yMax)]
                   for x in range(xMax)]
    else:
        return [[[0 for z in range(zMax)]
                    for y in range(yMax)]
                    for x in range(xMax)]

##########
# print the ini file
##########

def print_ini(f,povName):
    print >>f,"Input_File_Name=%s" % povName
    print >>f
    print >>f,"quality=9"
    print >>f
    print >>f,"+A"
    print >>f,"Sampling_Method=2"
    print >>f,"Antialias_Depth=2"
    print >>f,"Jitter_Amount=0.0"
    print >>f

    if (useClock):
        print >>f,"Initial_Frame=1"
        print >>f,"Final_Frame=%d" % numPieces
        print >>f,"Initial_Clock=1.01"
        print >>f,"Final_Clock=%d.01" % (numPieces+1)
        print >>f
        print >>f,"Cyclic_Animation=on"
        print >>f,"Pause_when_Done=off"

##########
# print the main part of the pov file, the part irrespective of the pieces
##########

def print_pov_main(f):
    print >>f,"//----------"
    print >>f,"// pov ray-tracing file, produced by solution_to_pov_pieces"
    print >>f,"//----------"
    print >>f
    print >>f,"#include \"colors.inc\""
    print >>f
    print >>f,"background { White }"
    print >>f

    print >>f,"//----------"
    print >>f,"// lighting"
    print >>f,"//----------"
    print >>f
    print >>f,"// a ring of lights to light everything"
    print >>f
    print >>f,"#declare numLights=31;"
    print >>f,"#declare ang=0;"
    print >>f,"#while (ang < 359)"
    print >>f,"\tlight_source { <0,0,9000> rgb .7/numLights rotate 40*x rotate ang*z }"
    print >>f,"\t#declare ang=ang+(360/numLights);"
    print >>f,"#end"
    print >>f
    print >>f,"// main light"
    print >>f
    print >>f,"light_source"
    print >>f,"\t{"
    print >>f,"\t<0,0,1000> .8*White"
    print >>f,"\tarea_light <100, 0, 0>, <100, 0, 5>, 7,7"
    print >>f,"\trotate 55*x rotate -70*z"
    print >>f,"\t}"
    print >>f

    if (xMax <= 5) and (yMax <= 5) and  (zMax <= 5):
        camDist = 12
        (x,y,z)    = ".8","-1","3.7"
    elif (xMax <= 7) and (yMax <= 7) and  (zMax <= 7):
        camDist = 17
        (x,y,z)    = "1.5","-1","5"
    elif (xMax <= 10) and (yMax <= 10) and  (zMax <= 10):
        camDist = 24
        (x,y,z)    = "2","-1","7"
    else: # (need to compute these values)
        camDist = 24
        (x,y,z)    = "2","-1","7"

    if (cameraDist != None):
        camDist = cameraDist

    if (cameraTrans != None):
        (x,y,z) = cameraTrans

    camAspect = ""
    if (imageAspect != None):
        camAspect = "*%s" % imageAspect

    print >>f,"//----------"
    print >>f,"// camera"
    print >>f,"//----------"
    print >>f
    print >>f,"#declare camera_elevation=25;"
    print >>f,"#declare camera_bearing=-120;"
    print >>f,"#declare camera_dist=%s;" % camDist
    print >>f,"#declare camera_ratio=1/2;"
    print >>f
    print >>f,"camera"
    print >>f,"\t{"
    print >>f,"\tperspective"
    print >>f,"\tdirection -x"
    print >>f,"\tright y*camera_ratio%s" % camAspect
    print >>f,"\tup z*camera_ratio"
    print >>f,"\tlocation camera_dist*x"
    print >>f,"\trotate -camera_elevation*y"
    print >>f,"\trotate camera_bearing*z"
    print >>f,"\ttranslate <%s,%s,%s>" % (x,y,z)
    print >>f,"\t}"
    print >>f

    print >>f,"//----------"
    print >>f,"// the floor"
    print >>f,"//----------"
    print >>f
    print >>f,"plane"
    print >>f,"\t{"
    print >>f,"\tz, -0.1"
    print >>f,"\tpigment { checker .95*White, 1*White scale 1 }"
    print >>f,"\tfinish  { ambient .2 diffuse .8 }"
    print >>f,"\t}"
    print >>f

    print >>f,"//----------"
    print >>f,"// cubie definition"
    print >>f,"//----------"
    print >>f
    print >>f,"#declare cubieCenter=<0.5,0.5,0.5>;"
    print >>f
    print >>f,"#declare r=0.05;"
    print >>f,"#declare ptSW=<r,r,r>;"
    print >>f,"#declare ptSE=<1-r,r,r>;"
    print >>f,"#declare ptNE=<1-r,r,1-r>;"
    print >>f,"#declare ptNW=<r,r,1-r>;"
    print >>f
    print >>f,"#declare threeEdges=union    // rounded edges"
    print >>f,"\t{"
    print >>f,"\tcylinder { ptSW,ptSE,r }   // southern edge (as viewed facing XZ plane)"
    print >>f,"\tsphere   { ptSE,r      }   // southeast corner"
    print >>f,"\tcylinder { ptSE,ptNE,r }   // eastern edge"
    print >>f,"\tsphere   { ptNE,r      }   // northeast corner"
    print >>f,"\tcylinder { ptNE,ptNW,r }   // northern edge"
    print >>f,"\t}"
    print >>f
    print >>f,"#declare cubie=union"
    print >>f,"\t{"
    print >>f,"\tbox { <r,r,0>,<1-r,1-r,1> }"
    print >>f,"\tbox { <r,r,0>,<1-r,1-r,1> translate -cubieCenter rotate 90*x translate cubieCenter }"
    print >>f,"\tbox { <r,r,0>,<1-r,1-r,1> translate -cubieCenter rotate 90*y translate cubieCenter }"
    print >>f,"\tobject { threeEdges }"
    print >>f,"\tobject { threeEdges translate -cubieCenter rotate  90*z translate cubieCenter }"
    print >>f,"\tobject { threeEdges translate -cubieCenter rotate 180*z translate cubieCenter }"
    print >>f,"\tobject { threeEdges translate -cubieCenter rotate 270*z translate cubieCenter }"
    print >>f,"\t}"
    print >>f

    print >>f,"//----------"
    print >>f,"// axes (for debug)"
    print >>f,"//----------"
    print >>f
    if (showAxes): print >>f,"#if(1)"
    else:         print >>f,"#if(0)"
    print >>f,"cylinder { <0,0,0>,<7,0,0>,0.1 texture { pigment { color rgb <1,0,0> } finish { ambient 1 } } }"
    print >>f,"cylinder { <0,0,0>,<0,7,0>,0.1 texture { pigment { color rgb <0,1,0> } finish { ambient 1 } } }"
    print >>f,"cylinder { <0,0,0>,<0,0,7>,0.1 texture { pigment { color rgb <0,0,1> } finish { ambient 1 } } }"
    print >>f,"#end"
    print >>f

    print >>f,"//----------"
    print >>f,"// bounding box"
    print >>f,"//----------"
    print >>f
    print >>f,"#declare r = 0.05;"
    print >>f,"#declare col = <1,0,0>;"
    print >>f

    if (showBox): print >>f,"#if(1)"
    else:         print >>f,"#if(0)"
    if (useClock): print >>f,"#if (clock < %d)" % numPieces
    cylinder = "cylinder { %s,r texture { pigment { color rgb col } finish { ambient 1 } } }"
    print >>f,cylinder % (" <0,0,0>,<%d,0,0>"     % (                xMax))
    print >>f,cylinder % (" <%d,0,0>,<%d,%d,0>"   % (xMax,           xMax,yMax     ))
    print >>f,cylinder % (" <%d,%d,0>,<0,%d,0>"   % (xMax,yMax,      xMax          ))
    print >>f,cylinder % (" <0,%d,0>,<0,0,0>"     % (     yMax                     ))
    print >>f,cylinder % (" <0,0,%d>,<%d,0,%d>"   % (          zMax, xMax,     zMax))
    print >>f,cylinder % (" <%d,0,%d>,<%d,%d,%d>" % (xMax,     zMax, xMax,yMax,zMax))
    print >>f,cylinder % (" <%d,%d,%d>,<0,%d,%d>" % (xMax,yMax,zMax,      yMax,zMax))
    print >>f,cylinder % (" <0,%d,%d>,<0,0,%d>"   % (     yMax,zMax,           zMax))
    print >>f,cylinder % (" <0,0,0>,<0,0,%d>"     % (                          zMax))
    print >>f,cylinder % (" <%d,0,0>,<%d,0,%d>"   % (xMax,           xMax,     zMax))
    print >>f,cylinder % (" <%d,%d,0>,<%d,%d,%d>" % (xMax,yMax,      xMax,yMax,zMax))
    print >>f,cylinder % (" <0,%d,0>,<0,%d,%d>"   % (     yMax,           yMax,zMax))
    if (useClock): print >>f,"#end"
    print >>f,"#end"
    print >>f

    print >>f,"//----------"
    print >>f,"// the polycubes"
    print >>f,"//----------"
    print >>f
    print >>f,"#macro make_piece(p,c)"
    print >>f,"\tunion"
    print >>f,"\t\t{"
    print >>f,"\t\t#declare ix = 0;"
    print >>f,"\t\t#while (ix < dimension_size(p,1))"
    print >>f,"\t\t\tobject { cubie translate p[ix] }"
    print >>f,"\t\t\t#declare ix = ix + 1;"
    print >>f,"\t\t\t#end"
    print >>f
    print >>f,"\t\tpigment { rgb c }"
    print >>f,"\t\tfinish  { phong 0.8 diffuse 0.6 ambient 0.2 }"
    print >>f,"\t\t}"
    print >>f,"#end"
    print >>f

##########
# miscellaneous print routines
##########

def print_cells(f,prefix,cells):
    if (prefix == None): prefix = ""

    w = max([max([max([len("%s" % cells[x][y][z]) for z in range(zMax)])
                                                  for y in range(yMax)])
                                                  for x in range(xMax)])

    for z in range(zMax-1,-1,-1):
        for y in range(yMax-1,-1,-1):
            print >>f,prefix+" ".join(["%-*s" % (w,cells[x][y][z]) for x in range(xMax)])
        if (z != 0): print >>f,prefix.rstrip()
    print >>f


def print_cells_across(f,prefix,cells):
    if (prefix == None): prefix = ""

    w = max([max([max([len("%s" % cells[x][y][z]) for z in range(zMax)])
                                                  for y in range(yMax)])
                                                  for x in range(xMax)])

    for y in range(yMax-1,-1,-1):
        rows = []
        for z in range(zMax-1,-1,-1):
            rows.append(" ".join(["%-*s" % (w,cells[x][y][z]) for x in range(xMax)]))
        print >>f,prefix+"   ".join(rows)


def print_obscured_across(f,prefix,cellObscured):
    xMax = len(cellObscured)
    yMax = len(cellObscured[0])
    zMax = len(cellObscured[0][0])

    for y in range(yMax-1,-1,-1):
        rows = []
        for z in range(zMax-1,-1,-1):
            rows.append(" ".join(["%2d" % cellObscured[x][y][z] for x in range(xMax)]))
        print >>f,prefix+"   ".join(rows)


def cubies_string(cubies):
    return ", ".join([cubie_string(x,y,z) for (x,y,z) in cubies])


def cubie_string(x,y,z):
    return "<%s,%s,%s>" % (x,y,z)


if __name__ == "__main__": main()