Skip to content

Instantly share code, notes, and snippets.

@remia

remia/ocio_icc.py

Created September 27, 2022 18:02
Show Gist options

  • Select an option

    Learn more about clone URLs
  • Save remia/4cab0fb18a794d699d8ae2d6b77091a4 to your computer and use it in GitHub Desktop.

Select an option

Learn more about clone URLs
Save remia/4cab0fb18a794d699d8ae2d6b77091a4 to your computer and use it in GitHub Desktop.
Download ZIP
ICC and ParametricCurve tests in Python
Raw
ocio_icc.py
import os
import re
import subprocess as sp
import tempfile
from warnings import warn
import matplotlib
matplotlib.use('Qt5Agg')
import matplotlib.pyplot as plt
import numpy as np
import colour
import PyOpenColorIO as OCIO
# print(f"OCIO {OCIO.__version__}")
#
# ICC Utilities
#
# Conversions to / from ICC s15Fixed16Number
def float_to_s15f16(x):
return np.around(x * 2**16).astype(np.int32)
def s15f16_to_float(x):
return x / 2**16
def s15f16_roundtrip(x):
return s15f16_to_float(float_to_s15f16(x))
def quantize_params(params):
""" For Python side computation, make sure we simulate imprecisions due
to the s15f16 encoding used for paraCurve parameters.
"""
for name, value in params.items():
if name != "t":
params[name] = s15f16_roundtrip(value)
return params
# The code below uses IccFromXML and IccApplyNamedCmm
# Requires DemoICCMax command line tools which you can download at:
# https://www.color.org/iccmax/index.xalter (Implementation section)
# Alternatively, build from https://github.com/InternationalColorConsortium/DemoIccMAX
# Expect some difficulties, including some RPATH issues when using CMake and macOS.
# DEMOICCMAX_BINS = "/Users/remi/ColorCode/DemoIccMAX/TestingBins-2.1.10"
DEMOICCMAX_BINS = "/Users/remi/ColorCode/DemoIccMAX/cmake_build/myinstall/bin"
def make_para_icc_profile(icc_profile_path, params):
""" Generate a MTX+TRC ICC profile with the requested ParametricCurve."""
para_type = params.get("t")
para_params = " ".join([str(p) for p in list(params.values())[1:]])
icc_profile = f"""<?xml version="1.0" encoding="UTF-8"?>
<IccProfile>
<Header>
<ProfileVersion>4.4</ProfileVersion>
<ProfileDeviceClass>mntr</ProfileDeviceClass>
<DataColourSpace>RGB</DataColourSpace>
<PCS>XYZ</PCS>
<CreationDateTime>2022-01-01T00:00:00</CreationDateTime>
<RenderingIntent>Perceptual</RenderingIntent>
<PCSIlluminant>
<XYZNumber X="0.96420288" Y="1.00000000" Z="0.82490540"/>
</PCSIlluminant>
<ProfileCreator>OCIO Project</ProfileCreator>
</Header>
<Tags>
<copyrightTag>
<textType>
<TextData><![CDATA[Copyright Contributors to the OpenColorIO Project. Generated by DemoICCMax's IccFromXML.]]></TextData>
</textType>
</copyrightTag>
<mediaWhitePointTag>
<XYZArrayType>
<XYZNumber X="0.95045471" Y="1.00000000" Z="1.08905029"/>
</XYZArrayType>
</mediaWhitePointTag>
<mediaBlackPointTag>
<XYZArrayType>
<XYZNumber X="0.00000000" Y="0.00000000" Z="0.00000000"/>
</XYZArrayType>
</mediaBlackPointTag>
<redColorantTag>
<XYZArrayType>
<XYZNumber X="0.43606567" Y="0.22248840" Z="0.01391602"/>
</XYZArrayType>
</redColorantTag>
<greenColorantTag>
<XYZArrayType>
<XYZNumber X="0.38514709" Y="0.71687317" Z="0.09707642"/>
</XYZArrayType>
</greenColorantTag>
<blueColorantTag>
<XYZArrayType>
<XYZNumber X="0.14306641" Y="0.06060791" Z="0.71409607"/>
</XYZArrayType>
</blueColorantTag>
<redTRCTag>
<parametricCurveType>
<ParametricCurve FunctionType="{para_type}">
{para_params}
</ParametricCurve>
</parametricCurveType>
</redTRCTag>
<greenTRCTag SameAs="redTRCTag"/>
<blueTRCTag SameAs="redTRCTag"/>
</Tags>
</IccProfile>
"""
# Create a new ICC profile with iccFromXML
with tempfile.NamedTemporaryFile(suffix=".xml", mode="w") as xml_path:
xml_path.writelines(icc_profile)
xml_path.flush()
try:
stdout = sp.check_output([
f"{DEMOICCMAX_BINS}/IccFromXML",
xml_path.name, # XML input file
icc_profile_path, # ICC output file
], text=True)
except sp.CalledProcessError as e:
print("Error", e.stdout)
CCS_D50 = colour.CCS_ILLUMINANTS["CIE 1931 2 Degree Standard Observer"]["D50"]
CCS_D65 = colour.CCS_ILLUMINANTS["CIE 1931 2 Degree Standard Observer"]["D65"]
def apply_icc_profile(x, icc_profile, invert=False, intent=1):
""" Apply the input ICC profile to the input ndarray.
The data will be interpreted as RGB device triplets for the forward
direction and XYZ (D50) PCS triplets otherwise. One dimensional input
arrays will be interpreted as grayscale R=G=B or D50 Y depending on the
direction requested.
Rendering intent parameter follow IccApplyNamedCmm enumeration.
"""
FP = 15 # Float precision for text data interchange
in_dims = len(x.shape)
with tempfile.NamedTemporaryFile(mode='w') as f:
# PCS to Device direction.
if invert == True:
# Lab or XYZ are supported as input, Lab seems to introduce a slight
# increase in error relative compared to Python.
# Just change the Data Format header from XYZ to Lab and apply
# the required conversion when choosing to input Lab triplets.
f.write("'XYZ ' ; Data Format\n")
f.write("icEncodeValue ; Encoding\n")
for i in range(x.shape[0]):
if in_dims == 2:
XYZ = x[i, :]
else:
XYZ = colour.xyY_to_XYZ([CCS_D50[0], CCS_D50[1], x[i]])
f.write(f"{(XYZ[0]):.{FP}f} {(XYZ[1]):.{FP}f} {(XYZ[2]):.{FP}f}\n")
# Device to PCS direction.
else:
f.write("'RGB ' ; Data Format\n")
f.write("icEncode8bit ; Encoding\n")
for i in range(x.shape[0]):
if in_dims == 2:
RGB = np.rint(x[i, :] * 255)
else:
RGB = np.tile(np.rint(x[i] * 255), 3)
f.write(f"{RGB[0]} {RGB[1]} {RGB[2]}\n")
f.flush()
# Print input data (debugging only)
# with open(f.name, "r") as f:
# print(f.read())
try:
stdout = sp.check_output([
f"{DEMOICCMAX_BINS}/IccApplyNamedCmm",
f.name, # Input data
f"0:{FP}", # Final data encoding : Precision (icEncodeValue)
"0", # Interpolation (linear)
os.path.abspath(icc_profile), # ICC profile
f"{intent}" # Rendering intent (relative)
], text=True)
# Print output data (debugging only)
# print(stdout)
d = r"\s*(-?\d+\.\d+)"
matchs = re.findall(fr"^{d}{d}{d}\s*;{d}{d}{d}\s*$", stdout, re.MULTILINE)
results = np.zeros((len(matchs), 6))
for idx, match in enumerate(matchs):
results[idx, :] = match[:]
return results[:, :3]
except sp.CalledProcessError as e:
print(e.stdout)
#
# ParametricCurve - OCIO implementation
#
def apply_icc_para_ocio(x, invert=False, **kw):
""" Apply the ParametricCurve using OCIO processing.
Construct an OCIO config with the corresponding transforms on the fly.
"""
# Increasing the LUT size seems to cause loss of precision in the inverse
# direction, could be related to the gamma curve being too flat near 0 and
# hitting float precision errors?
OCIO_LUT_SIZE = 2**10
with tempfile.NamedTemporaryFile(suffix=".spi1d") as lut_path:
config = OCIO.Config.CreateRaw()
config.addColorSpace(OCIO.ColorSpace(name="Linear"))
# Simple ExponentTransform for type 0
if kw.get("t", None) == 0:
g = kw.get("g")
config.addColorSpace(OCIO.ColorSpace(
name="ParaCurve",
toReference=OCIO.ExponentTransform(value=[g, g, g, 1])))
# Generate 1DLUT for type 1-4
else:
lut = colour.LUT1D(size=OCIO_LUT_SIZE)
lut.table = apply_iccpara_python(lut.table, **kw)
colour.write_LUT(lut, lut_path.name)
config.addColorSpace(OCIO.ColorSpace(
name="ParaCurve",
toReference=OCIO.FileTransform(src=lut_path.name)))
# Pixel processing
if invert:
proc = config.getProcessor("Linear", "ParaCurve")
else:
proc = config.getProcessor("ParaCurve", "Linear")
x_rgb = np.tile(x[..., np.newaxis], (1, 3)).astype(np.float32)
proc.getDefaultCPUProcessor().applyRGB(x_rgb)
return x_rgb[..., 1]
#
# ParametricCurve - Python implementation
#
def qb(v, n=8):
""" Quantize the input float ndarray to n bits."""
return np.around(v * (2**n - 1)).astype(np.int32) / (2**n - 1)
def check_paracurve_params(params):
""" ParametricCurve specification does not impose constraints on parameters,
some combinations can produce ill-defined curves with a range of issues.
Implements the checks recommended by the ICC (p.15):
https://www.color.org/whitepapers/ICC_White_Paper35-Use_of_the_parametricCurveType.pdf
"""
t = params.get("t")
g = params.get("g")
a = params.get("a")
b = params.get("b")
c = params.get("c")
d = params.get("d")
e = params.get("e")
f = params.get("f")
###########################################################################
# Expected number of parameters per paraCurve type
###########################################################################
num_params = {
0: 1,
1: 3,
2: 4,
3: 5,
4: 7,
}
if len(params) - 1 != num_params[t]:
raise ValueError("Invalid number of parameter for paraCurve")
###########################################################################
# Monotically non-decreasing (flat segments permitted) check
###########################################################################
# g > 0 force the power law to be monotically non-decreasing
if g <= 0.0:
raise ValueError(f"ParaCurve {t} expect 0 < x <= 1 but got {g}")
# a > 0 force the argument to the power law to be an increasing function of x
if (t != 0) and a <= 0.0:
raise ValueError(f"ParaCurve {t} expect a > 0 but got {a}")
# c >= 0 force the linear segment to be flat or increasing
if (t == 3 or t == 4) and c < 0:
raise ValueError(f"ParaCurve {t} expect c >= 0 but got {c}")
# check for negative discontinuity at the linear segment boundary
if (t == 3) and not qb(c * d) <= qb(np.power(a * d + b, g)):
raise ValueError(f"ParaCurve {t} expect cd <= (ad + b)^g")
# check for negative discontinuity at the linear segment boundary
if (t == 4) and not qb(c * d + f) <= qb(np.power(a * d + b, g) + e):
raise ValueError(f"ParaCurve {t} expect cd + f <= (ad+b)^g + e")
###########################################################################
# No complex / imaginary numbers check
###########################################################################
# negative arguments to the power can produce imaginary numbers
if (t == 3 or t == 4) and (a * d + b) < 0:
raise ValueError(f"ParaCurve {t} invalid parameters (imaginary numbers not allowed)")
###########################################################################
# Boundary warnings
###########################################################################
# Breakpoint is at x = -b/a, assuming a > 0, b should be negative so that
# the breakpoint occurs at positive x values.
if (t == 1 or t == 2) and b >= 0.0:
warn(f"ParaCurve {t} expect b >= 0 but got {b}")
if t == 1 and qb(a + b) != 1:
warn(f"ParaCurve {t} does not reach maximum at (1,1)")
if t == 2 and qb(np.power(a + b, g) + c) != 1:
print((np.power(a + b, g) + c))
warn(f"ParaCurve {t} does not reach maximum at (1,1)")
###########################################################################
# Continuity warnings
###########################################################################
# Note that type 0, 1, 2 are continuous by definition
if (t == 3) and not qb(c * d) == qb(np.power(a * d + b, g)):
warn(f"ParaCurve {t} is not continuous")
if (t == 4) and not qb(c * d + f) == qb(np.power(a * d + b, g) + e):
warn(f"ParaCurve {t} is not continuous")
###########################################################################
# Smoothness warnings
###########################################################################
# Note that type 0 is smooth by definition
if (t == 1 or t == 2) and g <= 1 and -b / a > 0:
warn(f"ParaCurve {t} is not smooth")
if (t == 3 or t == 4) and not qb(c) == qb(a * g * np.power(a * d + b, g - 1)):
print(qb(c), qb(a * g * np.power(a * d + b, g - 1)))
warn(f"ParaCurve {t} is not smooth")
def subsitute_paracurve_params(params):
""" ParametricCurve specification does not impose constraints on parameters,
some combinations can produce ill-defined curves with a range of issues.
Implements the subsitutions recommended by the ICC (p.16):
https://www.color.org/whitepapers/ICC_White_Paper35-Use_of_the_parametricCurveType.pdf
"""
t = params.get("t")
g = params.get("g")
a = params.get("a")
b = params.get("b")
c = params.get("c")
d = params.get("d")
e = params.get("e")
f = params.get("f")
if g <= 0:
g = 1
if t != 0 and a <= 0:
a = 1
if (t == 3 or t == 4) and (a * d + b) < 0:
d = -b / a
if 0 < d < 1:
if t == 3 and (c * d) > np.power(a * d + b, g):
c = (1 / d) * np.power(a * d + b, g)
elif t == 4 and (c * d + f) > (np.power(a * d + b, g) + e):
if f > (np.power(a * d + b, g) + e):
f = np.power(a * d + b, g) + e
c = (1 / d) * (np.power(a * d + b, g) + e - f)
for key, val in zip(params.keys(), [t,g,a,b,c,d,e,f]):
if params[key] != val:
params[key] = val
return params
def apply_iccpara_python(x, invert=False, simulate_s15f16n=False, **kw):
""" Apply the ParametricCurve using Numpy.
Clipping is performed on both input / output side as recommended by the
ICC, even when not strictly needed.
"""
t = kw.get("t")
g = kw.get("g")
a = kw.get("a")
b = kw.get("b")
c = kw.get("c")
d = kw.get("d")
e = kw.get("e")
f = kw.get("f")
# Para 0
if t == 0 and not invert:
x = np.clip(x, 0, 1)
return np.clip(np.power(x, g), 0, 1)
elif t == 0 and invert:
x = np.clip(x, 0, 1)
return np.clip(np.power(x, 1 / g), 0, 1)
# Para 1
elif t == 1 and not invert:
x = np.clip(x, 0, 1)
return np.clip(
np.where(
x >= -b / a,
np.power(a * x + b, g),
np.zeros_like(x)
)
, 0, 1
)
# Flat segment inversion is ambigous, any value in the range [0, -b/a] is
# valid, here we match OCIO 1DLUT behaviour by returning -b/a. Note that
# DemoICCMax's IccApplyNamedCmm returns 0 here.
elif t == 1 and invert:
x = np.clip(x, 0, 1)
return np.clip((np.power(x, 1 / g) - b) / a, 0, 1)
# Para 2
elif t == 2 and not invert:
x = np.clip(x, 0, 1)
return np.clip(
np.where(
x >= -b / a,
np.power(a * x + b, g) + c,
c
)
, 0, 1
)
elif t == 2 and invert:
x = np.clip(x, 0, 1)
return np.clip(
np.where(
x > c,
(np.power(x - c, 1 / g) - b) / a,
c
)
, 0, 1
)
# Para 3
elif t == 3 and not invert:
x = np.clip(x, 0, 1)
return np.clip(
np.where(
x >= d,
np.power(a * x + b, g),
c * x
)
, 0, 1
)
elif t == 3 and invert:
x = np.clip(x, 0, 1)
p = np.power(a * d + b, g)
return np.clip(
np.where(
x >= p,
(np.power(x, 1 / g) - b) / a,
x / c
)
, 0, 1
)
# Para 4
elif t == 4 and not invert:
x = np.clip(x, 0, 1)
return np.clip(
np.where(
x >= d,
np.power(a * x + b, g) + e,
c * x + f
)
, 0, 1
)
elif t == 4 and invert:
x = np.clip(x, 0, 1)
p = np.power(a * d + b, g) + e
return np.clip(
np.where(
x >= p,
(np.power(x - e, 1 / g) - b) / a,
(x - f) / c
)
, 0, 1
)
raise NotImplementedError
#
# Testing
#
def check_python_roundtrip(params):
x = np.linspace(0, 1, 256)
y = apply_iccpara_python(apply_iccpara_python(x, **params), invert=True, **params)
np.testing.assert_allclose(x, y, rtol=1e-15, atol=1e-15)
def check_python_icc_match(params):
with tempfile.NamedTemporaryFile(suffix=".icc", mode="w") as icc_path:
make_para_icc_profile(icc_path.name, params)
# Apply with apply_icc
x = np.linspace(0, 1, 256)
atol = 1e-4
rtol = 1e-4
# Test the PCS to device direction
y_python = apply_iccpara_python(x, **params, invert=True)
y_icc = apply_icc_profile(x, icc_path.name, invert=True)[:, 1]
np.testing.assert_allclose(y_python, y_icc, rtol=rtol, atol=atol)
# Test the Device to PCS direction
# Note that we get XYZ back here because it is our profile PCS
y_python = apply_iccpara_python(x, **params, invert=False)
y_icc = apply_icc_profile(x, icc_path.name, invert=False)[:, 1]
np.testing.assert_allclose(y_python, y_icc, rtol=rtol, atol=atol)
def check_python_ocio_match(params):
x = np.linspace(0, 1, 256)
atol = 1e-4
rtol = 1e-4
y_python = apply_iccpara_python(x, **params, invert=True)
y_ocio = apply_icc_para_ocio(x, **params, invert=True)
np.testing.assert_allclose(y_python, y_ocio, rtol=rtol, atol=atol)
y_python = apply_iccpara_python(x, **params, invert=False)
y_ocio = apply_icc_para_ocio(x, **params, invert=False)
np.testing.assert_allclose(y_python, y_ocio, rtol=rtol, atol=atol)
def plot_curve(title, params):
plt.title(title)
x = np.linspace(0, 1, 2**14)
y = apply_iccpara_python(x, **params)
y_inv = apply_iccpara_python(x, invert=True, **params)
y_ocio = apply_icc_para_ocio(x, **params)
y_ocio_inv = apply_icc_para_ocio(x, invert=True, **params)
plt.grid(True)
plt.plot(x, y, label="Python")
plt.plot(x, y_inv, '--', label="Python Inv")
plt.plot(x, y_ocio, label="OCIO")
plt.plot(x, y_ocio_inv, '--', label="OCIO Inv")
plt.legend()
plt.show()
def ocio_unittest_helper(profile_name, params):
with np.printoptions(edgeitems=30, linewidth=100000):
make_para_icc_profile(profile_name, params)
x = np.array([-1.0, 0.0, 0.18, 0.5, 0.75, 1.0, 2.0])
y = apply_icc_profile(x, profile_name)[:, 1]
y_py = apply_iccpara_python(x, **params)
y_inv = apply_icc_profile(x, profile_name, invert=True)[:, 1]
y_py_inv = apply_iccpara_python(x, invert=True, **params)
y_py_bck = apply_iccpara_python(y_py_inv, **params)
# For type 1 and 2, OCIO match python behviour for flat segment inversion
print(profile_name, "FWD: ", y)
print(profile_name, "INV: ", y_inv)
print(profile_name, "FWD (py):", y_py)
print(profile_name, "INV (py):", y_py_inv)
print(profile_name, "BCK (py):", y_py_bck)
#
# Script entry point
#
if __name__ == "__main__":
# Check para type 0
type0_params = {"t": 0, "g": 2.4}
# plot_curve("Type 0", type0_params)
check_paracurve_params(type0_params)
check_python_roundtrip(type0_params)
check_python_ocio_match(type0_params)
check_python_icc_match(type0_params)
# Check para type 1
type1_params = {"t": 1, "g": 2.4, "a": 1.1, "b": -0.1}
# plot_curve("Type 1", type1_params)
check_paracurve_params(type1_params)
# Cannot roundtrip with flat segment
# check_python_roundtrip(type1_params)
check_python_ocio_match(type1_params)
# Inverting flat segment ambiguity, must be commented
# check_python_icc_match(type1_params)
ocio_unittest_helper("icc-test-pc1.icc", type1_params)
# Check para type 2
type2_params = {"t": 2, "g": 2.4, "a": 1.057049452, "b": -0.1, "c": 0.1}
# plot_curve("Type 2", type2_params)
check_paracurve_params(type2_params)
# Cannot roundtrip with flat segment
# check_python_roundtrip(type2_params)
# Inverting flat segment ambiguity, must be commented
# OCIO 1DLUT inversion returns a value that seem to be not matching either
# our Python nor ICC implementation
# check_python_ocio_match(type2_params)
# Inverting flat segment ambiguity, must be commented
# check_python_icc_match(type2_params)
ocio_unittest_helper("icc-test-pc2.icc", type2_params)
# Check para type 3
# rec709_params = {"t": 3, "g": 1/0.45, "a": 1/1.099, "b": 0.099/1.099, "c": 1/4.5, "d": 0.081}
srgb_params = {"t": 3, "g": 2.4, "a": 1/1.055, "b": 0.055/1.055, "c": 1/12.92, "d": 0.04045}
check_paracurve_params(srgb_params)
check_python_roundtrip(srgb_params)
check_python_ocio_match(srgb_params)
check_python_icc_match(srgb_params)
# plot_curve("Type 3 sRGB", srgb_params)
ocio_unittest_helper("icc-test-pc3.icc", srgb_params)
# Check para type 4
# Naive attempt to get continuous and smooth curve parameters (starting from sRGB parameters)
g = 2.4
b = 0.055/1.055
d = 0.04045
e = 0.1
a = np.power(1 - e, 1 / g) - b
c = a * g * np.power(a * d + b, g - 1)
f = np.power(a * d + b, g) + e - c * d
type4_params = {"t": 4, "g": g, "a": a, "b": b, "c": c, "d": d, "e": e, "f": f}
# plot_curve("Type 4", type4_params)
check_paracurve_params(type4_params)
check_python_roundtrip(type4_params)
check_python_ocio_match(type4_params)
check_python_icc_match(type4_params)
ocio_unittest_helper("icc-test-pc4.icc", type4_params)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment