#! /usr/bin/python
### -*- coding: utf-8 -*-
###
### Rubik cube model and rendering
###
### (c) 2018 Mark Wooding
###

###----- Licensing notice ---------------------------------------------------
###
### This program is free software; you can redistribute it and/or modify
### it under the terms of the GNU General Public License as published by
### the Free Software Foundation; either version 2 of the License, or
### (at your option) any later version.
###
### This program is distributed in the hope that it will be useful,
### but WITHOUT ANY WARRANTY; without even the implied warranty of
### MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
### GNU General Public License for more details.
###
### You should have received a copy of the GNU General Public License
### along with this program; if not, write to the Free Software Foundation,
### Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.

###--------------------------------------------------------------------------
### External dependencies.

import math as M; TAU = 2*M.pi
import os as OS
import sys as SYS

import cairo as CR
import numpy as NP

###--------------------------------------------------------------------------
### The main cube model.

## Names for the faces and corresponding colours.
U, D, L, R, F, B, UNK = 0, 1, 2, 3, 4, 5, 6
FACE = ['U', 'D', 'L', 'R', 'F', 'B']

## Number of faces.
NFACE = 6

## Colours of the faces.
COLOUR = ['#', '+', '%', '/', '=', '-', '·']
RGB = [(1, 0, 0), (1, 0.75, 0),
       (1, 1, 0), (1, 1, 1),
       (0, 0, 1), (0, 1, 0),
       (0.4, 0.4, 0.4)]
CMAP = { 'r': U, 'o': D, 'y': L, 'w': R, 'b': F, 'g': B, '?': UNK }

## Geometry table.
##
## For 0 <= k < 6 and 0 <= d < 4, GEOM[k][d] holds a triple (k', ix, iy) with
## the following meaning.  If you start at face k, move in the direction d
## (up, down, left, or right), and treat the edge that you crossed as the x
## axis of a right-handed coordinate system, increasing clockwise, then you
## are on face k', and (ix, iy) is a unit vector in the coordinate system of
## the state array pointing in the positive x direction.
GEOM = [[(B, -1,  0), (F, -1,  0), (L, -1,  0), (R, -1,  0)],
        [(F, +1,  0), (B, +1,  0), (L, +1,  0), (R, +1,  0)],
        [(U,  0, -1), (D,  0, -1), (B,  0, +1), (F,  0, -1)],
        [(U,  0, +1), (D,  0, +1), (F,  0, +1), (B,  0, -1)],
        [(U, +1,  0), (D, -1,  0), (L,  0, +1), (R,  0, -1)],
        [(U, -1,  0), (D, +1,  0), (R,  0, +1), (L,  0, -1)]]

def convert_side_coords(w, orig, face, x, y):
  """
  Convert relative coordinates (X, Y) on a side face into raw coordinates.

  Suppose we start on face ORIG and move to the adjacent FACE (one of `U',
  `D', `L', `R'), treating the edge we crossed as the x-axis of a right-hand
  coordinate system, increasing clockwise; return a triple (K, U, V), such
  that m[K][U][V] identifies the slot in the model corresponding to the given
  coordinates.  The argument W is the side length of the cube.
  """
  k, ix, iy = GEOM[orig][face]
  jx, jy = -iy, ix
  if ix < 0 or jx < 0: ox = w - 1
  else: ox = 0
  if iy < 0 or jy < 0: oy = w - 1
  else: oy = 0
  return k, ox + x*ix + y*jx, oy + x*iy + y*jy

def slice_orbits(w, face, dp):
  """
  Return the effect of a slice turn on the side stickers.

  Returns a list of lists of (K, X, Y) model coordinates describing the cycle
  decomposition of the permutation induced on the side stickers of a
  clockwise turn of the DPth slice of a size-W cube as viewed from FACE.
  """
  oo = []
  for x in xrange(w):
    oo.append([convert_side_coords(w, face, k, x, dp)
               for k in [U, R, D, L]])
  return oo

