Quick and Dirty X3D edition tools.

x3ditor.py

import math
import numpy
import BeautifulSoup  # pip install --user BeautifulSoup

# The purpose of this lib is to fix X3D files exported via the `vtkExporter`,
# for example using mayavi, in order to feed them to `x3dom.js`.
# Tasks :
#   - `add_axes()`
#     Create axes with labels and ticks (works best if you have an outline)

# We need separate methods for separate tasks.
# Example tasks still to do -- no deadline :
# - Set / Remove Background
# - Set / Remove Lights
# - Set embedded viewpoint to use an orthographic projection


class BeautifulX3dSoup(BeautifulSoup.BeautifulStoneSoup):
    """
    Configure BeautifulSoup for the X3D syntax.
    """
    def __init__(self, *args, **kwargs):
        if not kwargs.has_key('smartQuotesTo'):
            kwargs['smartQuotesTo'] = self.HTML_ENTITIES  # is this relevant?
        kwargs['isHTML'] = False
        BeautifulSoup.BeautifulStoneSoup.__init__(self, *args, **kwargs)

    SELF_CLOSING_TAGS = BeautifulSoup.buildTagMap(None, (
        'meta',
        'background',
        'proximitysensor',
        'directionallight',
        'navigationinfo',
        'route',
        'fontstyle',
        'viewpoint',
        'normal',
        'color',
        'coordinate',
        'material',
    ))


# TASKS ########################################################################

# todo: Pipe this into a __main__, for CLI scripting
# Well, right now we're debugging, so it's here, import it from your own python

# TASK : CREATE AXES ###########################################################

