Determine sRGB section in xyY space for specific Y value¶

I'm trying to determine the intersection of the sRGB "cube" in the xyY colorspace with a plane at a specific $Y$-value.

Find this project on GitHub.

When transformed to the xyY colorspace the RGB cube changes its shape:

In [1]:
Show code 
import colour
from colour import xyY_to_XYZ, xy_to_xyY
from colour.adaptation import matrix_chromatic_adaptation_VonKries
from colour.plotting import plot_RGB_colourspace_section
from colour.plotting import plot_RGB_colourspaces_gamuts

from IPython.display import display, IFrame, HTML, Latex, Markdown

import matplotlib.pyplot as plt

import sympy
from sympy import Matrix, Eq
from sympy import Max, Min, cancel, together, solve
from sympy import Piecewise, oo
from sympy import MatrixSymbol
from sympy.printing import latex as sympy_latex
from sympy.codegen.ast import Assignment
In [2]:
Show code 
with colour.utilities.suppress_warnings(colour_usage_warnings=True):
    plot_RGB_colourspaces_gamuts('sRGB', 'CIE xyY',segments=16)

Better Visualization on Wikipedia: https://en.wikipedia.org/wiki/File:SRGB_gamut_within_CIExyY_color_space_mesh.webm

Show Video from Wikipedia

I want to analytically calculate the shape of the section for a specific $Y$ value. The colour module can do it, but it uses the trimesh library to calculate the shape from a mesh. I want to calculate the points using math.

In [3]:
Show code 
with colour.utilities.suppress_warnings(colour_usage_warnings=True):
    plot_RGB_colourspace_section('sRGB', 'CIE xyY', Y = 0.3)

The interactive result¶

In [4]:
Show code 
display(IFrame('./tool.html', "100%", "900px"))

Requirements¶

for my fresh Miniconda3 Install:

conda create -n sRGB_Y_section_in_xyY_space python=3.9 numpy sympy networkx trimesh
activate sRGB_Y_section_in_xyY_space
conda install -c conda-forge colour-science trimesh
conda install -n sRGB_Y_section_in_xyY_space ipykernel --update-deps

Determining the transformation¶

We determine the $rgb$ to $xyY$ transformation using sympy

First we transform $rgb$ to $XYZ$

In [5]:
Show code 
sRGB = colour.models.rgb.RGB_COLOURSPACE_sRGB
r, g, b = sympy.symbols('r g b', real=True, positive=True)
x, y, Y = sympy.symbols('x y Y', real=True, positive=True)
v_rgb = Matrix([[r], [g], [b]])

# From colour.sRGB_to_XYZ
matrix_RGB_to_XYZ = Matrix(sRGB.matrix_RGB_to_XYZ)

M_CAT = matrix_chromatic_adaptation_VonKries(
    xyY_to_XYZ(xy_to_xyY(sRGB.whitepoint)),
    xyY_to_XYZ(xy_to_xyY(sRGB.whitepoint)),
    transform="CAT02",
)
M_CAT = Matrix(M_CAT)

v_xyz = matrix_RGB_to_XYZ * v_rgb
v_xyz = M_CAT * v_xyz
display(Latex(r"""$$
\begin{align*}
    \begin{pmatrix}X\\Y\\Z\end{pmatrix} &= M_{CAT} \cdot M_{RGB\rightarrow XYZ} \cdot \begin{pmatrix}r\\g\\b\end{pmatrix}\\""" + 
    f"&= {sympy_latex(M_CAT)} \cdot {sympy_latex(matrix_RGB_to_XYZ)}" + r"\cdot \begin{pmatrix}r\\g\\b\end{pmatrix} \\"
    f"&= {sympy_latex(M_CAT * matrix_RGB_to_XYZ)}" + r"\cdot \begin{pmatrix}r\\g\\b\end{pmatrix} \\"
    f"&= {sympy_latex(v_xyz)}" + 
r"""\end{align*}$$
"""))
$$ \begin{align*} \begin{pmatrix}X\\Y\\Z\end{pmatrix} &= M_{CAT} \cdot M_{RGB\rightarrow XYZ} \cdot \begin{pmatrix}r\\g\\b\end{pmatrix}\\&= \left[\begin{matrix}1.0 & -3.29597460435593 \cdot 10^{-17} & -2.77555756156289 \cdot 10^{-17}\\-1.49348849259878 \cdot 10^{-17} & 1.0 & 1.38777878078145 \cdot 10^{-17}\\-1.30104260698261 \cdot 10^{-18} & 0 & 1.0\end{matrix}\right] \cdot \left[\begin{matrix}0.4124 & 0.3576 & 0.1805\\0.2126 & 0.7152 & 0.0722\\0.0193 & 0.1192 & 0.9505\end{matrix}\right]\cdot \begin{pmatrix}r\\g\\b\end{pmatrix} \\&= \left[\begin{matrix}0.4124 & 0.3576 & 0.1805\\0.2126 & 0.7152 & 0.0722\\0.0193 & 0.1192 & 0.9505\end{matrix}\right]\cdot \begin{pmatrix}r\\g\\b\end{pmatrix} \\&= \left[\begin{matrix}0.1805 b + 0.3576 g + 0.4124 r\\0.0722 b + 0.7152 g + 0.2126 r\\0.9505 b + 0.1192 g + 0.0193 r\end{matrix}\right]\end{align*}$$

