Thursday, 31 January 2013

The Moment You Didn't Know You Were Waiting For

Following on from the export side of the Soma Cube work, here's the import side, in Python, to be run as a Blender script. This is more of an exercise in finding one's way around Blender's object model than anything else, so there's not much point in a deep commentary on it. Each shape of the cube is loaded into a separate Mesh object, and each solution of the cube becomes a set of keyframes for those objects.

import math
import bpy
import io

file = io.open

# RGB components are in the range [0,1].
#
# NMesh.materials is a list of Material objects. NMesh.faces[i].mat is the
# index into the mesh's material list of the material to use for the i'th
# face.

# When assigning to a face's vertex index list, a final zero is ignored.
# So, when creating a square face, rotate its vertex list if necessary
# to make sure vertex zero can be specified.
def RawVertices(verts):
  if len(verts) == 4 and verts[-1] == 0:
    return [0] + verts[:-1]
  else:
    return verts

def ReadShapes():
  meshes = eval(file('shapes.txt').read())

  materials = {}
  for (_, (_, faces)) in meshes:
    for i, (_, mat_name) in enumerate(faces):
      if mat_name not in materials:
        n = len(materials)
        material = bpy.data.materials.new(mat_name)
        material.diffuse_color = [(n & (1<<j)) and 1.0 or 0.3 for j in range(3)]
        materials[mat_name] = material

  for i, (name, (verts, faces)) in enumerate(meshes):
    poly = bpy.data.meshes.new(name)

    poly.vertices.add(len(verts))
    for j, v in enumerate(verts):
      poly.vertices[j].co = v

    mesh_materials = []
    for (_, mat_name) in faces:
      if mat_name not in mesh_materials:
        mesh_materials.append(mat_name)
    for mat_name in mesh_materials:
      poly.materials.append(bpy.data.materials[mat_name])

    poly.faces.add(len(faces))
    for j, (verts, mat_name) in enumerate(faces):
      f = poly.faces[j]
      f.vertices_raw = RawVertices(verts)
      f.material_index = mesh_materials.index(mat_name)

    poly.update()
    obj = bpy.data.objects.new(name, poly)
    bpy.context.scene.objects.link(obj)

def ReadSolutions():
  solutions = eval(file('pysolutions.txt').read())
  frames_per_solution = 10

  for i, solution in enumerate(solutions + solutions[:1]):
    for name, (euler, location) in solution:
      props = {'rotation_euler': euler, 'location': location }

      obj = bpy.data.objects[name]

      anim_data = obj.animation_data or obj.animation_data_create()
      action = anim_data.action
      if action is None:
        action = bpy.data.actions.new('%s_Action' % name)
        for prop in 'location', 'rotation_euler':
          for index in range(3):
            action.fcurves.new(prop, index)
        anim_data.action = action

      for curve in action.fcurves:
        value = props[curve.data_path][curve.array_index]
        curve.keyframe_points.insert(i * frames_per_solution, value)

  bpy.context.scene.frame_end = i * frames_per_solution

if __name__ == '__main__':
  ReadShapes()
  ReadSolutions()


... And if you've read this far, here's a little reward - a crude Blender render of the resulting animation. Nothing fancy yet, but it represents progress.



Throwing Shapes

The representation of solutions in the Cube module having been reworked, it's now possible to export data about the shapes and their positions in each solution, in a form that Blender can then load via a short Python script.

Here's the Haskell side of things:

> module Mesh where

> import Data.Map (fromList, keys, lookupIndex)
> import Data.Maybe (catMaybes, fromJust, listToMaybe)

> import System.Process (runCommand)

> import Cube

All the logic for moving a shape around within the cube has equivalents for moving a vertex or face around inside the cube. This will be useful for determining which faces of a shape end up on which faces of the cube in a particular solution.

> vApplyR :: Rotation -> (Int, Int, Int) -> (Int, Int, Int)
> vApplyR (Rotation a n) = times n (r a)
>  where
>   r X (x, y, z) = (x, 3-z, y)
>   r Y (x, y, z) = (z, y, 3-x)
>   r Z (x, y, z) = (3-y, x, z)