def add_axes(input_x3d, output_x3d=None,
             x_min=None, x_max=None, x_txt="X",
             y_min=None, y_max=None, y_txt="Y",
             z_min=None, z_max=None, z_txt="Z",
             x_val_min=None, x_val_max=None,
             y_val_min=None, y_val_max=None,
             z_val_min=None, z_val_max=None,
             ticks_cnt=None, ticks_fmt=None,):
    """
    TASK : CREATE AXES

    This adds text objects to the x3d file at the end of the scene, the labels
    and ticks of the axes.

    The `x_min`, `x_man`, etc. are the coordinates extrema of the ouline, in the
    local space of the scene. They will be used to position the labels and ticks
    in the scene, along the axes of the outline.

    The `x_val_min`, `x_val_max`, etc. are the displayed extrema of the ticks.
    You can set these values to the numbers of your choosing, and they will be
    defaulted to their respective coordinates values when not set.

    input_x3d: str
        Filepath to the input x3ds file.
        /!\ The input x3d file will be mutated if `output_x3d` is not set!
    output_x3d: str
        Optional filepath to a different output path, so as to leave the input
        file unchanged.
    x_txt: str | None
        The (single-line) text label showed on the X axes.
        When set to None, no label will be shown.
    y_txt: str | None
        The (single-line) text label showed on the Y axes.
        When set to None, no label will be shown.
    z_txt: str | None
        The (single-line) text label showed on the Z axes.
        When set to None, no label will be shown.
    ticks_cnt: int | int[3]
        How many ticks per axis, 5 on each of them by default.
        When zero, does not show ticks for this axis.
        When invalid (strictly negative or 1), this defaults to 2.
        You can customize each axis independently by specifying a 3-tuple or a
        list in the `numpy` order : `(z, y, x)`.
    ticks_fmt: str | str[3]
        The format of the ticks per axis, in `%`-compatible syntax.
        It is defaulted to `%.2f`.
        You can customize each axis independently by specifying a 3-tuple or a
        list in the `numpy` order : `(z, y, x)`.

    """

    DEFAULT_TICKS_CNT = 5
    DEFAULT_TICKS_FMT = "%.2f"

    X_AXIS = (1., 0., 0.)
    Y_AXIS = (0., 1., 0.)
    Z_AXIS = (0., 0., 1.)

    FULL_TURN = 6.283185307179586  # math.pi * 2
    HALF_TURN = FULL_TURN / 2.
    QUAR_TURN = HALF_TURN / 2.
    TRQU_TURN = HALF_TURN + QUAR_TURN

    # OPEN THE X3D WITH XML SCRAPER ############################################

    with open(input_x3d) as x3d_file:
        soup = BeautifulX3dSoup(markup=x3d_file)

    scene = soup.find("scene")

    # SANITIZE #################################################################

    output_x3d = output_x3d if output_x3d is not None else input_x3d

    x_txt = _escape_html_attribute(x_txt)
    y_txt = _escape_html_attribute(y_txt)
    z_txt = _escape_html_attribute(z_txt)

    # todo: guess extrema from X3D coordinates tags' attributes?
    x_min = x_min if x_min is not None else 0.0
    x_max = x_max if x_max is not None else 1.0
    y_min = y_min if y_min is not None else 0.0
    y_max = y_max if y_max is not None else 1.0
    z_min = z_min if z_min is not None else 0.0
    z_max = z_max if z_max is not None else 1.0
    
    x_val_min = x_val_min if x_val_min is not None else x_min
    y_val_min = y_val_min if y_val_min is not None else y_min
    z_val_min = z_val_min if z_val_min is not None else z_min
    x_val_max = x_val_max if x_val_max is not None else x_max
    y_val_max = y_val_max if y_val_max is not None else y_max
    z_val_max = z_val_max if z_val_max is not None else z_max

    x_mid = (x_max - x_min) / 2.0
    y_mid = (y_max - y_min) / 2.0
    z_mid = (z_max - z_min) / 2.0

    ticks_cnt = ticks_cnt if ticks_cnt is not None else DEFAULT_TICKS_CNT
    if not hasattr(ticks_cnt, '__len__') or isinstance(ticks_cnt, basestring):
        ticks_cnt = (int(ticks_cnt), int(ticks_cnt), int(ticks_cnt))
    if not len(ticks_cnt) == 3:
        raise ValueError("Provide one tick count or three ticks counts.")

    ticks_cnt = ((ticks_cnt[0] if ticks_cnt[0] > 1 or ticks_cnt[0] == 0 else 2),
                 (ticks_cnt[1] if ticks_cnt[1] > 1 or ticks_cnt[1] == 0 else 2),
                 (ticks_cnt[2] if ticks_cnt[2] > 1 or ticks_cnt[2] == 0 else 2))

    z_tik, y_tik, x_tik = ticks_cnt

    ticks_fmt = ticks_fmt if ticks_fmt is not None else DEFAULT_TICKS_FMT
    if not hasattr(ticks_fmt, '__len__') or isinstance(ticks_fmt, basestring):
        ticks_fmt = (ticks_fmt, ticks_fmt, ticks_fmt)
    if not len(ticks_fmt) == 3:
        raise ValueError("Provide one tick format or three ticks formats.")

    z_fmt, y_fmt, x_fmt = ticks_fmt

    x_stp = (x_max - x_min) / (x_tik-1)
    y_stp = (y_max - y_min) / (y_tik-1)
    z_stp = (z_max - z_min) / (z_tik-1)

    x_val_stp = (x_val_max - x_val_min) / (x_tik-1)
    y_val_stp = (y_val_max - y_val_min) / (y_tik-1)
    z_val_stp = (z_val_max - z_val_min) / (z_tik-1)

    s_tik = 6  # todo tick text size as input?

    # Text is faced by default like so, in a right-handed 3D referential :
    #     Y
    #     ^         (Z points towards the eye)
    #     |BLABLA
    #     +---> X

    # todo
    # It would be ideal to have a look_at(direction, up) in text_node() to
    # automatically compute the rotation axis and angle from default orientation
    # but the setting of the angle becomes tricky.

    # Show labels
    # This is the BAD part.
    # The position and rotation of the labels has been hand-crafted.
    # I'm not satisfied by it : it is messy and not elegant. But it works.

    # Pads for the labels, to put them below the ticks
    # Only some (about 2) axes require special padding.
    # todo: figure out how to guess the paddings
    pad = 4.
    sad = pad * 0.25  # small pad
    bad = pad * 1.25  # big pad

    # Lots of hard-code below.
    # It positions and rotates the labels and ticks, relatively to their initial
    # position. It's hard to maintain, and not ideal. But it works. (kinda)
    # This may be dried up using config matrices, and a smarter text_node()

    # Eye direction is the direction the eye of the viewer is facing,
    # and the axis for which we create a label and ticks.

    # Eye direction: -Z, Axis: X
    if x_txt is not None:
        scene.append(text_node(x_txt, x=x_mid, y=y_min-pad))
    for i in range(x_tik):
        # Pull the extremums to prevent ticks from different axes to overlap.
        p = x_min+i*x_stp if i > 0 else x_min+sad
        p = p if i < x_tik-1 else x_max-sad
        # Display value is not position value, so we format it.
        v = x_fmt % (x_val_min+i*x_val_stp)
        scene.append(text_node(v, size=s_tik, x=p))
    # Eye direction: -Z, Axis: Y
    ax, an = [(Z_AXIS,), (QUAR_TURN,)]
    if y_txt is not None:
        scene.append(text_node(y_txt, y=y_mid, x=x_min-pad, axis=ax, angle=an))
    for i in range(y_tik):
        p = y_min+i*y_stp if i > 0 else y_min+sad
        p = p if i < y_tik-1 else y_max-sad
        v = y_fmt % (y_val_min+i*y_val_stp)
        scene.append(text_node(v, size=s_tik, y=p, axis=ax, angle=an))
    # Eye direction: +Z, Axis X
    ax, an = [(Y_AXIS,), (HALF_TURN,)]
    if x_txt is not None:
        scene.append(text_node(x_txt, x=x_mid, y=y_min-pad, z=z_max, axis=ax, angle=an))
    for i in range(x_tik):
        p = x_min+i*x_stp if i > 0 else x_min+sad
        p = p if i < x_tik-1 else x_max-sad
        v = x_fmt % (x_val_min+i*x_val_stp)
        scene.append(text_node(v, size=s_tik, x=p, z=z_max, axis=ax, angle=an))
    # Eye direction: +Z, Axis Y
    ax, an = [(Y_AXIS, Z_AXIS), (HALF_TURN, QUAR_TURN)]
    if y_txt is not None:
        scene.append(text_node(y_txt, y=y_mid, z=z_max, x=x_max + bad, axis=ax, angle=an))
    for i in range(y_tik):
        p = y_min+i*y_stp if i > 0 else y_min+sad
        p = p if i < y_tik-1 else y_max-sad
        v = y_fmt % (y_val_min+i*y_val_stp)
        scene.append(text_node(v, size=s_tik, y=p, z=z_max, x=x_max + sad, axis=ax, angle=an))
    # Eye direction: -X, Axis Y
    ax, an = [(Y_AXIS, Z_AXIS), (QUAR_TURN, QUAR_TURN)]
    if y_txt is not None:
        scene.append(text_node(y_txt, y=y_mid, z=z_max+bad, axis=ax, angle=an))
    for i in range(y_tik):
        p = y_min+i*y_stp if i > 0 else y_min+sad
        p = p if i < y_tik-1 else y_max-sad
        v = y_fmt % (y_val_min+i*y_val_stp)
        scene.append(text_node(v, size=s_tik, y=p, z=z_max+sad, axis=ax, angle=an))
    # Eye direction: -X, Axis Z
    ax, an = [(Y_AXIS,), (QUAR_TURN,)]
    if z_txt is not None:
        scene.append(text_node(z_txt, z=z_mid, y=y_min-pad, axis=ax, angle=an))
    for i in range(z_tik):
        p = z_min+i*z_stp if i > 0 else z_min+sad
        p = p if i < z_tik-1 else z_max-sad
        v = z_fmt % (z_val_min+i*z_val_stp)
        scene.append(text_node(v, size=s_tik, z=p, axis=ax, angle=an))
    # Eye direction: +X, Axis Y
    ax, an = [(X_AXIS, Y_AXIS), (TRQU_TURN, TRQU_TURN)]
    if y_txt is not None:
        scene.append(text_node(y_txt, y=y_mid, x=x_max, z=z_min-pad, axis=ax, angle=an))
    for i in range(y_tik):
        p = y_min+i*y_stp if i > 0 else y_min+sad
        p = p if i < y_tik-1 else y_max-sad
        v = y_fmt % (y_val_min+i*y_val_stp)
        scene.append(text_node(v, size=s_tik, y=p, x=x_max, axis=ax, angle=an))
    # Eye direction: +X, Axis Y
    ax, an = [(Y_AXIS,), (TRQU_TURN,)]
    if z_txt is not None:
        scene.append(text_node(z_txt, z=z_mid, x=x_max, y=y_min-pad, axis=ax, angle=an))
    for i in range(z_tik):
        p = z_min+i*z_stp if i > 0 else z_min+sad
        p = p if i < z_tik-1 else z_max-sad
        v = z_fmt % (z_val_min+i*z_val_stp)
        scene.append(text_node(v, size=s_tik, z=p, x=x_max, axis=ax, angle=an))
    # Eye direction: -Y, Axis X
    ax, an = [(X_AXIS,), (-QUAR_TURN,)]
    if x_txt is not None:
        scene.append(text_node(x_txt, x=x_mid, y=y_max, z=z_min-pad, axis=ax, angle=an))
    for i in range(x_tik):
        p = x_min+i*x_stp if i > 0 else x_min+sad
        p = p if i < x_tik-1 else x_max-sad
        v = x_fmt % (x_val_min+i*x_val_stp)
        scene.append(text_node(v, size=s_tik, x=p, y=y_max, axis=ax, angle=an))
    # Eye direction: -Y, Axis Z
    ax, an = [(Y_AXIS, X_AXIS), (QUAR_TURN, -QUAR_TURN)]
    if z_txt is not None:
        scene.append(text_node(z_txt, z=z_mid, y=y_max, x=x_min-pad, axis=ax, angle=an))
    for i in range(z_tik):
        p = z_min+i*z_stp if i > 0 else z_min+sad
        p = p if i < z_tik-1 else z_max-sad
        v = z_fmt % (z_val_min+i*z_val_stp)
        scene.append(text_node(v, size=s_tik, z=p, y=y_max, axis=ax, angle=an))
    # Eye direction: +Y, Axis X
    ax, an = [(X_AXIS,), (QUAR_TURN,)]
    if x_txt is not None:
        scene.append(text_node(x_txt, x=x_mid, y=y_max, z=z_min-pad, axis=ax, angle=an))
    for i in range(x_tik):
        p = x_min+i*x_stp if i > 0 else x_min+sad
        p = p if i < x_tik-1 else x_max-sad
        v = x_fmt % (x_val_min+i*x_val_stp)
        scene.append(text_node(v, size=s_tik, x=p, y=y_max, axis=ax, angle=an))
    # Eye direction: +Y, Axis Z
    ax, an = [(Y_AXIS, X_AXIS), (-QUAR_TURN, QUAR_TURN)]
    if z_txt is not None:
        scene.append(text_node(z_txt, z=z_mid, y=y_max, x=x_min-pad, axis=ax, angle=an))
    for i in range(z_tik):
        p = z_min+i*z_stp if i > 0 else z_min+sad
        p = p if i < z_tik-1 else z_max-sad
        v = z_fmt % (z_val_min+i*z_val_stp)
        scene.append(text_node(v, size=s_tik, z=p, y=y_max, axis=ax, angle=an))
    # OK WE ARE DONE WITH THE LABELS AND TICKS! -- phew...

    # SAVE !
    with open(output_x3d, "w") as x3d_file:
        x3d_file.write(soup.prettify())


