https://www.itu.int/dms_pubrec/itu-r/rec/bt/R-REC-BT.2124-0-201901-I!!PDF-E.pdf

delta_e_itp.py

"""See https://www.itu.int/dms_pubrec/itu-r/rec/bt/R-REC-BT.2124-0-201901-I!!PDF-E.pdf."""

import numpy as np
import numpy.typing as npt

_bgr2lms_matrix = (
    np.array(
        [
            [262, 462, 3688],
            [2146, 2951, 309],
            [1688, 683, 99],
        ]
    )
    / 4096
)


def _bgr2lms(bgr: npt.NDArray) -> npt.NDArray:
    """Convert BGR to LMS."""
    return np.dot(bgr, _bgr2lms_matrix)


_m1 = 2610 / 16384
_m2 = 2523 / 4096 * 128
_c1 = 3424 / 4096
_c2 = 2413 / 4096 * 32
_c3 = 2392 / 4096 * 32


def _pq_inverse(f: npt.NDArray) -> npt.NDArray:
    """The perceptual quantizer inverse electro-optical transfer function."""
    y = f / 10000
    return ((_c1 + _c2 * y**_m1) / (1 + _c3 * y**_m1)) ** _m2


_lms2itp_matrix = (
    np.array(
        [
            [2048, 6610 / 2, 17933],
            [2048, -13613 / 2, -17390],
            [0, 7003 / 2, -543],
        ]
    )
    / 4096
)


def _lms2itp(lms: npt.NDArray) -> npt.NDArray:
    """Convert LMS to ITP."""
    return np.dot(lms, _lms2itp_matrix)


def _delta_e_itp(itp1: npt.NDArray, itp2: npt.NDArray) -> npt.NDArray:
    """Calculate the delta E ITP between two ITPs."""
    return 720 * np.sqrt(np.sum((itp1 - itp2) ** 2, axis=-1))


def delta_e_itp(bgr1: npt.NDArray, bgr2: npt.NDArray) -> npt.NDArray:
    """Calculate the delta E ITP between two BGRs."""
    return _delta_e_itp(
        _lms2itp(_pq_inverse(_bgr2lms(bgr1))),
        _lms2itp(_pq_inverse(_bgr2lms(bgr2))),
    )