Then we transform $XYZ$ to $xy$ for fixed $Y$

$$ \begin{pmatrix}x\\y\\Y\end{pmatrix} = \begin{pmatrix}\frac{X}{X + Y +Z}\\\frac{Y}{X + Y + Z}\\Y\end{pmatrix} $$
In [6]:
Show code 
v_xyY = Matrix([
    [v_xyz[0] / (v_xyz[0] + v_xyz[1] + v_xyz[2])],
    [Y / (v_xyz[0] + v_xyz[1] + v_xyz[2])],
    [Y]
])

#v_xyY = v_xyY.subs(v_xyz[1], Y) #Maybe replace the Y value?

#We only care about xy, because Y will be given
v_xy = v_xyY[0:2,0]
Eq(Matrix([[x],[y]]),v_xy)
Out[6]:
$\displaystyle \left[\begin{matrix}x\\y\end{matrix}\right] = \left[\begin{matrix}\frac{0.1805 b + 0.3576 g + 0.4124 r}{1.2032 b + 1.192 g + 0.6443 r}\\\frac{Y}{1.2032 b + 1.192 g + 0.6443 r}\end{matrix}\right]$
In [7]:
Show code 
Eq(Y,v_xyz[1])
Out[7]:
$\displaystyle Y = 0.0722 b + 0.7152 g + 0.2126 r$

Calculate intersections of $RGB$-Cube in $xyY$ with $Y$-Plane¶

We calculate the intersections between all 6 sides of the $RGB$-Cube ($r=0,r=1,g=0,g=1,b=0,b=1$). Then we use the previous $Y$ equation to replace one of the remaining variables.

Because I don't know any better I will try all combinations of sides and variables:

In [8]:
Show code 
variations = [
    ((b,0), g, r), #Read as: set b=0, solve for g, remaining variable is r
        ((b,0), r, g),
    ((b,1), g, r),
        ((b,1), r, g),

    ((r,0), g, b),
        ((r,0), b, g),
    ((r,1), g, b),
        ((r,1), b, g),

    ((g,0), r, b),
        ((g,0), b, r),
    ((g,1), r, b),
        ((g,1), b, r),
]
In [9]:
Show code 
def show_step(doshow, step):
    if doshow:
        display(step())

results = {}
for i, (side_sub, solve_for, unsolved_for) in enumerate(variations):
    show_steps = i == 0
    show_step(show_steps, lambda: Latex(r'Starting with $$\begin{pmatrix}x\\y\end{pmatrix} = ' + f'{sympy_latex(v_xy)} \\text{{ with }} {sympy_latex(Eq(v_xyz[1],Y))}$$'))
    v_xy_side = v_xy.subs([side_sub])
    solve_eq = Eq(v_xyz[1],Y).subs([side_sub])
    show_step(show_steps, lambda: Latex(f'Setting ${side_sub[0]}={side_sub[1]}$: $${sympy_latex(v_xy_side)} \\text{{ with }} {sympy_latex(solve_eq)}$$'))
    solved_for_var = solve(solve_eq, solve_for)
    if (solved_for_var):
        show_step(show_steps, lambda: Latex(f'Solving ${sympy_latex(solve_eq)}$ for ${solve_for}$: $${sympy_latex(Eq(solve_for,solved_for_var[0]))} \\text{{ with }} g,r\\in [0,1]$$'))

        f = v_xy_side.subs(solve_for, solved_for_var[0])

        show_step(show_steps, lambda: Latex(f"""Inserting ${solve_for}$ into the $(x,y)$ yields:
        $$\\begin{{pmatrix}}x\\\\y\\end{{pmatrix}} = {sympy_latex(f)} $$
        """))

        f0 = solve(Eq(solved_for_var[0], 0), unsolved_for)[0]
        f1 = solve(Eq(solved_for_var[0], 1), unsolved_for)[0]

        min_other_var = together(cancel(Min(f0, f1)))
        max_other_var = together(cancel(Max(f0, f1)))

        show_step(show_steps, lambda: Latex(f"""To find the start and end of the intersection line, we solve for ${solve_for}=0$ and ${solve_for}=1$:
        $$\\begin{{align*}}
            {solve_for}=0 &\implies {unsolved_for} = {f0} \\\\
            {solve_for}=1 &\implies {unsolved_for} = {f1} \\\\
                &\implies {unsolved_for} \in C = [{sympy_latex(min_other_var)},{sympy_latex(max_other_var)}]
        \\end{{align*}}$$
        """))

        # Solution is only valid if other var is between 0 and 1
        condition_big = (min_other_var <= 1) & (max_other_var >= 0)
        condition = sympy.simplify(condition_big)

        show_step(show_steps, lambda: Latex(f"""If $C$ does not overlap with $[0,1]$ there is no intersection and no point can be calculated.
        $C$ overlaps with $[0,1]$ and an intersection can be calculted if
        $$\\begin{{align*}}
            &({sympy_latex(min_other_var)} <= 1) \\land ({sympy_latex(max_other_var)} >= 0) \\\\
            \implies&{sympy_latex(condition_big)} \\\\
            \implies& {sympy_latex(condition)} \\\\
        \\end{{align*}}$$
        """))

        # Alternatively one could use a Piecewise here and fold it later
        clamp0 = Max(0, min_other_var)
        clamp1 = Min(1, max_other_var)

        show_step(show_steps, lambda: Latex(f"""The start and end of the intersection line are
        $$ r={sympy_latex(clamp0)} \\text{{ and }} r={sympy_latex(clamp1)} $$
        """))

        p1 = together(f.subs(unsolved_for, clamp0))
        p2 = together(f.subs(unsolved_for, clamp1))

        show_step(show_steps, lambda: Latex(f"""This yields the start and end point of the intersection line:
        $$\\begin{{align*}}
            p_1 &= {sympy_latex(p1)} \\text{{ and }} \\\\
            p_2 &= {sympy_latex(p2)} \\\\
            &\\text{{ if }} {sympy_latex(condition)}
        \\end{{align*}}$$
        """))


        results[(side_sub, solve_for)] = [
             p1, p2, condition
        ]
