Film Like Curve Bezier Stuffs

stuffs.py

#!/usr/bin/env python
# -*- coding: utf-8 -*-

import numpy
import math


# Return value as scalar or array. Shamelessly borrowed from
# Colour-science.org.
def as_numeric(obj, as_type=numpy.float64):
    try:
        return as_type(obj)
    except TypeError:
        return obj


def shape_OCIO_matrix(numpy_matrix):
    # Shape the RGB to XYZ array for OpenColorIO
    ocio_matrix = numpy.pad(
        numpy_matrix,
        [(0, 1), (0, 1)],
        mode='constant'
    )
    ocio_matrix = ocio_matrix.flatten()
    ocio_matrix[-1] = 1.

    return ocio_matrix


# Convert relative EV to radiometric linear value.
def calculate_ev_to_rl(
    in_ev,
    rl_middle_grey=0.18
):
    in_ev = numpy.asarray(in_ev)

    return as_numeric(numpy.power(2., in_ev) * rl_middle_grey)


# Convert radiometric linear value to relative EV
def calculate_rl_to_ev(
    in_rl,
    rl_middle_grey=0.18
):
    in_rl = numpy.asarray(in_rl)

    return as_numeric(numpy.log2(in_rl) - numpy.log2(rl_middle_grey))


# Calculate the Y intercept based on slope and an
# X / Y coordinate pair.
def calculate_line_y_intercept(in_x, in_y, slope):
    # y = mx + b, b = y - mx
    return (in_y - (slope * in_x))


# Calculate the Y of a line given by slope and X coordinate.
def calculate_line_y(in_x, y_intercept, slope):
    # y = mx + b
    return (slope * in_x) + y_intercept


# Calculate the X of a line given by the slope and Y intercept.
def calculate_line_x(in_y, y_intercept, slope):
    # y = mx + b, x = (y - b) / m
    return (in_y - y_intercept) / slope


# Calculate the slope of a line given by the Y intercept
# and an X / Y coordinate pair.
def calculate_line_slope(in_x, in_y, y_intercept):
    # y = mx + b, m = (y - b) / x
    return (in_y - y_intercept) / in_x


# Calculate the slope of a line given by two coordinates.
def calculate_line_slope(in_x1, in_y1, in_x2, in_y2):
    # x = x1 - x2, y = y1 - y2
    # y = mx + b, m = (y - b) / x
    return (in_x1 - in_x2) / (in_y1 - in_y2)


# Calculate the slope of a line given by the coordinates of the line's
# X / Y coordinate pairs.
def calculate_line_slope_coord(a_xy, b_xy):
    difference = a_xy - b_xy
    return difference[1] / difference[0]


# Calculate the control points for an S quadratic curve given four
# sets of values that dictate the minimum, low break point, high break point,
# and maximum. This is a workhorse tool for calculating quadratic, single
# tension curves.
#
# Parameters include:
#   minimum:
#       The minimum range of the total curve inputs, as specified as an input /
#       output pair.
#   pivot:
#       The fulcrum point for the linear portion of the curve, specified as an
#       input / output pair.
#   maximum:
#       The maximum range of the total curve inputs, as specified as an input /
#       output pair.
#   slope_degrees:
#       The slope of the linear portion of the curve, given in degrees where
#       valid values are between 0 and 90 degrees.
#   tension_low:
#       Controls the tension of the curved portion of the curve below
#       the low breakpoint into the linear portion. Maintains the slope of the
#       angle of the linear portion.
#   tension_high:
#       Controls the tension of the curved portion of the curve above
#       the high breakpoint into the linear portion. Maintains the slope of the
#       angle of the linear portion.

