Cubic root solving

Update: Autodesk released their interpolation code for Maya animation curves, weighted tangents on animation curves do exactly this.
Refer to (and use!) https://github.com/Autodesk/animx instead of the code below. I know it’s not python, but it does work where I found bugs in my version below.

I really need to get back to this, but I mocked up this bit of code and finally got it to work.
It is about animation curve evaluation, due to the parametric nature of curves it is very complicated to get the value of a curve on at an arbitrary time, and it involves finding the parameter T for a value X.

This demo shows how the parameter (left to right) relates to the time of the animcurve (top to bottom), so the red curve shows linear parameter increasing results in non linear time samples.
It then computes backwards from the found value to the parameter again, showing we can successfully reconstruct the linear input given the time output.

I intend to clean up this code and make a full 2D example, rendering a 2D cubic spline segment as normal and then overlaying an evaluation based on the X coordinate, but wanted to dump the result nonetheless. Knowing how bad I am at getting back to things…
Using QT purely for demonstration, code itself is pure python…

from PyQt4.QtCore import *
from PyQt4.QtGui import *


def cubicArgs(x0, x1, x2, x3):
    a = x3 + (x1 - x2) * 3.0 - x0
    b = 3.0 * (x2 - 2.0 * x1 + x0)
    c = 3.0 * (x1 - x0)
    d = x0
    return a, b, c, d


def cubicEvaluate(x0, x1, x2, x3, p):
    # convert points to a cubic function & evaluate at p
    a, b, c, d = cubicArgs(x0, x1, x2, x3)
    return a * p * p * p + b * p * p + c * p + d


class CurveDebug(QWidget):
    def __init__(self):
        super(CurveDebug, self).__init__()
        self.t = QTimer()
        self.t.timeout.connect(self.update)
        self.t.start(16)
        self.ot = time.time()

    def paintEvent(self, event):
        painter = QPainter(self)

        life = time.time() - self.ot

        padding = 100
        w = self.width() - padding * 2
        h = self.height() - padding * 2
        painter.translate(padding, padding)

        # zig zag 2D bezier
        x0, y0 = 0, 0
        x1, y1 = (sin(life) * 0.5 + 0.5) * w, 0
        x2, y2 = 0, h
        x3, y3 = w, h

        # draw hull
        # painter.setPen(QColor(100, 220, 220))
        # painter.drawLine(x0, y0, x1, y1)
        # painter.drawLine(x1, y1, x2, y2)
        # painter.drawLine(x2, y2, x3, y3)

        for i in xrange(w):
            p = i / float(w - 1)

            # draw curve
            # painter.setPen(QColor(220, 100, 220))
            # x = cubicEvaluate(x0, x1, x2, x3, p)
            # y = cubicEvaluate(y0, y1, y2, y3, p)
            # painter.drawPoint(x, y)

            # draw X as function of P
            painter.setPen(QColor(220, 100, 100))
            x = cubicEvaluate(x0, x1, x2, x3, p)
            painter.drawPoint(i, x)

            # now let's evaluate the curve at x and see if we can get the original p back
            # make cubic with offset
            a, b, c, d = cubicArgs(x0 - x, x1 - x, x2 - x, x3 - x)

            # find roots
            # http://www.1728.org/cubic2.htm
            f = ((3.0 * c / a) - ((b * b) / (a * a))) / 3.0
            g = (((2.0 * b * b * b) / (a * a * a)) - ((9.0 * b * c) / (a * a)) + ((27.0 * d) / a)) / 27.0
            _h = ((g * g) / 4.0) + ((f * f * f) / 27.0)
            root0, root1, root2 = None, None, None
            if _h <= 0.0:
                # we have 3 real roots
                if f == 0 and g == 0:
                    # all roots are real & equal
                    _i = d / a
                    root0 = -copysign(pow(abs(_i), 0.3333), _i)
                else:
                    _i = sqrt((g * g / 4.0) - _h)
                    j = pow(_i, 0.3333333)
                    k = acos(-(g / (2.0 * _i)))
                    m = cos(k / 3.0)
                    n = sqrt(3.0) * sin(k / 3.0)
                    _p = b / (3.0 * a)
                    root0 = 2.0 * j * m - _p
                    root1 = -j * (m + n) - _p
                    root2 = -j * (m - n) - _p
            else:
                # we have only 1 real root
                R = -(g / 2.0) + sqrt(_h)
                S = copysign(pow(abs(R), 0.3333333), R)
                T = -(g / 2.0) - sqrt(_h)
                U = copysign(pow(abs(T), 0.3333333), T)
                root0 = (S + U) - (b / (3.0 * a))

            painter.setPen(QColor(100, 100, 220))
            painter.drawPoint(i, root0 * h)
            if root1:
                painter.drawPoint(i, root1 * h)
                painter.drawPoint(i, root2 * h)


app = QApplication([])
cvd = CurveDebug()
cvd.show()
cvd.resize(300, 300)
app.exec_()

Leave a Reply

Your email address will not be published. Required fields are marked *