Starting with $$\begin{pmatrix}x\\y\end{pmatrix} = \left[\begin{matrix}\frac{0.1805 b + 0.3576 g + 0.4124 r}{1.2032 b + 1.192 g + 0.6443 r}\\\frac{Y}{1.2032 b + 1.192 g + 0.6443 r}\end{matrix}\right] \text{ with } 0.0722 b + 0.7152 g + 0.2126 r = Y$$
Setting $b=0$: $$\left[\begin{matrix}\frac{0.3576 g + 0.4124 r}{1.192 g + 0.6443 r}\\\frac{Y}{1.192 g + 0.6443 r}\end{matrix}\right] \text{ with } 0.7152 g + 0.2126 r = Y$$
Solving $0.7152 g + 0.2126 r = Y$ for $g$: $$g = 1.39821029082774 Y - 0.297259507829978 r \text{ with } g,r\in [0,1]$$
Inserting $g$ into the $(x,y)$ yields: $$\begin{pmatrix}x\\y\end{pmatrix} = \left[\begin{matrix}\frac{0.5 Y + 0.3061 r}{1.66666666666667 Y + 0.289966666666667 r}\\\frac{Y}{1.66666666666667 Y + 0.289966666666667 r}\end{matrix}\right] $$
To find the start and end of the intersection line, we solve for $g=0$ and $g=1$: $$\begin{align*} g=0 &\implies r = 4.70366886171213*Y \\ g=1 &\implies r = 4.70366886171213*Y - 3.36406396989651 \\ &\implies r \in C = [4.70366886171213 Y - 3.36406396989651,4.70366886171213 Y] \end{align*}$$
If $C$ does not overlap with $[0,1]$ there is no intersection and no point can be calculated. $C$ overlaps with $[0,1]$ and an intersection can be calculted if $$\begin{align*} &(4.70366886171213 Y - 3.36406396989651 <= 1) \land (4.70366886171213 Y >= 0) \\ \implies&4.70366886171213 Y - 3.36406396989651 \leq 1 \\ \implies& Y \leq 0.927800000000001 \\ \end{align*}$$
The start and end of the intersection line are $$ r=\max\left(0, 4.70366886171213 Y - 3.36406396989651\right) \text{ and } r=\min\left(1, 4.70366886171213 Y\right) $$
This yields the start and end point of the intersection line: $$\begin{align*} p_1 &= \left[\begin{matrix}\frac{0.5 Y + 0.3061 \max\left(0, 4.70366886171213 Y - 3.36406396989651\right)}{1.66666666666667 Y + 0.289966666666667 \max\left(0, 4.70366886171213 Y - 3.36406396989651\right)}\\\frac{Y}{1.66666666666667 Y + 0.289966666666667 \max\left(0, 4.70366886171213 Y - 3.36406396989651\right)}\end{matrix}\right] \text{ and } \\ p_2 &= \left[\begin{matrix}\frac{0.5 Y + 0.3061 \min\left(1, 4.70366886171213 Y\right)}{1.66666666666667 Y + 0.289966666666667 \min\left(1, 4.70366886171213 Y\right)}\\\frac{Y}{1.66666666666667 Y + 0.289966666666667 \min\left(1, 4.70366886171213 Y\right)}\end{matrix}\right] \\ &\text{ if } Y \leq 0.927800000000001 \end{align*}$$

We repeat these steps for all sides and variables, which gives us 8 pairs of points that are corner points of the section - if their condition is valid.

Show solutions for different $Y$ values¶

Shown alongside the plot_RGB_colourspace_section result from the colour module, which uses a mesh-intersection algorithm to get the same result

In [10]:
Ys = [0.3, 0.5, 0.75, 0.9]
fig, axs = plt.subplots(2,2)

markers = {0: '+', 1:'x'}
colors = {b: 'blue', g: 'green', r:'yellow'}
linestyles = {r: 'solid', g: 'dashed', b: 'dotted'}