def face_orbits(w, face, inv):
  """
  Return the effect of a face turn on the face's stickers.

  Returns a list of lists of (K, X, Y) model coordinates describing the cycle
  decomposition of the permutation induced on the face stickers of a
  clockwise (or anticlockwise, if INV is true) turn of a FACE of a size-W
  cube.
  """
  if inv: ta, tb, tc, td, dx, dy = 0, -1, 1, 0, w - 1, 0
  else: ta, tb, tc, td, dx, dy = 0, 1, -1, 0, 0, w - 1
  oo = []
  for x in xrange((w + 1)//2, w):
    for y in xrange(w//2, x + 1):
      o = []
      u, v = x, y
      for i in xrange(4):
        o.append((face, u, v))
        u, v = ta*u + tb*v + dx, tc*u + td*v + dy
      oo.append(o)
  return oo

def orbits(w, face, dp):
  """
  Return the effect of a slice turn on the cube's stickers.

  Returns a list of lists of (K, X, Y) model coordinates describing the cycle
  decompositon of the permutation induced on the cube's stickers of a
  clockwise turn of the DPth slice of a size-W cube, as viewed from FACE.
  """
  oo = slice_orbits(w, face, dp)
  if dp == 0: oo += face_orbits(w, face, False)
  elif dp == w - 1: oo += face_orbits(w, face ^ 1, True)
  return oo

def maybe_parens(parens, s):
  """If S, then return the string `(S)'; otherwise just return S."""
  if parens: return '(' + s + ')'
  else: return s

class BaseMove (object):
  """
  Base class for move objects.

  Subclasses must implement the following methods.

    * apply(PUZ): apply the move to a puzzle PUZ
    * invert(): return the inverse of the move
    * _str(lv): return a string representation, in a priority-LV context
  """
  def repeat(me, m):
    """Return the move repeated M times."""
    if m == 0: return NULLMOVE
    else: return RepeatedMove(me, m)
  def flatten(me):
    """Return a list of constituent moves."""
    return [me]
  def __str__(me):
    """Return a string representation of the move."""
    return me._str(9)

class NullMove (BaseMove):
  """The null move, which does nothing at all."""
  def apply(me, puz):
    """Apply the move to the puzzle PUZ."""
    pass
  def invert(me):
    """Return the inverse of the move."""
    return me
  def repeat(me, m):
    """Return the move repeated M times."""
    return me
  def flatten(me):
    """Return a list of constituent moves."""
    return []
  def _str(me, lv):
    """Return a string representation, in a priority-LV context."""
    return '#<null>'
NULLMOVE = NullMove()

class AtomicMove (BaseMove):
  """A turn of a single slice."""

  MAP = { }
  for k in xrange(NFACE):
    MAP[FACE[k]] = (k, 0)
    MAP[FACE[k].lower()] = (k, 1)

  def __init__(me, face, dp = 0, n = 1):
    """
    Construct an atomic move.

    Returns a move representing N clockwise turns of the DPth slice, as
    viewed from FACE.
    """
    me.face = face
    me.dp = dp
    me.n = n
  def apply(me, puz):
    """Apply the move to the puzzle PUZ."""
    if me.dp >= 0: puz.turn(me.face, me.dp, me.n)
    else: puz.turn_multi(me.face, -me.dp, me.n)
  def invert(me):
    """Return the inverse of the move."""
    return AtomicMove(me.face, me.dp, -me.n)
  def repeat(me, m):
    """Return the move repeated M times."""
    mm = (me.n*m)%4
    if not mm: return NULLMOVE
    else: return AtomicMove(me.face, me.dp, mm)
  def _str(me, lv):
    """Return a string representation, in a priority-LV context."""
    f = FACE[me.face]
    if me.n%4 == 0: d = '0'
    elif me.n%4 == 1: d = ''
    elif me.n%4 == 2: d = '2'
    elif me.n%4 == 3: d = "'"
    if me.dp == 0: return f + d
    elif me.dp == 1: return f.lower() + d
    elif me.dp < 0: return '%s%s_%d' % (f, d, -me.dp)

class SpinMove (BaseMove):
  """A spin of the entire cube."""
  MAP = { 'X': R, 'Y': U, 'Z': F }
  UNMAP = { }
  for k, v in MAP.iteritems(): UNMAP[v] = k

  def __init__(me, face, n = 1):
    """Construct a spin move.

    Returns a move representing N clockwise spins the entire cube DPth slice,
    as viewed from FACE.
    """
    me.face = face
    me.n = n
  def apply(me, puz):
    """Apply the move to the puzzle PUZ."""
    puz.turn_multi(me.face, puz.w, me.n)
  def invert(me):
    """Return the inverse of the move."""
    return SpinMove(me.face, -me.n)
  def repeat(me, m):
    """Return the move repeated M times."""
    mm = (me.n*m)%4
    if not mm: return NULLMOVE
    else: return SpinMove(me.face, mm)
  def _str(me, lv):
    """Return a string representation, in a priority-LV context."""
    f = me.MAP[me.face]
    if me.n%4 == 0: d = '0'
    elif me.n%4 == 1: d = ''
    elif me.n%4 == 2: d = '2'
    elif me.n%4 == 3: d = "'"
    return f + d

class MoveSequence (BaseMove):
  """A sequence of moves, applied in left-to-right order."""
  def __init__(me, moves):
    """Construct a sequence of moves from the list MOVES."""
    me.moves = reduce(lambda x, y: x + y, [mv.flatten() for mv in moves])
  def apply(me, puz):
    """Apply the move to the puzzle PUZ."""
    for mv in me.moves: mv.apply(puz)
  def invert(me):
    """Return the inverse of the move."""
    return MoveSequence(list(reversed([mv.invert() for mv in me.moves])))
  def _str(me, lv):
    """Return a string representation, in a priority-LV context."""
    return maybe_parens(lv < 5, ' '.join([mv._str(5) for mv in me.moves]))

class RepeatedMove (BaseMove):
  """A move applied multiple times."""
  def __init__(me, move, m):
    """Construct the move which consists of MOVE repeated M times."""
    me.move = move
    me.m = m
  def apply(me, puz):
    """Apply the move to the puzzle PUZ."""
    for i in xrange(me.m): me.move.apply(puz)
  def invert(me):
    """Return the inverse of the move."""
    return RepeatedMove(me.move.invert(), me.m)
  def repeat(me, m):
    """Return the move repeated M times."""
    if m == 0: return NULLMOVE
    else: return RepeatedMove(me.move, m*me.m)
  def _str(me, lv):
    """Return a string representation, in a priority-LV context."""
    return maybe_parens(lv < 3, '%s%d' % (me.move._str(2), me.m))

def skip_spaces(s, i):
  """Return the smallest index J >= I such that S[J] is not whitespace."""
  while i < len(s) and s[i].isspace(): i += 1
  return i

def parse_number(s, i):
  """
  Parse an integer from S, starting at index I.

  The substring S[I:] consists of zero or more digits, optionally followed by
  a nondigit character and any other characters.  Convert the digit sequence
  into an integer N and return the pair (J, N), where J is the index of the
  first non-whitespace character following the digit sequence, or the length
  of S if there is no such character.
  """
  n = 0
  while i < len(s) and s[i].isdigit():
    n = 10*n + ord(s[i]) - ord('0')
    i += 1
  i = skip_spaces(s, i)
  return i, n

def parse_parens(s, i, delim):
  i = skip_spaces(s, i)
  i, mv = parse_conjugate(s, i)
  if s[i] != delim: raise SyntaxError('missing %s' % delim)
  i = skip_spaces(s, i + 1)
  return i, mv

def parse_primary(s, i):
  if i >= len(s): return i, None
  elif s[i] == '(': return parse_parens(s, i + 1, ')')
  elif s[i] == '[': return parse_parens(s, i + 1, ']')
  elif s[i] in SpinMove.MAP:
    k = SpinMove.MAP[s[i]]
    i = skip_spaces(s, i + 1)
    return i, SpinMove(k, 1)
  elif s[i] in AtomicMove.MAP:
    k, dp = AtomicMove.MAP[s[i]]
    if dp == 0 and i + 2 < len(s) and s[i + 1] == '_' and s[i + 2].isdigit():
      i, n = parse_number(s, i + 2)
      if n == 1: dp = 0
      else: dp = -n
    else:
      i = skip_spaces(s, i + 1)
    mv = AtomicMove(k, dp, 1)
    return i, mv
  else: return i, None

def parse_move(s, i):
  i, mv = parse_primary(s, i)
  if mv is None: return i, None
  while True:
    if i < len(s) and s[i] == "'":
      i = skip_spaces(s, i + 1)
      mv = mv.invert()
    elif i < len(s) and s[i].isdigit():
      i, m = parse_number(s, i)
      mv = mv.repeat(m)
    elif i + 1 < len(s) and s[i] == '^' and s[i + 1].isdigit():
      i, m = parse_number(s, i + 1)
      mv = mv.repeat(m)
    else:
      break
  return i, mv

def parse_sequence(s, i):
  seq = []
  while True:
    i, mv = parse_move(s, i)
    if mv is None: break
    seq.append(mv)
  if len(seq) == 0: raise SyntaxError('unrecognized move')
  elif len(seq) == 1: return i, seq[0]
  else: return i, MoveSequence(seq)

def parse_commutator(s, i):
  i, mv = parse_sequence(s, i)
  if i < len(s) and s[i] == ';':
    i = skip_spaces(s, i + 1)
    i, mv1 = parse_sequence(s, i)
    mv = MoveSequence([mv, mv1, mv.invert(), mv1.invert()])
  return i, mv

def parse_conjugate(s, i):
  i, mv = parse_commutator(s, i)
  if i < len(s) and (s[i] == '.' or s[i] == '*'):
    i = skip_spaces(s, i + 1)
    i, mv1 = parse_commutator(s, i)
    mv = MoveSequence([mv, mv1, mv.invert()])
  return i, mv

def parse_moves(s, i = 0):
  i = skip_spaces(s, i)
  i, mv = parse_conjugate(s, i)
  i = skip_spaces(s, i)
  if i < len(s): raise SyntaxError('junk at eol')
  return mv

def xltmatrix(dx, dy, dz):
  return NP.matrix([[1, 0, 0, dx],
                    [0, 1, 0, dy],
                    [0, 0, 1, dz],
                    [0, 0, 0, 1]], NP.float64)

def rotmatrix(theta, axis):
  R = NP.matrix(NP.identity(4, NP.float64))
  i, j = 0, 1
  if i >= axis: i += 1
  if j >= axis: j += 1
  c, s = M.cos(theta), M.sin(theta)
  R[i, i] = R[j, j] = c
  R[i, j] = s; R[j, i] = -s
  return R

def scalematrix(sx, sy, sz):
  return NP.matrix([[sx, 0, 0, 0],
                    [0, sy, 0, 0],
                    [0, 0, sz, 0],
                    [0, 0, 0, 1]], NP.float64)

class Transform (object):
  def __init__(me, theta_x, theta_y, delta_z):
    me.theta_x = theta_x
    me.theta_y = theta_y
    me.delta_z = delta_z
    me._T = xltmatrix(0, 0, delta_z)*\
            rotmatrix(theta_x, 0)*rotmatrix(theta_y, 1)
  def project(me, x, y, z):
    u, v, w, _ = (me._T*NP.array([x, y, z, 1])[:, NP.newaxis]).A1
    return me.delta_z*u/w, me.delta_z*v/w

def vector(x, y, z):
  return NP.array([x, y, z], NP.float64)

def unit_vector(v):
  mag = M.sqrt(NP.dot(v, v))
  return v/mag

def shorten(p, q, delta):
  u = unit_vector(q - p)
  return p + delta*u, q - delta*u

def along(f, p, q):
  return p + f*(q - p)

class BasePainter (object):
  def __init__(me, xf, cr):
    me._cr = cr
    me._xf = xf
    dx, _ = me._cr.device_to_user_distance(0.4, 0)
    me._cr.set_line_width(dx)
  def set_colour(me, r, g, b):
    me._cr.set_source_rgb(r, g, b)
    return me
  def path(me, pts):
    firstp = True
    for x, y, z in pts:
      u, v = me._xf.project(x, y, z)
      if firstp: me._cr.move_to(u, v)
      else: me._cr.line_to(u, v)
      firstp = False
    return me
  def _rounded_path_corner(me, r, q1, p1, p2):
    q2, q3 = shorten(p1, p2, r)
    c1, c2 = along(0.7, q1, p1), along(0.7, q2, p1)
    x0, y0 = me._xf.project(*c1)
    x1, y1 = me._xf.project(*c2)
    x2, y2 = me._xf.project(*q2)
    me._cr.curve_to(x0, y0, x1, y1, x2, y2)
    return p2, q3
  def _rounded_path_step(me, r, q1, p1, p2):
    p2, q3 = me._rounded_path_corner(r, q1, p1, p2)
    x3, y3 = me._xf.project(*q3)
    me._cr.line_to(x3, y3)
    return p2, q3
  def rounded_polygon(me, r, pts):
    o0 = o1 = None
    for x, y, z in pts:
      p2 = vector(x, y, z)
      if o0 is None: o0 = p2
      elif o1 is None:
        o1 = p1 = p2
        q0, q1 = shorten(o0, p1, r)
        x0, y0 = me._xf.project(*q1)
        me._cr.move_to(x0, y0)
      else: p1, q1 = me._rounded_path_step(r, q1, p1, p2)
    p1, q1 = me._rounded_path_step(r, q1, p1, o0)
    p1, q1 = me._rounded_path_corner(r, q1, p1, o1)
    me._cr.close_path()
    return me
  def polygon(me, pts):
    me.path(pts)
    me._cr.close_path()
    return me

class Painter (BasePainter):
  def stroke(me):
    me._cr.stroke()
    return me
  def fill(me):
    me._cr.fill()
    return me

class MeasuringPainter (BasePainter):
  def __init__(me, *args, **kw):
    super(MeasuringPainter, me).__init__(*args, **kw)
    me.x0 = me.y0 = me.x1 = me.y1 = None
  def _update_bbox(me, x0, y0, x1, y1):
    if me.x0 is None or me.x0 > x0: me.x0 = x0
    if me.y0 is None or me.y0 > y0: me.y0 = y0
    if me.x1 is None or me.x1 < x1: me.x1 = x1
    if me.y1 is None or me.y1 < y1: me.y1 = y1
  def stroke(me):
    x0, y0, x1, y1 = me._cr.stroke_extents()
    me._cr.new_path()
    me._update_bbox(x0, y0, x1, y1)
    return me
  def fill(me):
    x0, y0, x1, y1 = me._cr.fill_extents()
    me._cr.new_path()
    me._update_bbox(x0, y0, x1, y1)
    return me

class Puzzle (object):

  def __init__(me, w):
    me.w = w
    me.m = NP.ndarray((NFACE, w, w), NP.uint8)

    for k in xrange(NFACE):
      for y in reversed(xrange(w)):
        for x in xrange(w):
          me.m[k, x, y] = k

    me.orbits = [[orbits(w, k, dp) for dp in xrange(w)]
                 for k in xrange(NFACE)]

  def turn(me, face, dp, n):
    for o in me.orbits[face][dp]:
      c = [me.m[k, x, y] for k, x, y in o]
      for i in xrange(4):
        k, x, y = o[(i + n)%4]
        me.m[k, x, y] = c[i]
    return me

  def turn_multi(me, face, dp, n):
    for i in xrange(dp): me.turn(face, i, n)

  def fill_face(me, s, i, k, ox, oy, ix, iy, jx, jy):
    i = skip_spaces(s, i)
    if i >= len(s): return i
    firstp = True
    for y in xrange(me.w):
      if firstp: pass
      elif i >= len(s) or s[i] != ',': raise SyntaxError('missing ,')
      else: i += 1
      for x in xrange(me.w):
        if i >= len(s): raise SyntaxError('unexpected eof')
        c = CMAP[s[i]]
        me.m[k, ox + x*ix + y*jx, oy + x*iy + y*jy] = c
        i += 1
      firstp = False
    if i < len(s) and not s[i].isspace():
      raise SyntaxError('face spec too long')
    return i

  def fill(me, s):
    i = 0
    me.m[:, :, :] = UNK
    i = me.fill_face(s, i, U, 0, 2, +1,  0,  0, -1)
    i = me.fill_face(s, i, F, 0, 2, +1,  0,  0, -1)
    i = me.fill_face(s, i, R, 0, 2, +1,  0,  0, -1)
    i = me.fill_face(s, i, B, 2, 2, -1,  0,  0, -1)
    i = me.fill_face(s, i, L, 2, 2, -1,  0,  0, -1)
    i = me.fill_face(s, i, D, 0, 0, +1,  0,  0, +1)
    i = skip_spaces(s, i)
    if i < len(s): raise SyntaxError('trailing junk')
    return me

  def show(me):
    for y in reversed(xrange(me.w)):
      print (2*me.w + 2)*' ' + \
        ' '.join([COLOUR[me.m[U, x, y]]
                  for x in xrange(me.w)])
    print
    for y in reversed(xrange(me.w)):
      print '   '.join([' '.join([COLOUR[me.m[k, x, y]]
                                  for x in xrange(me.w)])
                        for k in [L, F, R, B]])
    print
    for y in reversed(xrange(me.w)):
      print (2*me.w + 2)*' ' + \
        ' '.join([COLOUR[me.m[D, x, y]]
                  for x in xrange(me.w)])
    print

  def render_face(me, paint, k, ox, oy, oz, ix, iy, iz, jx, jy, jz):
    w = 2.0/me.w
    g = w/16.0
    for x in xrange(me.w):
      for y in xrange(me.w):
        paint.set_colour(*RGB[me.m[k, x, y]])
        paint.rounded_polygon(0.04,
                              [(ox + (x*w + g)*ix + (y*w + g)*jx,
                                oy + (x*w + g)*iy + (y*w + g)*jy,
                                oz + (x*w + g)*iz + (y*w + g)*jz),
                               (ox + ((x + 1)*w - g)*ix + (y*w + g)*jx,
                                oy + ((x + 1)*w - g)*iy + (y*w + g)*jy,
                                oz + ((x + 1)*w - g)*iz + (y*w + g)*jz),
                               (ox + ((x + 1)*w - g)*ix + ((y + 1)*w - g)*jx,
                                oy + ((x + 1)*w - g)*iy + ((y + 1)*w - g)*jy,
                                oz + ((x + 1)*w - g)*iz + ((y + 1)*w - g)*jz),
                               (ox + (x*w + g)*ix + ((y + 1)*w - g)*jx,
                                oy + (x*w + g)*iy + ((y + 1)*w - g)*jy,
                                oz + (x*w + g)*iz + ((y + 1)*w - g)*jz)]).fill()

  def render(me, paint, reflectp):

    if reflectp:
      paint.set_colour(0, 0, 0)
      paint.rounded_polygon(0.04, [(-1, +1, +4), (+1, +1, +4),
                                   (+1, -1, +4), (-1, -1, +4)]).fill()
      paint.rounded_polygon(0.04, [(-1, -4, +1), (+1, -4, +1),
                                   (+1, -4, -1), (-1, -4, -1)]).fill()
      paint.rounded_polygon(0.04, [(-4, -1, +1), (-4, +1, +1),
                                   (-4, +1, -1), (-4, -1, -1)]).fill()
      me.render_face(paint, D, -1, -4, +1, +1,  0,  0,  0,  0, -1)
      me.render_face(paint, B, +1, -1, +4, -1,  0,  0,  0, +1,  0)
      me.render_face(paint, L, -4, -1, +1,  0,  0, -1,  0, +1,  0)

    paint.set_colour(0, 0, 0)
    paint.rounded_polygon(0.04, [(-1, +1, -1), (-1, +1, +1),
                                 (+1, +1, +1), (+1, -1, +1),
                                 (+1, -1, -1), (-1, -1, -1)]).fill()

    ##paint.set_colour(0.3, 0.3, 0.3)
    ##paint.path([(-1, +1, -1), (+1, +1, -1), (+1, +1, +1)]).stroke()
    ##paint.path([(+1, +1, -1), (+1, -1, -1)]).stroke()

    me.render_face(paint, U, -1, +1, -1, +1,  0,  0,  0,  0, +1)
    me.render_face(paint, F, -1, -1, -1, +1,  0,  0,  0, +1,  0)
    me.render_face(paint, R, +1, -1, -1,  0,  0, +1,  0, +1,  0)

st = Puzzle(3)
invertp = False
reflectp = True
scale = 32
for a in SYS.argv[1:]:
  if not a: pass
  elif a[0] == '-' or a[0] == '+':
    w = a[1:]
    sense = a[0] == '-'
    if w.isdigit():
      n = int(w)
      if sense: st = Puzzle(n)
      else: scale = n
    elif w == 'invert': invertp = sense
    elif w == 'reflect': reflectp = sense
    else: raise SyntaxError('unknown option')
  elif a[0] == '=': st.fill(a[1:])
  elif a[0] == ':': parse_moves(a[1:]).apply(st)
  elif a == '?': st.show()
  else:
    p = a
    _, e = OS.path.splitext(p)
    surf = CR.ImageSurface(CR.FORMAT_ARGB32, 100, 100)
    xf = Transform(TAU/10, TAU/10, 8)
    cr = CR.Context(surf)
    if invertp: cr.scale(scale, scale)
    else: cr.scale(scale, -scale)
    m = MeasuringPainter(xf, cr)
    st.render(m, reflectp)
    x0, y0 = cr.user_to_device(m.x0, m.y0)
    x1, y1 = cr.user_to_device(m.x1, m.y1)
    if x0 > x1: x0, x1 = x1, x0
    if y0 > y1: y0, y1 = y1, y0
    w, h = map(lambda x: int(M.ceil(x)), [x1 - x0, y1 - y0])
    if e == '.pdf': surf = CR.PDFSurface(p, w + 2, h + 2)
    elif e == '.png': surf = CR.ImageSurface(CR.FORMAT_ARGB32, w + 2, h + 2)
    elif e == '.ps':
      surf = CR.PSSurface(p, w + 2, h + 2)
      surf.set_eps(True)
    else: raise SyntaxError('unrecognized extension')
    cr = CR.Context(surf)
    cr.translate(1 - x0, 1 - y0)
    if invertp: cr.scale(scale, scale)
    else: cr.scale(scale, -scale)
    paint = Painter(xf, cr)
    st.render(paint, reflectp)
    surf.show_page()
    if e == '.png': surf.write_to_png(p)
    surf.finish()
    del cr, surf