def calculate_control_points_from_descriptor(descriptor=None):
    try:
        minimum = numpy.asarray(descriptor["minimum"])
        if not(
                (minimum.dtype == numpy.float64) and
                (minimum.size == 2) and
                numpy.all(minimum > 0.0)
        ):
            raise ValueError("Minimum does not meet critera.")
        maximum = numpy.asarray(descriptor["maximum"])
        if not(
                (maximum.dtype == numpy.float64) and
                (maximum.size == 2)
        ):
            raise ValueError("Maximum does not meet critera.")
        pivot = numpy.asarray(descriptor["pivot"])
        if not(
                (pivot.dtype == numpy.float64) and
                (pivot.size == 2) and
                numpy.all(pivot < maximum)
        ):
            raise ValueError("Pivot does not meet critera.")
        slope_degrees = numpy.asarray(descriptor["slope_degrees"])
        if not(
                (slope_degrees.dtype == numpy.float64) and
                (slope_degrees.size == 1) and
                (slope_degrees >= 45.0) and
                (slope_degrees <= 90.0)
        ):
            raise ValueError("Slope does not meet critera.")
        tension_low = numpy.asarray(descriptor["tension_low"])
        if not(
                (tension_low.dtype == numpy.float64) and
                (tension_low.size == 1) and
                (tension_low > 0.0) and
                (tension_low < 1.0)
        ):
            raise ValueError("Low tension does not meet critera.")
        tension_high = numpy.asarray(descriptor["tension_high"])
        if not(
                (tension_high.dtype == numpy.float64) and
                (tension_high.size == 1) and
                (tension_high > 0.0) and
                (tension_high < 1.0)
        ):
            raise ValueError("High tension does not meet critera.")

    except Exception:
        raise

    control_points = numpy.zeros((3, 3, 2))
    input_linear_range = numpy.zeros(2)
    input_breakpoints = numpy.zeros(2)

    total_range = (maximum - minimum)
    total_ev = (
        calculate_rl_to_ev(total_range, pivot) -
        calculate_rl_to_ev(minimum, pivot)
    )

    print("TOTAL_EV: {}".format(total_ev))

    # Calculate the linear range from an EV based value, and calculate
    # the output from it.
    #
    # TODO: Remove the hard coded linear range and offset, and substitue in
    # variables via the descriptor.

    # Base linear range off of approximately 70% of the total EVs covered.
    linear_ev = 0.7 * total_ev

    # Arbitrarily centre the offset at 50%.
    offset_ev = 0.5 * linear_ev

    print("OFFSET_EV: {}".format(offset_ev))
    # Calculate the low input offset in terms of linear nits.
    input_breakpoints[0] = calculate_ev_to_rl(
        0.0 - offset_ev[0], pivot[0]
    )

    # Calculate the high input offset in terms of linear nits.
    input_breakpoints[1] = calculate_ev_to_rl(
        0.0 + offset_ev[1], pivot[0]
    )

    slope = math.tan(math.radians(slope_degrees))

    # Calculate the Y intercept given the slope and the pivot.
    y_intercept = calculate_line_y_intercept(pivot[0], pivot[1], slope)
    print("PIVOT: {}, Y_INTERCEPT: {}".format(pivot, y_intercept))

    print("Y_INTERCEPT: {}".format(y_intercept))
    control_low = [
        calculate_line_x(minimum[1], y_intercept, slope),
        minimum[1]
    ]

    control_high = [
        calculate_line_x(maximum[1], y_intercept, slope),
        maximum[1]
    ]
    print("CONTROL_LOW: {}, CONTROL_HIGH: {}".format(
        control_low, control_high))

    control_points = numpy.array([
        [
            minimum,
            control_low + (1.0 - tension_low) * (low_break -
                                                 control_low),
            low_break
        ],
        [
            low_break,
            pivot,
            high_break
        ],
        [
            high_break,
            control_high - (1.0 - tension_high) * (control_high -
                                                   high_break),
            maximum
        ]
    ])

    return control_points