for Yset, ax in zip(Ys, axs.flat):
    print(f"Y = {Yset}:")
    plot_RGB_colourspace_section('sRGB', 'CIE xyY', origin=Yset, normalise=False, standalone=False, axes=ax)

    for (side_sub, solve_for),(p1, p2, condition) in results.items():
        print(f"{side_sub}, {solve_for}: ", end="")
        if not condition.subs(Y, Yset):
            print("Skipping! Outside interval")
            continue
        x1, y1 = p1.subs(Y, Yset)
        x2, y2 = p2.subs(Y, Yset)
        print(f"({x1}, {y1}) - ({x2}, {y2})")
        ax.plot([x1,x2],[y1,y2], label=f"{side_sub}, {solve_for}",
            marker = markers[side_sub[1]], color=colors[side_sub[0]],
            linestyle=linestyles[solve_for]
        )
plt.show()

Y = 0.3:
(b, 0), g: (0.300000000000000, 0.600000000000000) - (0.577366133592135, 0.379762859192371)
(b, 0), r: (0.577366133592134, 0.379762859192371) - (0.299999999999999, 0.600000000000001)
(b, 1), g: (0.185991660699996, 0.189529545550267) - (0.320637180742191, 0.160185102785441)
(b, 1), r: (0.320637180742191, 0.160185102785441) - (0.185991660699995, 0.189529545550267)
(r, 0), g: (0.300000000000000, 0.600000000000000) - (0.185991660699996, 0.189529545550267)
(r, 0), b: (0.185991660699996, 0.189529545550268) - (0.300000000000003, 0.600000000000009)
(r, 1), g: (0.577366133592135, 0.379762859192371) - (0.320637180742191, 0.160185102785441)
(r, 1), b: (0.320637180742192, 0.160185102785442) - (0.577366133592139, 0.379762859192375)
(g, 0), r: Skipping! Outside interval
(g, 0), b: Skipping! Outside interval
(g, 1), r: Skipping! Outside interval
(g, 1), b: Skipping! Outside interval
Y = 0.5:
(b, 0), g: (0.300000000000000, 0.600000000000000) - (0.495059200569750, 0.445117065788302)
(b, 0), r: (0.495059200569750, 0.445117065788303) - (0.299999999999999, 0.600000000000000)
(b, 1), g: (0.205824026719549, 0.260933096753992) - (0.317519075319181, 0.226637455616832)
(b, 1), r: (0.317519075319181, 0.226637455616832) - (0.205824026719549, 0.260933096753992)
(r, 0), g: (0.300000000000000, 0.600000000000000) - (0.205824026719549, 0.260933096753992)
(r, 0), b: (0.205824026719549, 0.260933096753994) - (0.300000000000002, 0.600000000000007)
(r, 1), g: (0.495059200569750, 0.445117065788302) - (0.317519075319181, 0.226637455616832)
(r, 1), b: (0.317519075319182, 0.226637455616833) - (0.495059200569754, 0.445117065788307)
(g, 0), r: Skipping! Outside interval
(g, 0), b: Skipping! Outside interval
(g, 1), r: Skipping! Outside interval
(g, 1), b: Skipping! Outside interval
Y = 0.75:
(b, 0), g: (0.327642853755811, 0.578050733894227) - (0.442282300482694, 0.487023528647806)
(b, 0), r: (0.442282300482694, 0.487023528647807) - (0.327642853755811, 0.578050733894227)
(b, 1), g: (0.222644528905781, 0.321492870002572) - (0.314735972548770, 0.285950308190888)
(b, 1), r: (0.314735972548770, 0.285950308190888) - (0.222644528905781, 0.321492870002572)
(r, 0), g: (0.250912034477843, 0.423265915111071) - (0.222644528905781, 0.321492870002572)
(r, 0), b: (0.222644528905781, 0.321492870002574) - (0.250912034477842, 0.423265915111070)
(r, 1), g: (0.442282300482694, 0.487023528647806) - (0.314735972548770, 0.285950308190888)
(r, 1), b: (0.314735972548771, 0.285950308190889) - (0.442282300482697, 0.487023528647811)
(g, 0), r: Skipping! Outside interval
(g, 0), b: Skipping! Outside interval
(g, 1), r: (0.327642853755811, 0.578050733894227) - (0.250912034477841, 0.423265915111068)
(g, 1), b: (0.250912034477842, 0.423265915111070) - (0.327642853755814, 0.578050733894233)
Y = 0.9:
(b, 0), g: (0.408706220886141, 0.513683956415632) - (0.422410100746755, 0.502802659267398)
(b, 0), r: (0.422410100746755, 0.502802659267398) - (0.408706220886142, 0.513683956415632)
(b, 1), g: (0.276461377969248, 0.328894161707499) - (0.313453617218774, 0.313279573011545)
(b, 1), r: (0.313453617218773, 0.313279573011545) - (0.276461377969248, 0.328894161707499)
(r, 0), g: Skipping! Outside interval
(r, 0), b: Skipping! Outside interval
(r, 1), g: (0.422410100746755, 0.502802659267398) - (0.313453617218774, 0.313279573011545)
(r, 1), b: (0.313453617218775, 0.313279573011547) - (0.422410100746758, 0.502802659267404)
(g, 0), r: Skipping! Outside interval
(g, 0), b: Skipping! Outside interval
(g, 1), r: (0.408706220886142, 0.513683956415632) - (0.276461377969248, 0.328894161707499)
(g, 1), b: (0.276461377969249, 0.328894161707502) - (0.408706220886145, 0.513683956415637)

Table showing results¶

From the text output it can be seen, that there are a lot of overlapping results.

