Initial commit

4 files changed, 233 insertions, 0 deletions

diff --git a/site/content/gsoc/otsu2018_1.typ b/site/content/gsoc/otsu2018_1.typ
new file mode 100644
index 0000000..f86db9b
--- /dev/null
+++ b/site/content/gsoc/otsu2018_1.typ
@@ -0,0 +1,58 @@
+#import "/common.typ": post, link-self
+
+#show: post.with(
+    title: "Otsu et al. (2018), part 1",
+    date: datetime(year: 2020, month: 7, day: 23),
+)
+
+Otsu et al. (2018) @Otsu2018 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) @Jakob2019, which I worked on in June and uses a mathematical model.
+
+At the moment the code is split in two parts. The #link("https://github.com/enneract/colour/tree/feature/otsu2018")[feature branch] in my fork implements the method using the dataset published with the article. Creating new datasets is done by #link("https://github.com/enneract/otsu2018")[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.
+
+= Reconstruction
+
+Spectral reconstruction in this method is very simple and straightforward, thus also remarkably fast to compute. Let $bold(s)$ be the reconstructed reflectance vector:
+$
+  bold(s) = w_1 bold(b)_1 + w_2 bold(b)_2 + w_3 bold(b)_3 + bold(b)_mu
+$
+where $bold(b)_i$ are three *basis vectors* and $bold(b)_mu$ is the dataset mean. The original code used spectra discretized at 36 points (at $Delta lambda = 10 "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:
+$
+  bold(s) =
+    mat(bold(b)_1, bold(b)_2, bold(b)_3)
+    bold(w) + bold(b)_mu
+$
+In a #link-self("gsoc/differentiating_code")[previous article] I've already shown how to compute CIE 1931 XYZ tristimulus values given a reflectance. Written in a matrix form, with $bold(c)$ being a vector of the tristimulus values:
+$
+  bold(c) = k Delta lambda bold(I)
+    mat(overline(bold(x))^T; overline(bold(y))^T; overline(bold(z))^T)
+    bold(s)
+$
+The $Delta lambda$ term appears because of discretization. Combining the two equations gives:
+$
+  bold(c) &= k Delta lambda bold(I)
+    mat(overline(bold(x))^T; overline(bold(y))^T; overline(bold(z))^T)
+    (mat(bold(b)_1, bold(b)_2, bold(b)_3)
+    bold(w) + bold(b)_mu) \
+  &= underbrace(
+    k Delta lambda bold(I)
+    mat(overline(bold(x))^T; overline(bold(y))^T; overline(bold(z))^T)
+    mat(bold(b)_1, bold(b)_2, bold(b)_3),
+    M,
+  )
+    bold(w) + bold(c)_mu
+$
+A quick note about the dimensions of $M$. The matrix of cmfs is three 36-dimensional row vectors stacked vertically, a 3x36 matrix, while the basis matrix is three column vectors stacked horizontally, a 36x3 matrix. Multiplying the two yields a 3x3 matrix, which can be inverted (degenerate cases aside).
+
+Solving this equation is very simple. Move $bold(c)_mu$ to the left side and left-multiply both sides by the inverse of $M$.
+$
+  M^(-1) (bold(c) - bold(c)_mu) = bold(w)
+$
+Knowing the weights is enough to do the reconstruction. An interesting thing to note is that $M^(-1)$ and $bold(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.
+
+= Notes
+
+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 also be 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.
+
+#bibliography("/refs.bib")
diff --git a/site/content/gsoc/otsu2018_2.typ b/site/content/gsoc/otsu2018_2.typ
new file mode 100644
index 0000000..4ae754b
--- /dev/null
+++ b/site/content/gsoc/otsu2018_2.typ
@@ -0,0 +1,79 @@
+#import "/common.typ": post, link-self, static, img
+
+#show: post.with(
+    title: "Otsu et al. (2018), part 2",
+    date: datetime(year: 2020, month: 7, day: 24),
+)
+
+The #link-self("gsoc/otsu2018_1")[previous post] mentions three basis vectors and a "dataset mean." Here I explain how to find a good basis and what "good" means in this context.
+
+The idea is to find a basis so that *any* spectral reflectivity can be expressed in the form:
+$
+  bold(s) = sum_i w_i bold(b)_i + bold(b)_mu
+$
+where $i$ runs over all dimensions of the space (36 using Otsu et al.'s spectral data @Otsu2018) and $bold(b)_mu$ is a constant term you can ignore for now. This is a very common problem and there are many general solutions. By far the most common one is the discrete Fourier transform, which turns arbitrary signals into weighted sums of sines and cosines.
+
+Unfortunately, recalling that the goal is to upsample a color triplet (e.g. XYZ tristimulus values or RGB data), there are only three degrees of freedom to work with and so at most three weights can be determined. Reformulating the problem,
+$
+  bold(s) = w_1 bold(b)_1 + w_2 bold(b)_2 + w_3 bold(b)_3 + bold(epsilon) + bold(b)_mu
+$
+where all but three terms are lumped together into an error term $bold(epsilon)$. I think it's clear that three sines or cosines (and a constant term) added together won't accurately reconstruct most spectra. There are far better purely mathematical models, such as Jakob and Hanika (2019) @Jakob2019.
+
+= Principal component analysis
+
+Otsu et al. came up with a *data-driven* solution. Given enough measured data, you can use statistics and linear algebra to create a small basis that carries as much information from the original dataset as possible. The authors used *principal component analysis* (PCA) and a set of $N = 1,279$ measured reflectances.
+
+Spectral reflectances are treated as random vectors, with each component a random variable. This method requires that their sample means are zero, so the measured data has to be adjusted first. Let
+$
+  bold(b)_mu = 1 / N sum_i bold(s)_i
+$
+which I called the "dataset mean" in the previous article, then simply subtract it from measured spectra:
+$
+  bold(s)'_i = bold(s) - bold(b)_mu
+$
+From these the so-called *data matrix*, in this case an N x 36 matrix, can be constructed:
+$
+  X = mat(bold(s)'_1, bold(s)'_2, dots.h, bold(s)_N)^T
+$
+Now the idea is to use this to derive a basis that maximizes its variance or "significance." I'll skip the proofs, as I'm not particularly well-versed in statistics and it would take me too much effort. The result is given by the eigensystem of the covariance matrix $Sigma$ of $X$. The eigenvector of $Sigma$ with the largest eigenvalue is the vector maximizing variance. The eigenvector with the second largest eigenvalue is perpendicular and has the second largest variance, and so on.
+
+Let's take a look at the code. Assume that the measured data, `sds`, is an N x 36 array.
+```
+basis_mean = np.mean(sds, axis=0)
+data_matrix = sds - basis_mean
+```
+A covariance matrix can be easily computed with `np.cov`, but a trick can be used. Because the column-wise means are zero, $X^T X$ is going to be equal to it up to a constant. The constant and any normalization can be ignored, as only eigenvectors are sought after.
+```
+covariance_matrix = np.dot(data_matrix.T, data_matrix)
+```
+The last step is the eigendecomposition. NumPy comes with a few algorithms, so it's not a problem. Because `cov` is Hermitian (it's a real, symmetric matrix), one can use:
+```
+eigenvalues, eigenvectors = np.linalg.eigh(covariance_matrix)
+```
+The bonus feature of `eigh` (aside from being more efficient) is that the eigenvalues are sorted in ascending order. The result is simply the last three columns of `eigenvectors`: the three basis vectors with the greatest variance.
+
+Because zero-mean spectra were used to create the basis, it best describes primed spectra:
+$
+  bold(s)' = w_1 bold(b)_1 + w_2 bold(b)_2 + w_3 bold(b)_3
+$
+Adding $bold(b)_mu$ to both sides,
+$
+  bold(s) = w_1 bold(b)_1 + w_2 bold(b)_2 + w_3 bold(b)_3 + bold(b)_mu
+$
+recovers the equation from the beginning of this post and explains why the fourth basis vector has kept appearing. Notice that the error term is missing: the reconstructed spectrum is slightly different. I'll discuss this in part 3.
+
+_Thanks to Sıddık Açil for #link("https://dev.to/akaame/implementing-simple-pca-using-numpy-3k0a")[his article] on performing PCA in NumPy and Hao Xie for letting me know about it._
+
+= Visualizing
+
+Below you can see plots of a few eigenvectors derived from Otsu et al.\'s dataset. Since the vectors are meaningful only up to a constant, the Y scale can be chosen arbitrarily (and differently for each vector). The usual way to visualize them is to make a "ladder" plot, with the vectors stacked vertically according to their number. Eigenvalues are written on the right. The script is available #link("https://github.com/enneract/otsu2018/blob/master/pca_decomposition.py")[here].
+
+The results look #link("https://en.wikipedia.org/wiki/Quantum_harmonic_oscillator")[oddly familiar]. Notice that every next eigenvector (starting from the top) has exactly one more zero crossing than the previous one. My guess is that this has to be the case for them to be perpendicular and linearly independent. The eigenvalues are interesting too. They decrease very rapidly, with the largest few being much, much larger than the others. This explains why just three basis functions are enough for a good reconstruction.
+
+#img(static("gsoc/otsu2018_eigenstates_1.png"))
+
+As expected, the eigenvectors corresponding to small eigenvalues are noisy and random-looking.
+
+#img(static("gsoc/otsu2018_eigenstates_2.png"))
+
+#bibliography("/refs.bib")
diff --git a/site/content/gsoc/otsu2018_3.typ b/site/content/gsoc/otsu2018_3.typ
new file mode 100644
index 0000000..aa1666e
--- /dev/null
+++ b/site/content/gsoc/otsu2018_3.typ
@@ -0,0 +1,79 @@
+#import "/common.typ": post, link-self, static, img
+
+#show: post.with(
+    title: "Otsu et al. (2018), part 3",
+    date: datetime(year: 2020, month: 7, day: 26),
+)
+
+In the #link-self("gsoc/otsu2018_2")[previous post] I explained how PCA can be used to create a convenient basis for representing real-world reflectivities. The procedure uses a set of measured spectra to compute an orthonormal basis, sorted by variance or "statistical significance." All but the three most significant vectors are then truncated. I called the difference between the actual spectrum and the reconstructed error term $bold(epsilon)$.
+
+Naturally the smaller the error, the better. Because $bold(epsilon)$ is a vector, the authors chose the Euclidean norm (2-norm or L2) to represent the (scalar) error:
+$
+  epsilon = sqrt(sum_i bold(epsilon)_i^2)
+$
+To get an idea of how well the method works in the general case, the authors would sum up reconstruction errors of *all* spectra in the dataset. I call this the total error $E$:
+$
+  E = sum_i epsilon |_(bold(s)_i)
+$
+where $epsilon |_(bold(s)_i)$ denotes the reconstruction error of $bold(s)_i$, the $i$-th measured spectrum.
+
+Using just one basis for everything isn't ideal. The picture below shows reconstructions of a few colorchecker spectra. While by definition XYZ tristimulus values (and perceived colors) match perfectly, the spectra aren't that close.
+
+#figure(
+  img(static("gsoc/otsu_clustering1.png")),
+  caption: [Reconstructions of a few colorchecker spectra.],
+)
+
+= Clustering
+
+The authors @Otsu2018 decided to subdivide the CIE 1931 xy space using a binary space partitioning tree. Every node that isn't a leaf divides its subspace in two with a horizontal or vertical line. Leaves, called clusters in the article, contain spectra that lie in their subspace. PCA is done for every cluster separately, creating a set of bases.
+
+#figure(
+  img(static("gsoc/chromaticity_diagram.png")),
+  caption: [CIE 1931 xy chromaticity diagram],
+)
+
+A quick note about CIE 1931 xy chromaticity coordinates. They're derived from CIE 1931 XYZ tristimulus values:
+$
+  x &= X / (X + Y + Z) \
+  y &= Y / (X + Y + Z)
+$
+They're often used where only the chromaticity matters and not the lightness $Y$. When plotted, they give rise to the familiar "color triangle" or "horseshoe" diagram, pictured on the right.
+
+The algorithm for creating an optimized clustering is described as follows.
+- Search every possible new partitioning in the tree. For every color of every leaf, try partitions with a horizontal and a vertical line going through the color, assigning the color itself to the "lesser" part (to the left or below the line).
+- Discard partitions that would create very small nodes. More on this later.
+- Compute how each partition would affect the total error of the entire tree (the sum of $E$ of every leaf). Discard those that would increase it.
+- Realize the partition that decreases the error the most.
+- Repeat the entire process.
+The algorithm terminates when no more partitions are possible or there are enough clusters. Authors chose 8 clusters, because of diminishing returns and dataset size concerns. You can look at #link("http://lightmetrica.org/h-otsu/download/rgb2spec.zip")[the original implementation] or #link("https://github.com/enneract/otsu2018/tree/d218ef0097174897634477bcc180b6268ba6c337")[my code]. The article also has many more plots and comparisons than I post here, so it's definitely worth looking at too.
+
+Below are the results of running my code. I used the entire dataset of 1,279 spectra and optimized for 8 clusters. This takes about 3 hours on my dated i7 3770.
+
+#figure(
+  img(static("gsoc/otsu_clustering2b.png")),
+  caption: [A visualization of the clustering. Each cluster has its own symbol for the colors that belong to it. (I actually used the cross symbol twice, my bad.)],
+)
+
+#figure(
+  img(static("gsoc/otsu_clustering2a.png")),
+  caption: [Reconstructions of the same colorchecker spectra as before, this time using clustering. The discrepancies are notably smaller and some reconstructed spectra look exactly like the measurements.],
+)
+
+= Minimum node size
+
+I mentioned that nodes cannot be allowed to be split too finely. First of all, there have to be at least three spectra in a node for PCA to even work. Fewer will result in a singular $M$ and the algorithm simply won't be able to continue. In practice, having few spectra is no good either. Below is a result of PCA decomposition done with just three spectra.
+
+#figure(
+  img(static("gsoc/otsu2018_eigenstates_3_colors.png")),
+  caption: [PCA decomposition of the first three spectra in Otsu et al.\'s dataset. An explanation of what is plotted here is in the #link-self("gsoc/otsu2018_2")[previous post].],
+)
+
+Only *two* eigenvectors are meaningful and the rest are basically noise. Trying to use this for reconstruction will create completely random-looking spectra, usually going way outside the [0, 1] range. Clearly a lot more data needs to be used to get useful results. In my program I (arbitrarily) set the minimum cluster size to the average cluster size divided by two and it seems to work pretty well.
+$
+  "min. cluster size" = floor(1 / 2 (("num. of spectra") / ("num. of clusters")))
+$
+
+But can reconstructed spectra exceed the allowed [0, 1] range even if the basis is good? Unfortunately yes, and this will be the topic of the next (and the last) post in this series.
+
+#bibliography("/refs.bib")
diff --git a/site/content/gsoc/otsu2018_4.typ b/site/content/gsoc/otsu2018_4.typ
new file mode 100644
index 0000000..a5d01db
--- /dev/null
+++ b/site/content/gsoc/otsu2018_4.typ
@@ -0,0 +1,17 @@
+#import "/common.typ": post, link-self, static, img
+
+#show: post.with(
+    title: "Otsu et al. (2018), part 4",
+    date: datetime(year: 2020, month: 8, day: 10),
+)
+
+In the #link-self("gsoc/otsu2018_3")[previous part] I mentioned that reflectances reconstructed using this method can exceed the [0, 1] range and violate conservation of energy. Indeed, I could simply ask it for an imaginary color and get a weird, unphysical spectrum.
+
+What's interesting is that this doesn't happen in practice, when upsampling real data. Out of almost 1300 reflectances I tested, I found only one bad reconstruction. Using the original clustering, the result for `magenta` from `BabelColor Average` colorchecker reflects 100,6% of light at 730 nm. Clipping this value to 100% creates a tiny color difference, orders of magnitude smaller than the #link("https://en.wikipedia.org/wiki/Just-noticeable_difference")[just-noticeable difference].
+
+#figure(
+  img(static("gsoc/otsu2018_error.png")),
+  caption: [You can see the small difference between the clipped and unclipped reflectances in the upper right corner of the figure. (Click to enlarge.)],
+)
+
+Many methods have this property of being able to create perfect metamers, so looking at ΔE's isn't very interesting. In a future post I'll take a look at MSE (mean square error), which quantifies differences between spectra themselves, and compare a few methods.