Project 1: Colorizing the
Prokudin-Gorskii Collection

E Harrison  —  CS180, Fall 2026

CS180 Projects  /  Home  /  Assignment

Overview

In the early 1900s, Sergei Prokudin-Gorskii photographed the Russian Empire in color long before color printing existed. He recorded each scene three times through a red, a green, and a blue filter onto a single glass plate. The digitized negatives — now held by the Library of Congress — stack the three exposures vertically (top to bottom: B, G, R).

The goal of this project is to take a stacked plate, split it into its three color channels, and align them into a single, artifact-free color photograph. The plate below shows the three raw exposures — notice how simply stacking them without alignment produces color fringing around edges.

Part 1: Naive Single-Scale Alignment

The simplest approach fixes the blue channel and searches for the (x, y) translation of the green and red channels that best matches it. I exhaustively try every displacement in a window of [-15, 15] pixels, score each candidate with an image-matching metric, and keep the best one.

I scored candidates with two metrics. The L2 (Euclidean) distance measures the raw pixel difference between the two channels, while normalized cross-correlation (NCC) takes the dot product of the two channels after mean-subtracting and normalizing them to unit length. NCC is more robust to overall brightness differences between channels, so I used it as my default metric.

import numpy as np
import cv2 as cv

# read in the image as grayscale and convert to float in [0,1]
im = cv.imread(imname, cv.IMREAD_GRAYSCALE)
im = im.astype(np.float32) / 255.0

# compute the height of each part (just 1/3 of total) and separate color channels
height = int(np.floor(im.shape[0] / 3.0))
b = im[:height]
g = im[height: 2*height]
r = im[2*height: 3*height]

def l2_dist(im1, im2):
    im1_vec = im1.flatten()
    im2_vec = im2.flatten()
    return -np.sqrt(np.sum((im1_vec - im2_vec) ** 2)) # negative since we want to minimize this score

def ncc(im1, im2):
    def normalize_and_scale(im):
        im_vec = im.flatten()
        im_shifted = im_vec - np.mean(im_vec)
        return im_shifted / np.linalg.norm(im_shifted)
    im1_norm, im2_norm = normalize_and_scale(im1), normalize_and_scale(im2)
    return np.dot(im1_norm, im2_norm)

def shift_image(im, shift):
    shifted_im = np.roll(im, shift[0], axis=0)
    shifted_im = np.roll(shifted_im, shift[1], axis=1)
    return shifted_im

def best_fit(im1, im2, matching_metric, shift_range):
    best_shift = (None, None)
    best_score = -np.inf
    for i in range(-shift_range, shift_range + 1):
        for j in range(-shift_range, shift_range + 1):
            shifted_im1 = shift_image(im1, (i, j))
            score = matching_metric(shifted_im1, im2)
            if score > best_score:
                best_shift = (i, j)
                best_score = score
    return best_shift

def align(im1, im2):
    shift = best_fit(im1, im2, matching_metric=ncc)
    return shift

# align G and R to the B channel, then stack into a color image
g_shift = align(g, b)
r_shift = align(r, b)
ag = shift_image(g, g_shift)
ar = shift_image(r, r_shift)
im_out = np.dstack([ar, ag, b])

This exhaustive search is cheap enough to run directly on the small, low-resolution JPEGs. My first attempt scored each candidate over the entire channel — and the results were not quite right.

The culprit is the ragged border of each glass-plate scan. The black, white, and discolored edges do not agree across the three channels, so when they are included in the matching metric they dominate the score and pull the search toward a bad displacement. monastery shows this most dramatically — the heavy red/cyan ghosting means the channels never locked onto the real content.

My fix is to score only the interior of each channel. Before evaluating the metric, I crop a centered window out of both channels, so the noisy borders are ignored and the metric reflects only the actual scene.

def crop_image(im, crop_size):
    left = (im.shape[0] - crop_size) // 2
    right = (im.shape[0] + crop_size) // 2
    bottom = (im.shape[1] - crop_size) // 2
    top = (im.shape[1] + crop_size) // 2
    return im[left:right, bottom:top]