For example the result for $f_{b=0,g}(0)$ is equal to $f_{b=0,r}(1)$.

These are the sets that grant equal results. Not all functions are valid for the whole $Y$-range, so the marked functions were manually selected to create a good coverage.

In [11]:
axes = {"x": 0, "y": 1}

allsets = [
    [
        ((b,0),g,0), ((b,0),r,1),
        ((r,0),g,0), ((r,0),b,1),
        ((g,1),r,0), ((g,1),b,1),
    ],
    [
        ((b,0),r,0), ((b,0),g,1),
        ((g,0),r,0), ((g,0),b,1),
        ((r,1),g,0), ((r,1),b,1),
    ],
    [
        ((b,1),r,0), ((b,1),g,1),
        ((g,0),b,0), ((g,0),r,1),
        ((r,1),b,0), ((r,1),g,1),
    ],
    [
        ((b,1),g,0), ((b,1),r,1),
        ((r,0),b,0), ((r,0),g,1),
        ((g,1),b,0), ((g,1),r,1),
    ],
]

selectedsets = [
    {
        ((b,0),g,0), ((g,1),r,0), 
    },
    {
        ((b,0),r,0), ((r,1),b,1),
    },
    {
        ((b,1),g,1), ((g,0),b,0),
    },
    {
        ((b,1),g,0), ((r,0),b,0), 
    },
]


Ys = [0, 0.01, 0.1, 0.2, 0.3, 0.7, 0.8, 0.9, 0.99,1] # 0.4, 0.5, 0.6 left out, because not close to any 
table_output = ["<table><tbody>"]
for i, (equalsets, selectedset) in enumerate(zip(allsets,selectedsets), 1):
    if i > 1:
        table_output.append(f'</tbody></table><table style="margin-top: 3em"><tbody>')
    table_output.append(f'<tr><th>Set {i}</th>')
    table_output.extend(f'<th>{Yset}</th>' for Yset in Ys)
    table_output.append('</tr>')
    for axisname, axis in axes.items():
        for s in equalsets:
            side_sub, solve_for, pnum = s
            if s in selectedset:
                table_output.append('<tr style="background-color: yellow;color: black">\n')
            else:
                table_output.append('<tr>\n')
            table_output.append(f'<th>\n$$f_{{{side_sub[0]}={side_sub[1]},{solve_for}}}({pnum}, {{{axisname}}})$$\n</th>')
            condition = results[(side_sub, solve_for)][2]
            for Yset in Ys:
                f = results[(side_sub, solve_for)][pnum][axis]
                f = Piecewise((f, condition), (oo, True))
                value = float(f.subs(Y, Yset))
                table_output.append(f'<td>{value:12.9f}</td>')
            table_output.append('</tr>\n')
table_output.append('</tbody></table>')