> vApplyT :: Translation -> (Int, Int, Int) -> (Int, Int, Int)
> vApplyT (Translation (x', y', z')) (x, y, z) = (x+x', y+y', z+z')

> vApplyRecipe (maj, min, shift) = vApplyR maj . vApplyR min . vApplyT shift

In effect, this function scans all the possible face positions, and determines where the cube is occupied by a form on one side along the X axis but not on the other. These will be the places where the form will need an X-face (i.e. one parallel to the YZ plane). (Locations for the Y- and Z-faces can be found by rotating the form appropriately first.)

> faceCoords :: Form -> [(Int, Int, Int)]
> faceCoords form = [(x, y, z) |
>                    (z, layer)   <- zip [0..] (unform form),
>                    (y, row)     <- zip [0..] layer,
>                    (x, isFace)  <- zip [0..] (markFaces row),
>                    isFace]
>  where markFaces row = zipWith (/=) (' ':row) (row ++ [' '])

A few new data types are needed, to record which side of the solved cube a face belongs to. The 'Show' instances for these have been carefully chosen to make parsing in Python easy.

> data Sign = Negative | Positive deriving (Show)
> data CubeFace = CubeFace Sign Axis
> instance Show CubeFace
>  where show (CubeFace sign axis) = [head $ show axis, head $ show sign]
> newtype FacePaint = FacePaint (Maybe CubeFace)
> instance Show FacePaint
>  where show (FacePaint (Just (CubeFace sign axis))) = [
>         '\"', head $ show axis, head $ show sign, '\"']
>        show (FacePaint Nothing) = "\"xx\""

> type Vert = (Int, Int, Int)

> type Face = (Int, Int, Int)

Now everything necessary is present, to be able to determine the locations of a form's faces for any axis. Note that applyR and vApplyR are used to represent a transformation and its inverse, respectively.

> axisFaces :: Axis -> Form -> [Face]
> axisFaces axis = map (vApplyR $ r axis 3) .
>                  faceCoords .
>                  applyR (r axis 1)
>  where
>   r X _ = Rotation X 0
>   r Y n = Rotation Z n
>   r Z n = Rotation Y n

Given a face on an axis, this function determines the locations of its vertices.

> faceVerts axis (x, y, z) = case axis of
>                              X -> [(x, y+y', z+z') | (y', z') <- square]
>                              Y -> [(x+x', y, z+z') | (x', z') <- square]
>                              Z -> [(x-x', y+y', z) | (x', y') <- square]
>  where square = [(0,0), (0,1), (1,1), (1,0)]

Knowing where a face's vertices lie (after transformation) is enough to determine which face of the cube (if any) it lies on.

> paintFace verts = listToMaybe . catMaybes $
>                   [paint vs axis | (vs, axis) <- [(xs, X), (ys, Y), (zs, Z)]]
>  where
>   (xs, ys, zs) = unzip3 verts
>   paint vs axis = if all (== 0) vs then Just $ CubeFace Negative axis
>                   else if all (==3) vs then Just $ CubeFace Positive axis
>                   else Nothing

Finally, all these things are brought together to give a way to describe how to position and colour all of a shape's faces to suit the space it occupies in a solution.

> axisInfo :: Recipe -> Axis -> Form -> [([Vert], FacePaint)]
> axisInfo recipe axis form = [(faceVerts axis face,
>                               FacePaint $ paintFace $ map (vApplyRecipe recipe) $ faceVerts axis face) |
>                              face <- axisFaces axis form]

> allFaces recipe f = concat [axisInfo recipe axis f | axis <- enumFrom X]

This little helper takes a list of faces (each having a list of vertices) and pulls all the vertices out into a separate (deduplicated) list, replacing them in the original data structure by indices into that vertex list.

> shareVerts :: Ord k => [([k], t)] -> ([k], [([Int], t)])
> shareVerts faces = (keys vertMap,
>                     [([fromJust (lookupIndex v vertMap) | v <- verts], x) |
>                      (verts, x) <- faces])
>  where vertMap = fromList [(v, ()) | v <- concat (map fst faces)]

Because the colouring of a shape depends on the choice of solution, the mesh can't be generated without specifying a recipe for the shape.

> mesh recipe = shareVerts . allFaces recipe . defaultForm
> allMeshes (Solution fs) = [(show shape, mesh recipe shape) | (shape, recipe, _) <- fs]

The file emitted here is one big list of mesh data. Luckily, 'show' for Haskell's built-in types and 'eval' in Python are near enough inverses that, with some careful implementations of 'Show' (as noted above) it's possible to avoid writing a parser for it at all.

> writeMeshes = writeFile "shapes.txt" $ show $ allMeshes (head solutions)

This rather ugly function describes the two rotations and a translation that compose a Recipe, in a form that can easily be parsed in Python.

> pyRecipe recipe@(maj, min, shift) = ((roll, pitch, heading), location)
>  where
>   location = vApplyRecipe recipe (0, 0, 0)
>   heading =
>     case maj of Rotation X _ -> error "X should only appear as a minor axis"
>                 Rotation Y _ -> 0
>                 Rotation Z n -> quarterTurns n
>   pitch =
>     case maj of Rotation X _ -> error "X should only appear as a minor axis"
>                 Rotation Y n -> quarterTurns n
>                 Rotation Z _ -> 0
>   roll =
>     case min of Rotation X n -> quarterTurns n
>                 _            -> error "Only X should appear as a minor axis"
>   quarterTurns n = (fromIntegral n) * pi/2

> pySolution (Solution fs) = [(show shape, pyRecipe recipe) |
>                             (shape, recipe, _) <- fs]

> writeSolutions = writeFile "pysolutions.txt" $
>                  show [pySolution solution | solution <- solutions]

Lastly, a test command glues everything together: writing out the data, then loading it into Blender.

> testMeshes = do
>   writeMeshes
>   writeSolutions
>   _ <- runCommand $ unwords ["blender", "-P", "blender.py"]
>   return ()

Tedious Groundwork

I've been at the cube again. Quite a lot has developed, but it's all been built upon some relatively small but important changes to the way that the original cube solving code models its work. Specifically, each solution now carries with it not only information about which shapes occupy which cells, but also information about how to translate and rotate each piece from its default position into its position within the solution.

This is crucial for tasks such as exporting solution geometry to a 3D modelling package such as Blender ... of which more later. For the moment, here's a very boring rehash of the existing code, for completeness.

> module Cube where

> import Data.Bits
> import Data.Function (on)
> import Data.List (intercalate, nubBy, sortBy, transpose, unfoldr)
> import Data.Monoid
> import Data.Ord (comparing)
> import Data.Word

The shapes remain the same. The last one has been slightly renamed to avoid conflict with the axis of the same name.

> data Shape = L | S | T | R | P | Q | Y'
>   deriving (Enum, Show)

> plan L = "XXX" ++
>          "X  "

> plan S = " XX" ++
>          "XX "

> plan T = "XXX" ++
>          " X "

> plan R = "XX " ++
>          "X  "

> plan P = "X  " ++
>          "X  " ++
>          "   " ++
>          "XX "

> plan Q = " X " ++
>          " X " ++
>          "   " ++
>          "XX "

> plan Y' = "XX " ++
>           "X  " ++
>           "   " ++
>           "X  "

> newtype Form = Form [[[Char]]]
> unform (Form x) = x

> instance Eq Form where
>   (Form x) == (Form y) = x == y

> instance Ord Form where
>   Form x `compare` Form y = x `compare` y

> withForm f = Form . f . unform

The order in which rows of a form are displayed has been reversed, so that Y increases up the screen. With X increasing to the right, this means that each layer of the cube is displayed as if viewed from above (looking towards -Z) in a right-handed coordinate system, so it fits Blender naturally.

> instance Show Form where
>   show (Form f) = unlines [showSlice s | s <- reverse $ transpose f]
>    where
>     showSlice s = intercalate "  " [ showRow r | r <- s ]
>     showRow r = "|" ++ r ++ "|"

> instance Monoid Form where
>   mempty = Form $ (replicate 3 . replicate 3 . replicate 3) ' '
>   mappend (Form x) (Form y) = Form $ (zipWith . zipWith . zipWith)
>                                      superpose x y
>    where
>     superpose ' ' ' ' = ' '
>     superpose  a  ' ' =  a
>     superpose ' '  b  =  b
>     superpose  _   _  = '#'  -- This should never happen!

The following splits a list into chunks, each of which (except maybe the last) has length 3. This would naturally generalise to chunks of length N if needed. This was previously defined within 'defaultForm' but is now more widely useful.

> inThrees = takeWhile (not . null) . unfoldr (Just . splitAt 3)

> flatten = concat . concat . unform

As well as flattening a form into a list, the inverse operation will also be necessary. This takes a flat list and produces a list of triples of triples.

> unflatten = inThrees . inThrees

> asBits :: Form -> Word32
> asBits f = sum [bit n | (n, x) <- zip [0..] (flatten f), x /= ' ']

The opposite of 'asBits', this takes a shape and a bit sequence representing its form, and reconstructs the form. The ability to do this means that computations can take place purely in terms of bit sequences, and the results can still be displayed in a comprehensible manner by expanding them afterwards.

> fromBits shape bs = Form $ unflatten [formChar (testBit bs n) | n <- [0..26]]
>  where
>   formChar False = ' '
>   formChar True  = head (show shape)

> size f = length [x | x <- flatten f, x /= ' ']

> defaultForm :: Shape -> Form
> defaultForm shape = Form $ unflatten fullPlan
>  where
>   fullPlan = [formChar c | c <- plan shape] ++
>              replicate (27 - length (plan shape)) ' '
>   formChar ' ' = ' '
>   formChar  _  = head (show shape)

To be able to produce geometrical descriptions of the solutions, it is useful to record not only which shapes occupy which cell of the cube, but also how the shapes' forms are composed. As mentioned before, each form can be constructed from a translation, followed by a rotation about the X axis and then a rotation about either the Y or Z axes.

These transformations could all be represented as matrices, but they're much easier to decipher if given a higher-level representation. This sequence of three transformations is referred to as a recipe.

> data Axis = X | Y | Z deriving (Enum, Show)
> 
> newtype Translation = Translation (Int, Int, Int) deriving Show
> 
> type QuarterTurns = Int
> data Rotation = Rotation Axis QuarterTurns deriving Show
> 
> type Recipe = (Rotation, Rotation, Translation)

Having given names to these transformations, they can now be listed in terms of those names instead of relying directly on function composition.

> majors = [Rotation Y n | n <- [0..3]] ++
>          [Rotation Z n | n <- [1, 3]]
> 
> minors = [Rotation X n | n <- [0..3]]
> 
> shifts shape = [t | t <- [Translation (x, y, z) |
>                           x <- [0..2],
>                           y <- [0..2],
>                           z <- [0..2]], inBounds t]
>  where
>   inBounds t = size (applyT t $ df) == size df
>   df = defaultForm shape

Translations can now be in either direction; this might prove useful later.

> shiftr def [a, b, _] = [def, a, b]
> shiftl def [_, b, c] = [b, c, def]

> withx op = (map . map) (op ' ')
> withy op = map (op "   ")
> withz op = op (replicate 3 "   ")
> tx  = withx shiftr
> tx' = withx shiftl
> ty  = withy shiftr
> ty' = withy shiftl
> tz  = withz shiftr
> tz' = withz shiftl

> cw = transpose . reverse
> ccw = reverse . transpose

Transformations being defined in terms of Recipes instead of raw functions, a little glue is needed to describe how a recipe is applied to a Form.

Note also that the rotation about X has been changed from clockwise to counter-clockwise. These rotations now match Blender's default settings, which will avoid complication further down the line.

> times 0 f = id
> times n f = f . times (n-1) f
> 
> applyR (Rotation a n) = withForm (times n (r a))
>  where
>   r X = cw
>   r Y = transpose . map ccw . transpose
>   r Z = map cw
> 
> applyT (Translation (x, y, z)) = withForm (times x tx . times y ty . times z tz)
> 
> applyRecipe (maj, min, shift) = applyR maj . applyR min . applyT shift

The special case for the forms of L have been described a little more neatly.

> allForms :: Shape -> [(Recipe, Form)]
> allForms shape = nubBy ((==) `on` snd) $ [(recipe, applyRecipe recipe df) |
>                                           recipe <- recipes shape]
>  where
>   df = defaultForm shape
>   recipes L = [(Rotation Y 0, Rotation X 0, t) |
>                t <- [s | s@(Translation (a, b, c)) <- shifts L, c /= 2]]
>   recipes shape =  [(maj, min, shift) |
>                     maj <- majors,
>                     min <- minors,
>                     shift <- shifts shape]

The list of possible forms for each shape can now carry the recipe for constructing each form, as well as the resulting bit pattern.

> combos :: [(Shape, [(Recipe, Word32)])]
> combos = sortBy (comparing (length . snd))
>                 [(shape, [(r, asBits f) | (r, f) <- allForms shape]) |
>                  shape <- enumFrom L]

> newtype Solution = Solution [(Shape, Recipe, Word32)] deriving (Show)

> solutions :: [Solution]
> solutions = [s | (s, _) <- foldl addForm [(Solution [], 0)] combos] where
>   addForm partialSolutions (shape, formData) = [
>       (Solution ((shape, recipe, bitmap):s), (pb .|. bitmap)) |
>       (Solution s, pb) <- partialSolutions,
>       (recipe, bitmap) <- formData,
>       pb .&. bitmap == 0
>       ]

> expandSolution (Solution s) = foldl1 mappend [fromBits shape bitmap |
>                                               (shape, _, bitmap) <- s]

Finally, this makes it possible to explain what sequence of rotations and transformations was used to place each shape in position in a solution.

> explainSolution (Solution s) = unlines [show (shape, recipe) |
>                                         (shape, recipe, _) <- s]

> main = putStrLn $ intercalate "\n" [show (expandSolution s) | s <- solutions]