def is_valid_control_points_quadratic(control_points):
    if control_points is None:
        return False
    else:
        control_points = numpy.asarray(control_points)
        # Total number of sets of control points test
        if control_points.shape[0] != 3:
            raise ValueError(
                "Invalid number of sections for the quadratic curve.\n"
                "\tQuadratic curves must have three sections."
            )

        # Total number of control points per set
        if control_points.shape[1] != 3:
            raise ValueError(
                "Invalid number of control points for a quadratic curve.\n"
                "\tQuadratic curves require three control points per "
                "section."
            )

        # Total number of coordinates per control point
        if control_points.shape[2] != 2:
            raise ValueError(
                "Invalid number of coordinates for a quadratic curve.\n"
                "\tQuadratic curves require two coordinates per point."
            )
    return True


def calculate_t_from_x_quadratic(in_x_values, control_points=None):
    total_entries = None
    if numpy.isscalar(in_x_values):
        total_entries = 1

    in_x_values = numpy.asarray(in_x_values)

    if total_entries is None:
        total_entries = in_x_values.shape[0]

    out_t = numpy.zeros(total_entries)

    try:
        if not is_valid_control_points_quadratic(control_points):
            control_points = calculate_control_points_from_descriptor()
    except Exception:
        raise

    # The following assumes that there is exactly one y solution for
    # each x input, which obviously is not the case for all Bezier
    # curves, and only works for the limited case of Bezier usage here.
    #
    # While there are several means to solve this, we will lean on the
    # handy functions in numpy to solve our coefficients. That is, we
    # have known x and y coordinates for each Bezier section, and we need
    # to isolate the t coefficient, then solve.
    #
    # First up, we take the basic Bezier quadratic formula, and multiply
    # it out.
    #
    # Basic Quadratic Bezier formula is:
    # x = (1. - t)**2 * cP0 + 2 * (1. - t) * t * cP1 + t**2 * cP2
    #
    # Equivalent to:
    # x =
    #   (1. - 2. * t + t**2) * cP0 +
    #   2. * t * cP1 - 2. * t**2 * cP1 +
    #   t**2 * cP2
    #
    # Equivalent to:
    # x =
    #   (cP0) - (2. * t * cP0) + (t**2 * cP0) +
    #   (2. * t * cP1) - (2. * t**2 * cP1) +
    #   (t**2 * cP2)
    #
    # Grouping together by orders of t:
    # x =
    #   (t**2 + cP0) - (t**2 * 2. * cP1) + (t**2 * cP2) -
    #   (t**1 * 2. * cP1) - (t**1 * 2. * cP0) +
    #   (t**0 * cP0)
    #
    # Set equal to zero and isolate out our t coefficients:
    # 0 =
    #   t**2 * (cP0 - (2. * cP1) + cP2) +
    #   t**1 * (2. * cP1) - (2. * cP0) +
    #   t**0 * (cP0 - x)
    #
    # Which gives us our coefficients to solve:
    # t_2 = (cP0 - 2. * cP1 + cP2)
    # t_1 = (2. * cP1) - (2. * cP0)
    # t_0 = (cP0 - x)
    for index in range(total_entries):
        for current_control_point in control_points:
            if current_control_point[0][0] <= in_x_values.item(index) \
               <= current_control_point[2][0]:
                coefficients = [
                    current_control_point[0][0] -
                    (2. * current_control_point[1][0]) +
                    current_control_point[2][0],
                    (2. * current_control_point[1][0]) -
                    (2. * current_control_point[0][0]),
                    current_control_point[0][0] - in_x_values.item(index)
                ]

                roots = numpy.roots(coefficients)
                correct_root = None

                for root in roots:
                    if numpy.isreal(root) and (0. <= root < 10000.):
                        correct_root = root
                    elif in_x_values.item(index) == 1.:
                        correct_root = 1.

                if correct_root is not None:
                    out_t[index] = correct_root
                else:
                    raise ValueError(
                        "No valid root found for coefficients. Invalid curve "
                        "or input x value."
                    )
    return as_numeric(out_t)