display(Markdown(''.join(table_output)))
Set 100.010.10.20.30.70.80.90.991
$$f_{b=0,g}(0, {x})$$ nan 0.300000000 0.300000000 0.300000000 0.300000000 0.300000000 0.360315454 0.408706221 inf inf
$$f_{b=0,r}(1, {x})$$ nan 0.300000000 0.300000000 0.300000000 0.300000000 0.300000000 0.360315454 0.408706221 inf inf
$$f_{r=0,g}(0, {x})$$ nan 0.300000000 0.300000000 0.300000000 0.300000000 0.300000000 inf inf inf inf
$$f_{r=0,b}(1, {x})$$ nan 0.300000000 0.300000000 0.300000000 0.300000000 0.300000000 inf inf inf inf
$$f_{g=1,r}(0, {x})$$ inf inf inf inf inf inf 0.360315454 0.408706221 0.322153757 inf
$$f_{g=1,b}(1, {x})$$ inf inf inf inf inf inf 0.360315454 0.408706221 0.322153757 inf
$$f_{b=0,g}(0, {y})$$ nan 0.600000000 0.600000000 0.600000000 0.600000000 0.600000000 0.552107696 0.513683956 inf inf
$$f_{b=0,r}(1, {y})$$ nan 0.600000000 0.600000000 0.600000000 0.600000000 0.600000000 0.552107696 0.513683956 inf inf
$$f_{r=0,g}(0, {y})$$ nan 0.600000000 0.600000000 0.600000000 0.600000000 0.600000000 inf inf inf inf
$$f_{r=0,b}(1, {y})$$ nan 0.600000000 0.600000000 0.600000000 0.600000000 0.600000000 inf inf inf inf
$$f_{g=1,r}(0, {y})$$ inf inf inf inf inf inf 0.552107696 0.513683956 0.344605315 inf
$$f_{g=1,b}(1, {y})$$ inf inf inf inf inf inf 0.552107696 0.513683956 0.344605315 inf
Set 200.010.10.20.30.70.80.90.991
$$f_{b=0,r}(0, {x})$$ nan 0.640074499 0.640074499 0.640074499 0.577366134 0.450422206 0.434978131 0.422410101 inf inf
$$f_{b=0,g}(1, {x})$$ nan 0.640074499 0.640074499 0.640074499 0.577366134 0.450422206 0.434978131 0.422410101 inf inf
$$f_{g=0,r}(0, {x})$$ nan 0.640074499 0.640074499 0.640074499 inf inf inf inf inf inf
$$f_{g=0,b}(1, {x})$$ nan 0.640074499 0.640074499 0.640074499 inf inf inf inf inf inf
$$f_{r=1,g}(0, {x})$$ inf inf inf inf 0.577366134 0.450422206 0.434978131 0.422410101 0.322153757 0.312715907
$$f_{r=1,b}(1, {x})$$ inf inf inf inf 0.577366134 0.450422206 0.434978131 0.422410101 0.322153757 inf
$$f_{b=0,r}(0, {y})$$ nan 0.329970511 0.329970511 0.329970511 0.379762859 0.480560196 0.492823261 0.502802659 inf inf
$$f_{b=0,g}(1, {y})$$ nan 0.329970511 0.329970511 0.329970511 0.379762859 0.480560196 0.492823261 0.502802659 inf inf
$$f_{g=0,r}(0, {y})$$ nan 0.329970511 0.329970511 0.329970511 inf inf inf inf inf inf
$$f_{g=0,b}(1, {y})$$ nan 0.329970511 0.329970511 0.329970511 inf inf inf inf inf inf
$$f_{r=1,g}(0, {y})$$ inf inf inf inf 0.379762859 0.480560196 0.492823261 0.502802659 0.344605315 0.329001481
$$f_{r=1,b}(1, {y})$$ inf inf inf inf 0.379762859 0.480560196 0.492823261 0.502802659 0.344605315 inf
Set 300.010.10.20.30.70.80.90.991
$$f_{b=1,r}(0, {x})$$ inf inf 0.182085716 0.269351508 0.320637181 0.315219531 0.314282195 0.313453617 0.312786018 inf
$$f_{b=1,g}(1, {x})$$ inf inf 0.182085716 0.269351508 0.320637181 0.315219531 0.314282195 0.313453617 0.312786018 0.312715907
$$f_{g=0,b}(0, {x})$$ nan 0.150016622 0.182085716 0.269351508 inf inf inf inf inf inf
$$f_{g=0,r}(1, {x})$$ nan 0.150016622 0.182085716 0.269351508 inf inf inf inf inf inf
$$f_{r=1,b}(0, {x})$$ inf inf inf inf 0.320637181 0.315219531 0.314282195 0.313453617 0.312786018 inf
$$f_{r=1,g}(1, {x})$$ inf inf inf inf 0.320637181 0.315219531 0.314282195 0.313453617 0.312786018 0.312715907
$$f_{b=1,r}(0, {y})$$ inf inf 0.077672922 0.125746040 0.160185103 0.275644812 0.295621112 0.313279573 0.327507306 inf
$$f_{b=1,g}(1, {y})$$ inf inf 0.077672922 0.125746040 0.160185103 0.275644812 0.295621112 0.313279573 0.327507306 0.329001481
$$f_{g=0,b}(0, {y})$$ nan 0.060006649 0.077672922 0.125746040 inf inf inf inf inf inf
$$f_{g=0,r}(1, {y})$$ nan 0.060006649 0.077672922 0.125746040 inf inf inf inf inf inf
$$f_{r=1,b}(0, {y})$$ inf inf inf inf 0.160185103 0.275644812 0.295621112 0.313279573 0.327507306 inf
$$f_{r=1,g}(1, {y})$$ inf inf inf inf 0.160185103 0.275644812 0.295621112 0.313279573 0.327507306 0.329001481
Set 400.010.10.20.30.70.80.90.991
$$f_{b=1,g}(0, {x})$$ inf inf 0.155578082 0.172574495 0.185991661 0.219778917 0.231176464 0.276461378 0.309419063 0.312715907
$$f_{b=1,r}(1, {x})$$ inf inf 0.155578082 0.172574495 0.185991661 0.219778917 0.231176464 0.276461378 0.309419063 inf
$$f_{r=0,b}(0, {x})$$ nan 0.150016622 0.155578082 0.172574495 0.185991661 0.219778917 inf inf inf inf
$$f_{r=0,g}(1, {x})$$ nan 0.150016622 0.155578082 0.172574495 0.185991661 0.219778917 inf inf inf inf
$$f_{g=1,b}(0, {x})$$ inf inf inf inf inf inf 0.231176464 0.276461378 0.309419063 inf
$$f_{g=1,r}(1, {x})$$ inf inf inf inf inf inf 0.231176464 0.276461378 0.309419063 inf
$$f_{b=1,g}(0, {y})$$ inf inf 0.080029878 0.141222991 0.189529546 0.311175651 0.328760112 0.328894162 0.328991721 0.329001481
$$f_{b=1,r}(1, {y})$$ inf inf 0.080029878 0.141222991 0.189529546 0.311175651 0.328760112 0.328894162 0.328991721 inf
$$f_{r=0,b}(0, {y})$$ nan 0.060006649 0.080029878 0.141222991 0.189529546 0.311175651 inf inf inf inf
$$f_{r=0,g}(1, {y})$$ nan 0.060006649 0.080029878 0.141222991 0.189529546 0.311175651 inf inf inf inf
$$f_{g=1,b}(0, {y})$$ inf inf inf inf inf inf 0.328760112 0.328894162 0.328991721 inf
$$f_{g=1,r}(1, {y})$$ inf inf inf inf inf inf 0.328760112 0.328894162 0.328991721 inf