def best_fit(im1, im2, matching_metric, shift_range, crop=True):
    best_shift = (None, None)
    best_score = -np.inf
    for i in range(-shift_range, shift_range + 1):
        for j in range(-shift_range, shift_range + 1):
            shifted = shift_image(im1, (i, j))
            a, b = shifted, im2
            if crop:
                a = crop_image(shifted, 200)
                b = crop_image(im2, 200)
            score = matching_metric(a, b)
            if score > best_score:
                best_shift = (i, j)
                best_score = score
    return best_shift

With the borders cropped out of the scoring, all three plates align cleanly. The displacement vectors listed below are the (row, col) shifts applied to the green and red channels.

The animation below walks through the entire single-scale process on cathedral: the plate is split into thirds and assigned to B, G, R, then the green channel is scanned over the fixed blue channel. The left panel is a red/cyan overlay (grey where the channels agree, colored fringing where they don't), and the right panel is the NCC score at every displacement, filling in as the search proceeds. The star marks the best match. The red channel is aligned the same way, and the animation finishes by sliding the three color filters together into the final image.

Part 2: Multi-Scale Alignment with an Image Pyramid

The full-resolution glass-plate scans are .tif files thousands of pixels on a side, and their true displacements can be well over 100 pixels. A brute-force search over such a large window would be hopelessly slow — the cost grows with the square of the search range and the number of pixels scored.

The fix is a coarse-to-fine image pyramid. I recursively downscale the image by a factor of two until it is small, estimate the alignment on the tiny coarse image (where a large displacement becomes a small one), then propagate that estimate down the pyramid — doubling it at each level and refining it with a small [-10, 10] local search. As in Part 1, I score only the cropped interior of each channel so the noisy borders don't corrupt the match.

import skimage as sk

def best_fit(im1, im2, matching_metric, shift_range, initial_shift, crop):
    best_shift = (None, None)
    best_score = -np.inf
    x_shift, y_shift = initial_shift
    for i in range(-shift_range + x_shift, shift_range + x_shift + 1):
        for j in range(-shift_range + y_shift, shift_range + y_shift + 1):
            shifted_im1 = shift_image(im1, (i, j))
            if crop:
                shifted_im1 = crop_image(shifted_im1, 200)
                im2 = crop_image(im2, 200)
            score = matching_metric(shifted_im1, im2)
            if score > best_score:
                best_shift = (i, j)
                best_score = score
    return best_shift

def align(im1, im2, threshold=300):
    # Image pyramid
    if im1.shape[0] >= threshold:
        small_im1 = sk.transform.rescale(im1, 0.5)
        small_im2 = sk.transform.rescale(im2, 0.5)
        i, j = align(small_im1, small_im2)
        i *= 2
        j *= 2
    else:
        i, j = (0, 0)

    shift = best_fit(
        im1, 
        im2, 
        matching_metric=ncc, 
        shift_range=10,  
        initial_shift=(i, j),
        crop=True
    )
    return shift

The animation below shows the coarse-to-fine search on harvesters. The image is first downscaled into a pyramid; the alignment is estimated on the tiny coarsest level (where a 120-pixel shift becomes a handful of pixels), and that estimate is then doubled and refined with a small local search at each finer level until the full resolution is reached.

With the pyramid, each high-resolution plate aligns in well under a minute. Below are my results on the provided example images, with the computed (row, col) displacement vectors listed beneath each one.

A Failure Case: Emir of Bukhara

One image does not align correctly with my basic pipeline. The green channel lines up, but the red channel lands at (24, -4) — less than half of the green channel's (48, 23), when it should be roughly double. The result is the heavy red ghosting seen below.

The problem is that my metrics compare raw pixel intensities. The Emir's bright silk robe reflects very differently through the red, green, and blue filters, so the same region has very different brightness in each channel. NCC and L2 therefore find a false optimum for the red channel — it matches the overall brightness pattern rather than the true structure. The fix would be to align on edges or gradients (which are consistent across channels) instead of raw intensities.

My Own Selections

Finally, I chose three additional plates from the Library of Congress Prokudin-Gorskii collection and ran them through the same pyramid pipeline.

Website template from Jon Barron.