# UTILS ########################################################################

X3D_TEXT_NODE = """
<Transform translation="{x} {y} {z}" scale="0.2 0.2 0.2" rotation="{rotationAxis} {rotationAngle}">
  <Shape>
    <Appearance>
      <Material diffuseColor="{diffuseColor}" emissiveColor="{emissiveColor}"/>
    </Appearance>
    <Text string="{title}">
      <FontStyle family='"TYPEWRITER"' horizontal="true" topToBottom="true" justify='"MIDDLE" "MIDDLE"' size="{size}"/>
    </Text>
  </Shape>
</Transform>
"""

def text_node(title,
              size=10,
              x=0., y=0., z=0.,
              color_diffuse=(0., 0., 0.),
              color_emissive=(0., 0., 0.),
              axis=None,   # 3-Tuple, only one or a list
              angle=None,  # float64, only one or a list
              ):  # #  :(
    """
    Create a BeautifulSoup-compatible text node to embed in a X3D file.
    When the axis looks at you, the rotation specified by the angle is CCW.
    """
    if axis is None:
        axis = (0.0, 0.0, 0.0)
    if angle is None:
        angle = 0.0
    axis, angle = _combine(axis, angle)
    return BeautifulX3dSoup(markup=X3D_TEXT_NODE.format(
        title=title,
        size=size,
        x=str(x), y=str(y), z=str(z),
        diffuseColor=' '.join([str(v) for v in color_diffuse]),
        emissiveColor=' '.join([str(v) for v in color_emissive]),
        rotationAxis=' '.join([str(v) for v in axis]),
        rotationAngle=angle,
    ))