def calculate_y_from_x_quadratic(
        in_x_values,
        in_t_values,
        control_points=None):
    total_entries = None
    if numpy.isscalar(in_x_values):
        total_entries = 1

    in_x_values = numpy.asarray(in_x_values)

    if total_entries is None:
        total_entries = in_x_values.shape[0]

    out_y = numpy.zeros(total_entries)

    try:
        if not is_valid_control_points_quadratic(control_points):
            control_points = calculate_control_points_from_descriptor()
    except Exception:
        raise

    for index in range(total_entries):
        for current_control_point in control_points:
            if current_control_point[0][0] <= in_x_values[index] \
               <= current_control_point[2][0]:

                # TODO: Vectorize these coefficients to remove the for and
                # if statements above?
                out_y[index] = (
                    ((1. - in_t_values[index])**2. *
                     current_control_point[0][1]) +
                    (2. * (1. - in_t_values[index]) * in_t_values[index] *
                     current_control_point[1][1]) +
                    (in_t_values[index]**2. * current_control_point[2][1]))

    return as_numeric(out_y)


def calculate_slope_from_xy_quadratic(
        in_xy_values,
        in_t_values,
        control_points):
    total_entries = None
    if numpy.isscalar(in_xy_values):
        total_entries = 1

    in_xy_values = numpy.asarray(in_xy_values)

    if total_entries is None:
        total_entries = in_xy_values.shape[0]

    out_slopes = numpy.zeros(total_entries)
    out_slopes_xy = numpy.zeros((total_entries, 2))

    try:
        if not is_valid_control_points_quadratic(control_points):
            control_points = calculate_control_points_from_descriptor()
    except Exception:
        raise

    # The first derivative, aka "slope" of the quadratic curve function
    # with respect to t is given by:
    # B'(t) = 2 * (1. - t) * (cP1 - cP0) + 2 * t * (cP2 - cP1)
    for index in range(total_entries):
        for current_control_point in control_points:
            if current_control_point[0, 0] <= in_xy_values[index, 0] \
               <= current_control_point[2, 0]:

                # TODO: Vectorize these coefficients to remove the for and
                # if statements above?
                out_slopes_xy[index] = (
                    2.0 * (1.0 - in_t_values[index]) * (
                        current_control_point[1] -
                        current_control_point[0]
                    ) + (2.0 * in_t_values[index] * (
                        current_control_point[2] -
                        current_control_point[1])
                    )
                )

                out_slopes = numpy.divide(
                    out_slopes_xy[:, 1], out_slopes_xy[:, 0],
                    where=out_slopes_xy[:, 0] != 0)
    return as_numeric(out_slopes)


def is_valid_control_points(control_points):
    if control_points is None:
        return False

    control_points = numpy.asarray(control_points)

    try:
        # Assert shape is roughly the shape of a graph curve.
        if not len(control_points.shape) == 3:
            raise ValueError(
                "Malformed number of entries to define a graph curve.\n"
                "\tNumber of nested entries: {}".format(
                    len(control_points.shape))
            )
        # Total number of sets of control points test, one minimum.
        if not control_points.shape[0] >= 1:
            raise ValueError(
                "Not enough control points to define a graph curve.\n"
                "\tNumber of control points: {}".format(
                    control_points.shape[0])
            )

        # Assert count of of control point section. A section must consist of
        # a degree and minimum two coordinates at a minimum.
        if not control_points.shape[1] >= 3:
            raise ValueError(
                "Not enough entries to define a section of a graph curve.\n"
                "\tNumber of entries: {}".format(control_points.shape[1])
            )

        # Assert the degree is within the current limits of OpenAgX. Degrees
        # supported are 1 and 2 currently.
        if not ((control_points.shape[1, 0] >= 1) and
                (control_points.shape[1, 0] >= 1)):
            raise ValueError(
                "Degree of graph curve section not supported.\n"
                "\tDegree: {}".format(control_points.shape[1, 0])
            )

        # Total number of coordinates per control point.
        if control_points.shape[2] != 2:
            raise ValueError(
                "One or more coordinate is not a coordinate pair.\n"
                "\tCoordinates: {}".format(control_points.shape[2])
            )
    except ValueError:
        raise
        # return False
    except IndexError as ie:
        raise IndexError(
            "Invalid number of entries to define a graph curve.\n"
            "\t{}".format(ie)
        )

    return True


