Haoyue Xiao

---

Introduction

The Prokudin-Gorskii collection, led by Russian photographer Sergei Mikhailovich Prokudin-Gorskii in the early 20th century, is a pioneering project documenting the Russian Empire in color. Using an innovative method commissioned by Tsar Nicholas II, Prokudin-Gorskii captured three black-and-white images of the same scene through red, green, and blue filters. When combined, these images produced vivid color photographs, offering a rare glimpse into life over a century ago. Now preserved by the Library of Congress, the collection includes 1,902 digitized triple-frame glass negatives.

Prokudin-Gorskii believed in the future of color photography, and his method involved creating triple-frame negatives for vibrant projections. The task today is to align these three color channels—red, green, and blue—accurately, recreating sharp and true-to-original color images. Proper alignment avoids blurring and ensures the photographs reflect Prokudin-Gorskii's vision.

---

Methology

The general approach to align the color channels will be firstly find a base channel (eg: the blue channel) and align the other two channels to it. Then, we stack all channels together in the order of r, g, b to create a colored image. Thus, the tricky part is how to align two images effectively and efficiently.

Alignment Metrics

To quantify how well two images align, we need a metric fucntion in the following form:

The SSD Norm

One easy-to-think-of way is to use the \(l_2\) distance of the brightness pixel map between the two images: Here, we add a negative sign to the sum of squares so that we can maximize it instead of minimizing.
This method is convenient to implement and effectively gives a degree of similarity between two monochromatic images. However, it has some draw backs such as it treats every pixel equally, and the \(l_2\) norm is sensitive to outliers, which can be problematic when the two film images have random dark spots (due to technology at that time) or irregular borders.

We can adopt a more robost metric using formula \(\mathrm{E_\text{SRD}}(A, B) = -\sum_i \sum_j \rho\left(A(i, j) - B(i,j)\right)\) where \(\rho(\cdot)\) is some robust norm function such as the \(l_1\) norm.

The NCC Norm

Another more interesting (and more effective in my implementation) method is to use the Normalized Cross-correlation (NCC) as the metric, with the formula:

where \(\bar{A} = \frac{1}{N}\sum_{i,j}A(i,j)\) is the mean of all pixel brightness in \(A\) and \(||\cdot||_F\) is the Frobinus norm of matrices, and the subtraction between \(A\) and \(\bar{A}\) is element-wise.

The NCC metric effectively measures the correlation between the images (in pixel brightness) and gives a degree of how similar two images are. However, both the \(SSD\) and \(NCC\) are applied to pixel brightness, which is technically a risky choice for alignment. Due to the filming conditions and the filters, the optimal alignment probbaly do not have the same pixel brightness across each channel, and there may be a worse alignment that achieves better score than the "human-judged" optimal alignment when using the metric on brightness.

Sobel Edge Detector

To solve this issue, we must analyze the invariant features accorss all color channels that can be used for pattern matching. In my solution, I decide to use edge of objects in the image to match, because when filming the same object, the edges are less likely to change during the short interval of filming.

To apply edge detection on images, we convolve the original image with the sobel operator both horizontally and vertically. The sobel operator is define as

where convolving with \(\mathrm{Sobel_{X}}\) gives detection on vertical edges and \(\mathrm{Sobel_{Y}}\) derects horizontal edges. Example results are:

and when combined, it gives the full edge map of the original image:

Using \(\mathrm{NCC}\) on the edge map, we are able to achieve much more accurate alignment than simple pixel brightness.

Naive Alignment

The naive alignment is to search every possible shift between two images.
This method will in theory give the best alignment according to the selected metric. Given that the images we are trying to align should be generally aligned with slight displacement, it's also desirable to search best shift between a selected range [-displacement, displacement]. In my naive implementation, I chose the displacement to be 15, and that means we search within the grid \([-15, 15]^2\). Here are some example aligned low-resolution images:

Pyramid Search (Hierarchical Motion Estimation)

However, for high resolution images, the pixel displacement can be much larger (pixels become "smaller") and it's very expensive to search all possible shifts even within a selected range. Therefore, we need more clever ways such as constructing a gaussian pyramid.