# PRIVATE UTILS ################################################################

def _escape_html_attribute(s):
    # There are probably plenty of other things to escape. Use a lib instead?
    return s.replace('"', '&#34;')


def _combine(axes, angles):
    """
    Reduce two (index-mapped) lists of axes and angles describing rotations,
    into the resulting axis and angle of the sequentially combined rotations.

    This leverages quaternion multiplication awesomeness.

    The axes are 3-Tuples, and the angles are floats.
    Provide either two lists of them, or two of them.

    Note: the lists order is important ! This operation is not commutative !
    Note: we make no input sanity checks whatsoever, so this method is private.
    Note: if you can ditch `transformations` and do this with sympy, do it !
    """
    if not hasattr(angles, '__len__'):
        return axes, angles

    count = len(angles)
    combined = quaternion_from_matrix(identity_matrix())
    for i in range(count):
        combined = quaternion_multiply(
            combined,
            quaternion_about_axis(angles[i], axes[i])
        )
    angle, axis, _ = rotation_from_matrix(quaternion_matrix(combined))

    return axis, angle


# EMBEDDED THIRD-PARTY LIBS ####################################################

# Homogeneous Transformation Matrices and Quaternions.
# The full library is available here :
# http://www.lfd.uci.edu/~gohlke/code/transformations.py.html

