Otsu et. al (2018) is the second spectral upsampling method I'm adding to Colour. It's a data-driven method, which not only produces reflectances matching the given color exactly, but also that are natural and resemble real, measured reflectances. This is in contrast with, for example, Jakob and Hanika (2019), which I worked on in June and uses a mathematical model.

At the moment the code is split in two parts. The feature branch in my fork implements the method using the dataset published with the article. Creating new datasets is done by this program and will be covered in depth in a later article. Work is mostly done, I just need to integrate the program into Colour and complete this series of blog posts.

Spectral reconstruction in this method is very simple and straightforward, thus also remarkably fast to compute. Let $\mathbf{s}$ be the reconstructed reflectance vector: $$\mathbf s = w_1\mathbf b_1 + w_2\mathbf b_2 + w_3\mathbf b_3 + \mathbf b_\mu$$ where $\mathbf b_i$ are three basis vectors and $\mathbf b_\mu$ is the dataset mean. The original code used spectra discretized at 36 points (at $\Delta \lambda = 10~\mathrm{nm}$ intervals), making the vectors 36-dimensional. Since the four vectors are known (precomputed ahead of time), reconstruction is simply a matter of computing the weights $w_i$. First rewrite the equation in a matrix form, treating vectors as column vectors: $$ \mathbf s = \begin{bmatrix}\mathbf b_1 & \mathbf b_2 & \mathbf b_3\end{bmatrix} \mathbf w + \mathbf b_\mu $$ In a previous article I've already shown how to compute CIE 1931 XYZ tristimulus values given a reflectance. Written in a matrix form, with $\mathbf c$ being a vector of the tristimulus values: $$ \mathbf c = k\Delta\lambda\mathbf I \begin{bmatrix}\mathbf{\bar x}^T
\mathbf{\bar y}^T
\mathbf{\bar z}^T\end{bmatrix} \mathbf s $$ The $\Delta\lambda$ term appears because of discretization. Combining the two equations gives: $$ \begin{aligned} \mathbf c &= k\Delta\lambda\mathbf I

\begin{bmatrix}\mathbf{\bar x}^T \\ \mathbf{\bar y}^T \\ \mathbf{\bar z}^T\end{bmatrix}
\left(\begin{bmatrix}\mathbf b_1 & \mathbf b_2 & \mathbf b_3\end{bmatrix}
\mathbf w + \mathbf b_\mu\right) \\
&= \underbrace{k\Delta\lambda\mathbf I
\begin{bmatrix}\mathbf{\bar x}^T \\ \mathbf{\bar y}^T \\ \mathbf{\bar z}^T\end{bmatrix}
\begin{bmatrix}\mathbf b_1 & \mathbf b_2 & \mathbf b_3\end{bmatrix}}_{M}
\mathbf w + \mathbf c_\mu

\end{aligned} $$ A quick note about the dimensions of $M$. The matrix of cmfs is three 36-dimensional row vectors stacked vertically, a 3×36 matrix, while the basis matrix is three column vectors stacked horizontally, a 36×3 matrix. Multiplying the two yields a 3×3 matrix, which can be inverted (degenerate cases aside).

Solving this equation is very simple. Move $\mathbf c_\mu$ to the left side and left-multiply both sides by the inverse of $M$. $$M^{-1}\left(\mathbf c - \mathbf c_\mu\right) = \mathbf{w}$$ Knowing the weights is enough to do the reconstruction. An interesting thing to note is that $M^{-1}$ and $\mathbf c_\mu$ could be pre-computed (for a given illuminant and cmfs), making the reconstruction process potentially extremely fast. Also, an RGB-to-XYZ matrix could be pre-multiplied into $M$ to speed up upsampling of RGB data.

This method works on discrete spectra. The range and resolution are fixed and determined by the pre-computed dataset. On the other hand, this method could theoretically be used also used for emissive spectra, so it can be quite versatile.

While reconstruction is very simple, coming up with good basis vectors is not. It's a complex and time consuming process that I'll cover in the next (at least) two parts.