def calculate_t_from_x_arbitrary_quadratic(in_x_values, control_points=None):
    total_entries = None
    if numpy.isscalar(in_x_values):
        total_entries = 1

    in_x_values = numpy.asarray(in_x_values)

    if total_entries is None:
        total_entries = in_x_values.shape[0]

    out_t = numpy.zeros(total_entries)

    try:
        if not is_valid_control_points(control_points):
            control_points = calculate_control_points_from_descriptor()
    except Exception:
        raise

    # The following assumes that there is exactly one y solution for
    # each x input, which obviously is not the case for all Bezier
    # curves, and only works for the limited case of Bezier usage here.
    #
    # While there are several means to solve this, we will lean on the
    # handy functions in numpy to solve our coefficients. That is, we
    # have known x and y coordinates for each Bezier section, and we need
    # to isolate the t coefficient, then solve.
    #
    # First up, we take the basic Bezier quadratic formula, and multiply
    # it out.
    #
    # Basic Quadratic Bezier formula is:
    # x = (1. - t)**2 * cP0 + 2 * (1. - t) * t * cP1 + t**2 * cP2
    #
    # Equivalent to:
    # x =
    #   (1. - 2. * t + t**2) * cP0 +
    #   2. * t * cP1 - 2. * t**2 * cP1 +
    #   t**2 * cP2
    #
    # Equivalent to:
    # x =
    #   (cP0) - (2. * t * cP0) + (t**2 * cP0) +
    #   (2. * t * cP1) - (2. * t**2 * cP1) +
    #   (t**2 * cP2)
    #
    # Grouping together by orders of t:
    # x =
    #   (t**2 + cP0) - (t**2 * 2. * cP1) + (t**2 * cP2) -
    #   (t**1 * 2. * cP1) - (t**1 * 2. * cP0) +
    #   (t**0 * cP0)
    #
    # Set equal to zero and isolate out our t coefficients:
    # 0 =
    #   t**2 * (cP0 - (2. * cP1) + cP2) +
    #   t**1 * (2. * cP1) - (2. * cP0) +
    #   t**0 * (cP0 - x)
    #
    # Which gives us our coefficients to solve:
    # t_2 = (cP0 - 2. * cP1 + cP2)
    # t_1 = (2. * cP1) - (2. * cP0)
    # t_0 = (cP0 - x)
    for index in range(total_entries):
        for current_control_point in control_points:
            if current_control_point[0][0] <= in_x_values.item(index) \
               <= current_control_point[2][0]:
                coefficients = [
                    current_control_point[0][0] -
                    (2. * current_control_point[1][0]) +
                    current_control_point[2][0],
                    (2. * current_control_point[1][0]) -
                    (2. * current_control_point[0][0]),
                    current_control_point[0][0] - in_x_values.item(index)
                ]

                roots = numpy.roots(coefficients)
                correct_root = None

                for root in roots:
                    if numpy.isreal(root) and (0. <= root < 1.):
                        correct_root = root
                    elif in_x_values.item(index) == 1.:
                        correct_root = 1.

                if correct_root is not None:
                    out_t[index] = correct_root
                else:
                    raise ValueError(
                        "No valid root found for coefficients. Invalid curve "
                        "or input x value."
                    )
    return as_numeric(out_t)


def calculate_gamut_map_colourimetric(
        source_colourspace,
        destination_colourspace):
    return False