# Copyright (c) 2006-2015, Christoph Gohlke
# Copyright (c) 2006-2015, The Regents of the University of California
# Produced at the Laboratory for Fluorescence Dynamics
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are met:
#
# * Redistributions of source code must retain the above copyright
#   notice, this list of conditions and the following disclaimer.
# * Redistributions in binary form must reproduce the above copyright
#   notice, this list of conditions and the following disclaimer in the
#   documentation and/or other materials provided with the distribution.
# * Neither the name of the copyright holders nor the names of any
#   contributors may be used to endorse or promote products derived
#   from this software without specific prior written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
# ARE DISCLAIMED.  IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
# LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
# CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
# SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
# INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
# CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
# POSSIBILITY OF SUCH DAMAGE.

# from transformations import identity_matrix
# from transformations import quaternion_about_axis as quat
# from transformations import quaternion_multiply
# from transformations import quaternion_from_matrix
# from transformations import rotation_from_matrix
# from transformations import quaternion_matrix

# epsilon for testing whether a number is close to zero
_EPS = numpy.finfo(float).eps * 4.0


def vector_norm(data, axis=None, out=None):
    """
    Return length, i.e. Euclidean norm, of ndarray along axis.

    >>> v = numpy.random.random(3)
    >>> n = vector_norm(v)
    >>> numpy.allclose(n, numpy.linalg.norm(v))
    True
    >>> v = numpy.random.rand(6, 5, 3)
    >>> n = vector_norm(v, axis=-1)
    >>> numpy.allclose(n, numpy.sqrt(numpy.sum(v*v, axis=2)))
    True
    >>> n = vector_norm(v, axis=1)
    >>> numpy.allclose(n, numpy.sqrt(numpy.sum(v*v, axis=1)))
    True
    >>> v = numpy.random.rand(5, 4, 3)
    >>> n = numpy.empty((5, 3))
    >>> vector_norm(v, axis=1, out=n)
    >>> numpy.allclose(n, numpy.sqrt(numpy.sum(v*v, axis=1)))
    True
    >>> vector_norm([])
    0.0
    >>> vector_norm([1])
    1.0
    """
    data = numpy.array(data, dtype=numpy.float64, copy=True)
    if out is None:
        if data.ndim == 1:
            return math.sqrt(numpy.dot(data, data))
        data *= data
        out = numpy.atleast_1d(numpy.sum(data, axis=axis))
        numpy.sqrt(out, out)
        return out
    else:
        data *= data
        numpy.sum(data, axis=axis, out=out)
        numpy.sqrt(out, out)