In my approach, I first compute the pyramid layer to be \(\left\lceil\log_2\left(\frac{\min(\text{H}, \text{W})}{\text{min edge pixels}}\right)\right\rceil\), where min edge pixels is the threshold that the edge of the smallest image in the pyramid should be larger than, eg. 25 pixles.

Then, from the original image, we continuously downsample the image in the previous layer to half of its size, until we reach the final layer.

To search for the best alignment, we move from top to the bottom using this approach:

1. Compute a full search on the lowest-resolution image to get an initial assumption of the position. This should be fast if the layer of the pyramid is high enough. In my implementation, by choosing min edge pixels=25, I get a pyramid of layer 7 and the final resolution can be as low as \(45\times 52\).



2. For the rest images in the pyramid, we start our search from the best alignment in the last layer. Because the resolition of the next image is doubled, if the best shift from last layer is (x, y), then this layer will search starting from (2x, 2y).

3. For layers larger than 1, we only search within the range \([-1, 1]\) around the starting point. If the best displacement from the start point is (a,b), then pass the location (2x + a, 2y + b) to the next layer.

This pyramid approach is much faster than the naive search even for low resolution images.

Before aligning two images, I apply some preprocessing as the following:

• Crop the image to the 90% center to avoid aligning irregular image borders.
• Optionally, we can apply Sobel operator to the images and align the edge of the image, which in practice, achieves much better performance than aligning pixel brightness directly.

Here are some aligned example photos using this approach(both with and without edge detection):

---

With NCC Only (Not Perfectly Aligned)

With NCC + Sobel Edge Detector

---

Fancy Stuff

Fourier-Based Fast NCC Computation

The Fourier transform approach is based on the fact that convolution in the spatial domain corresponds to multiplication in the Fourier domain, namely:

where \(A_N = \frac{A - \overline{A}}{\|A - \overline{A}\|_F}\) is the normalized image.

Fourier Runtime Analysis

Suppose all images have the same dimension \(m \times n\), let's do a runtime analysis on the naive approach, the pyramid approach, and the FFT approach:

1. To search for the optimal shift using the naive approach, we need to search every possible combination in the \(x\)-\(y\) directions, and that's \(\mathcal{O}(mn)\) combinations. Even if we only search a selected range, as the image size grows, this range should also scale proportionally. And at each possible shift, we evaluate the metric which involves all pixels, and that's at least \(\mathcal{O}(mn)\). Therefore, the runtime of the naive approach is \(\mathcal{O}(m^2n^2)\).


2. In the pyramid approach, we always downsample the original image to the lowest resolution above some threshold (layer 1). Therefore, the exhaustive search at this level is a constant. For the rest of the layers, we only search in a \(3 \times 3\) grid, and each evaluation takes \(\mathcal{O}\left(\frac{mn}{4^i}\right)\) time. Therefore, the total runtime is
   $$ \begin{equation*} \begin{aligned} T&(m, n)\\ &= \underbrace{\mathcal{O}(mn) + \mathcal{O}\left(\frac{mn}{4}\right) + \mathcal{O}\left(\frac{mn}{4^2}\right) + \cdots + \mathcal{O}(1)}_{\log_2{\min{(m,n)}} \text{ terms}} \\ &= \mathcal{O}(mn) + C \end{aligned} \end{equation*} $$
   where \(C\) is the exhaustive search constant.


3. In the Fourier-based approach, using the fast Fourier transform, we can compute the Fourier transform of the images in \(\mathcal{O}(mn \log(mn))\). The conjugate multiplication takes \(\mathcal{O}(mn)\), and the inverse Fourier transform takes \(\mathcal{O}(mn \log(mn))\). Therefore, we can get the NCC map of every possible shift in \(\mathcal{O}(mn \log(mn))\), and finding the argmax is \(\mathcal{O}(mn)\). Finally, the \(\mathcal{O}(mn \log(mn))\) term dominates.

Even though asymptotically, the pyramid approach is better than the Fourier approach, but in practice, the pyramid approach has a huge overhead and behaves worse than the FFT. Here's a graph showing the relationship between runtime and input size:

Here's also a table summarizing the runtime and best shift