Sadly I don't know of a better way to do this instead of the manual sorting and guessing. If you know: Let me know.

Resulting points¶

We combine the selected results into a function for each point. Showing $p_0$ as an example.

In [12]:
def generate_functions_for_set(s):
    axes = {"x", "y"}
    result = []
    inf_matrix = Matrix([[oo], [oo]])
    for i, equalset in enumerate(s):
        fs = sympy.piecewise_exclusive(Piecewise(*[
            (
                results[(side_sub, solve_for)][pnum],
                results[(side_sub, solve_for)][2],
            )
            for side_sub, solve_for, pnum in equalset
        ], (inf_matrix, True)))
        #Readd the "otherwise" case that was removed by piecewise_exclusive
        fs = Piecewise(*[(arg[0], sympy.simplify(arg[1])) for arg in fs.args],(inf_matrix, True))
        result.append(fs)
    return result
point_functions = generate_functions_for_set(selectedsets)

display(Latex(f'$$p_{0} = {sympy.printing.latex(point_functions[0])}$$'))
$$p_0 = \begin{cases} \left[\begin{matrix}\frac{1.93979303857008 Y + 0.0404469426152398 \max\left(0, 13.8504155124654 Y - 12.8504155124654\right) - 1.02973998118532}{3.03057384760113 Y + 0.984392568203199 \max\left(0, 13.8504155124654 Y - 12.8504155124654\right) - 0.975466415804328}\\\frac{Y}{3.03057384760113 Y + 0.984392568203199 \max\left(0, 13.8504155124654 Y - 12.8504155124654\right) - 0.975466415804328}\end{matrix}\right] & \text{for}\: Y \geq 0.715199999999999 \wedge Y \leq 0.999999999999999 \\\left[\begin{matrix}\frac{0.5 Y + 0.3061 \max\left(0, 4.70366886171213 Y - 3.36406396989651\right)}{1.66666666666667 Y + 0.289966666666667 \max\left(0, 4.70366886171213 Y - 3.36406396989651\right)}\\\frac{Y}{1.66666666666667 Y + 0.289966666666667 \max\left(0, 4.70366886171213 Y - 3.36406396989651\right)}\end{matrix}\right] & \text{for}\: Y \leq 0.927800000000001 \wedge \left(Y > 0.999999999999999 \vee Y < 0.715199999999999\right) \\\left[\begin{matrix}\infty\\\infty\end{matrix}\right] & \text{otherwise} \end{cases}$$

Generated code¶

In [13]:
for pf, p_mat_sym in zip(point_functions, (MatrixSymbol(f'p{i}', 2, 1) for i in range(4))):
    refactored_function = Piecewise(*[(Assignment(p_mat_sym, arg[0]), arg[1]) for arg in pf.args])
    code = sympy.jscode(refactored_function)
    display(Markdown(f'```javascript\n{code}\n```'))

if (Y >= 0.715199999999999 && Y <= 0.999999999999999) {
   p0[0] = (1.93979303857008*Y + 0.0404469426152398*Math.max(0, 13.8504155124654*Y - 12.8504155124654) - 1.02973998118532)/(3.03057384760113*Y + 0.984392568203199*Math.max(0, 13.8504155124654*Y - 12.8504155124654) - 0.975466415804328);
   p0[1] = Y/(3.03057384760113*Y + 0.984392568203199*Math.max(0, 13.8504155124654*Y - 12.8504155124654) - 0.975466415804328);
}
else if (Y <= 0.927800000000001 && (Y > 0.999999999999999 || Y < 0.715199999999999)) {
   p0[0] = (0.5*Y + 0.3061*Math.max(0, 4.70366886171213*Y - 3.36406396989651))/(1.66666666666667*Y + 0.289966666666667*Math.max(0, 4.70366886171213*Y - 3.36406396989651));
   p0[1] = Y/(1.66666666666667*Y + 0.289966666666667*Math.max(0, 4.70366886171213*Y - 3.36406396989651));
}
else {
   p0[0] = Number.POSITIVE_INFINITY;
   p0[1] = Number.POSITIVE_INFINITY;
}

if (Y >= 0.2126 && Y <= 0.999999999999998) {
   p1[0] = (2.5*Y - 1.4304*Math.min(1, 1.39821029082774*Y - 0.297259507829978) - 0.1191)/(16.6648199445983*Y - 10.7266792243767*Math.min(1, 1.39821029082774*Y - 0.297259507829978) - 2.89864072022161);
   p1[1] = Y/(16.6648199445983*Y - 10.7266792243767*Math.min(1, 1.39821029082774*Y - 0.297259507829978) - 2.89864072022161);
}
else if (Y <= 0.927799999999999 && (Y > 0.999999999999998 || Y < 0.2126)) {
   p1[0] = (1.93979303857008*Y - 1.02973998118532*Math.max(0, 1.39821029082774*Y - 0.297259507829978))/(3.03057384760113*Y - 0.975466415804328*Math.max(0, 1.39821029082774*Y - 0.297259507829978));
   p1[1] = Y/(3.03057384760113*Y - 0.975466415804328*Math.max(0, 1.39821029082774*Y - 0.297259507829978));
}
else {
   p1[0] = Number.POSITIVE_INFINITY;
   p1[1] = Number.POSITIVE_INFINITY;
}