def identity_matrix():
    """
    Return 4x4 identity/unit matrix.

    >>> I = identity_matrix()
    >>> numpy.allclose(I, numpy.dot(I, I))
    True
    >>> numpy.sum(I), numpy.trace(I)
    (4.0, 4.0)
    >>> numpy.allclose(I, numpy.identity(4))
    True
    """
    return numpy.identity(4)


def quaternion_about_axis(angle, axis):
    """
    Return quaternion for rotation about axis.

    >>> q = quaternion_about_axis(0.123, [1, 0, 0])
    >>> numpy.allclose(q, [0.99810947, 0.06146124, 0, 0])
    True
    """
    q = numpy.array([0.0, axis[0], axis[1], axis[2]])
    qlen = vector_norm(q)
    if qlen > _EPS:
        q *= math.sin(angle/2.0) / qlen
    q[0] = math.cos(angle/2.0)
    return q


def quaternion_multiply(quaternion1, quaternion0):
    """
    Return multiplication of two quaternions.

    >>> q = quaternion_multiply([4, 1, -2, 3], [8, -5, 6, 7])
    >>> numpy.allclose(q, [28, -44, -14, 48])
    True
    """
    w0, x0, y0, z0 = quaternion0
    w1, x1, y1, z1 = quaternion1
    return numpy.array([-x1*x0 - y1*y0 - z1*z0 + w1*w0,
                         x1*w0 + y1*z0 - z1*y0 + w1*x0,
                        -x1*z0 + y1*w0 + z1*x0 + w1*y0,
                         x1*y0 - y1*x0 + z1*w0 + w1*z0], dtype=numpy.float64)


def quaternion_from_matrix(matrix, isprecise=False):
    """
    Return quaternion from rotation matrix.

    If isprecise is True, the input matrix is assumed to be a precise rotation
    matrix and a faster algorithm is used.

    >>> q = quaternion_from_matrix(numpy.identity(4), True)
    >>> numpy.allclose(q, [1, 0, 0, 0])
    True
    >>> q = quaternion_from_matrix(numpy.diag([1, -1, -1, 1]))
    >>> numpy.allclose(q, [0, 1, 0, 0]) or numpy.allclose(q, [0, -1, 0, 0])
    True
    """
    m = numpy.array(matrix, dtype=numpy.float64, copy=False)[:4, :4]
    if isprecise:
        q = numpy.empty((4, ))
        t = numpy.trace(m)
        if t > m[3, 3]:
            q[0] = t
            q[3] = m[1, 0] - m[0, 1]
            q[2] = m[0, 2] - m[2, 0]
            q[1] = m[2, 1] - m[1, 2]
        else:
            i, j, k = 1, 2, 3
            if m[1, 1] > m[0, 0]:
                i, j, k = 2, 3, 1
            if m[2, 2] > m[i, i]:
                i, j, k = 3, 1, 2
            t = m[i, i] - (m[j, j] + m[k, k]) + m[3, 3]
            q[i] = t
            q[j] = m[i, j] + m[j, i]
            q[k] = m[k, i] + m[i, k]
            q[3] = m[k, j] - m[j, k]
        q *= 0.5 / math.sqrt(t * m[3, 3])
    else:
        m00 = m[0, 0]
        m01 = m[0, 1]
        m02 = m[0, 2]
        m10 = m[1, 0]
        m11 = m[1, 1]
        m12 = m[1, 2]
        m20 = m[2, 0]
        m21 = m[2, 1]
        m22 = m[2, 2]
        # symmetric matrix K
        K = numpy.array([[m00-m11-m22, 0.0,         0.0,         0.0],
                         [m01+m10,     m11-m00-m22, 0.0,         0.0],
                         [m02+m20,     m12+m21,     m22-m00-m11, 0.0],
                         [m21-m12,     m02-m20,     m10-m01,     m00+m11+m22]])
        K /= 3.0
        # quaternion is eigenvector of K that corresponds to largest eigenvalue
        w, V = numpy.linalg.eigh(K)
        q = V[[3, 0, 1, 2], numpy.argmax(w)]
    if q[0] < 0.0:
        numpy.negative(q, q)
    return q