if (Y <= 0.2848) {
   p2[0] = (2.5*Y - 0.1191*Math.max(0, 4.70366886171214*Y - 0.339604891815616))/(16.6648199445983*Y - 2.89864072022161*Math.max(0, 4.70366886171214*Y - 0.339604891815616));
   p2[1] = Y/(16.6648199445983*Y - 2.89864072022161*Math.max(0, 4.70366886171214*Y - 0.339604891815616));
}
else if (Y <= 1.0) {
   p2[0] = (0.5*Y + 0.3061*Math.min(1, 4.70366886171213*Y - 0.339604891815616) + 0.1444)/(1.66666666666667*Y + 0.289966666666667*Math.min(1, 4.70366886171213*Y - 0.339604891815616) + 1.08286666666667);
   p2[1] = Y/(1.66666666666667*Y + 0.289966666666667*Math.min(1, 4.70366886171213*Y - 0.339604891815616) + 1.08286666666667);
}
else {
   p2[0] = Number.POSITIVE_INFINITY;
   p2[1] = Number.POSITIVE_INFINITY;
}

if (Y >= 0.0722000000000001 && Y <= 1.0) {
   p3[0] = (0.5*Y + 0.3061*Math.max(0, 4.70366886171213*Y - 3.70366886171213) + 0.1444)/(1.66666666666667*Y + 0.289966666666667*Math.max(0, 4.70366886171213*Y - 3.70366886171213) + 1.08286666666667);
   p3[1] = Y/(1.66666666666667*Y + 0.289966666666667*Math.max(0, 4.70366886171213*Y - 3.70366886171213) + 1.08286666666667);
}
else if (Y <= 0.787399999999999 && (Y > 1.0 || Y < 0.0722000000000001)) {
   p3[0] = (2.5*Y - 1.4304*Math.max(0, 1.39821029082774*Y - 0.100950782997763))/(16.6648199445983*Y - 10.7266792243767*Math.max(0, 1.39821029082774*Y - 0.100950782997763));
   p3[1] = Y/(16.6648199445983*Y - 10.7266792243767*Math.max(0, 1.39821029082774*Y - 0.100950782997763));
}
else {
   p3[0] = Number.POSITIVE_INFINITY;
   p3[1] = Number.POSITIVE_INFINITY;
}

Lazy code generation¶

As an alternative to the manual selection one could also just use all the functions and let sympy.piecewise_exclusive do the job. That does not get rid of manually grouping the functions into groups, though, and the generated functions are quite elaborate:

In [14]:
point_functions_unfiltered = generate_functions_for_set(allsets)

display(Latex(f'$$p_{0} = {sympy.printing.latex(point_functions_unfiltered[0])}$$'))
$$p_0 = \begin{cases} \left[\begin{matrix}\frac{0.5 Y + 0.3061 \max\left(0, 4.70366886171213 Y - 3.36406396989651\right)}{1.66666666666667 Y + 0.289966666666667 \max\left(0, 4.70366886171213 Y - 3.36406396989651\right)}\\\frac{Y}{1.66666666666667 Y + 0.289966666666667 \max\left(0, 4.70366886171213 Y - 3.36406396989651\right)}\end{matrix}\right] & \text{for}\: Y \leq 0.927800000000001 \\\left[\begin{matrix}\frac{1.93979303857008 Y + 0.0404469426152398 \max\left(0, 13.8504155124654 Y - 12.8504155124654\right) - 1.02973998118532}{3.03057384760113 Y + 0.984392568203199 \max\left(0, 13.8504155124654 Y - 12.8504155124654\right) - 0.975466415804328}\\\frac{Y}{3.03057384760113 Y + 0.984392568203199 \max\left(0, 13.8504155124654 Y - 12.8504155124654\right) - 0.975466415804328}\end{matrix}\right] & \text{for}\: Y \leq 0.999999999999999 \wedge Y \neq 0.787399999999999 \wedge Y \neq 0.7874 \wedge Y \neq 0.927799999999999 \\\left[\begin{matrix}\frac{2.5 Y - 0.1191 \min\left(1, 4.70366886171214 Y - 3.36406396989652\right) - 1.4304}{16.6648199445983 Y - 2.89864072022161 \min\left(1, 4.70366886171214 Y - 3.36406396989652\right) - 10.7266792243767}\\\frac{Y}{16.6648199445983 Y - 2.89864072022161 \min\left(1, 4.70366886171214 Y - 3.36406396989652\right) - 10.7266792243767}\end{matrix}\right] & \text{for}\: Y \leq 0.999999999999998 \wedge \left(Y > 0.999999999999999 \vee Y < 0.715199999999999\right) \wedge Y \neq 0.787399999999999 \wedge Y \neq 0.7874 \wedge Y \neq 0.927799999999999 \\\left[\begin{matrix}\infty\\\infty\end{matrix}\right] & \text{otherwise} \end{cases}$$
In [ ]:
for pf, p_mat_sym in zip(point_functions_unfiltered, (MatrixSymbol(f'p{i}', 2, 1) for i in range(4))):
    refactored_function = Piecewise(*[(Assignment(p_mat_sym, arg[0]), arg[1]) for arg in pf.args])
    code = sympy.jscode(refactored_function)
    display(Markdown(f'```javascript\n{code}\n```'))