def rotation_from_matrix(matrix):
    """
    Return rotation angle and axis from rotation matrix.
    """
    r = numpy.array(matrix, dtype=numpy.float64, copy=False)
    r33 = r[:3, :3]
    # direction: unit eigenvector of r33 corresponding to eigenvalue of 1
    w, ww = numpy.linalg.eig(r33.T)
    i = numpy.where(abs(numpy.real(w) - 1.0) < 1e-8)[0]
    if not len(i):
        raise ValueError("no unit eigenvector corresponding to eigenvalue 1")
    direction = numpy.real(ww[:, i[-1]]).squeeze()
    # point: unit eigenvector of r33 corresponding to eigenvalue of 1
    w, q = numpy.linalg.eig(r)
    i = numpy.where(abs(numpy.real(w) - 1.0) < 1e-8)[0]
    if not len(i):
        raise ValueError("no unit eigenvector corresponding to eigenvalue 1")
    point = numpy.real(q[:, i[-1]]).squeeze()
    point /= point[3]
    # rotation angle depending on direction
    cosa = (numpy.trace(r33) - 1.0) / 2.0
    if abs(direction[2]) > 1e-8:
        sina = (r[1, 0] + (cosa-1.0)*direction[0]*direction[1]) / direction[2]
    elif abs(direction[1]) > 1e-8:
        sina = (r[0, 2] + (cosa-1.0)*direction[0]*direction[2]) / direction[1]
    else:
        sina = (r[2, 1] + (cosa-1.0)*direction[1]*direction[2]) / direction[0]
    angle = math.atan2(sina, cosa)
    return angle, direction, point



def quaternion_matrix(quaternion):
    """
    Return homogeneous rotation matrix from quaternion.

    >>> M = quaternion_matrix([1, 0, 0, 0])
    >>> numpy.allclose(M, numpy.identity(4))
    True
    >>> M = quaternion_matrix([0, 1, 0, 0])
    >>> numpy.allclose(M, numpy.diag([1, -1, -1, 1]))
    True
    """
    q = numpy.array(quaternion, dtype=numpy.float64, copy=True)
    n = numpy.dot(q, q)
    if n < _EPS:
        return numpy.identity(4)
    q *= math.sqrt(2.0 / n)
    q = numpy.outer(q, q)
    return numpy.array([
        [1.0-q[2, 2]-q[3, 3],     q[1, 2]-q[3, 0],     q[1, 3]+q[2, 0], 0.0],
        [    q[1, 2]+q[3, 0], 1.0-q[1, 1]-q[3, 3],     q[2, 3]-q[1, 0], 0.0],
        [    q[1, 3]-q[2, 0],     q[2, 3]+q[1, 0], 1.0-q[1, 1]-q[2, 2], 0.0],
        [                0.0,                 0.0,                 0.0, 1.0]])

just an FYI but this borks under BeautifulSoup4's many, many non-backward compatible changes...