Scanning Negatives the Hard Way


Lately I've been getting back into film photography. Some years ago I took a fair number of medium format pictures using a Holga pinhole camera. I never digitized any of my negatives, mostly due to the lack of a decent scanner. I still lack a proper scanner, having only a sub-par Brother model meant for documents. However, I may be able to use a few mathematical tricks to make it work far past its intended potential.

scan of film negative
A poor-quality film scan.

As seen in my first scan, the main problem is I can't control the exposure time on my scanner and negatives end up much too dark. While some information can be recovered by correcting the contrast and gamma of each image, this is overwhelmed by the noise produced by the scanner. Moreover, I was lazy so the film was obscured by dust, hairs, and fingerprints. Fortunately, the random noise can be removed by averaging multiple scans. However, I also wanted to try removing various sources of systematic noise (i.e., prints, hairs, film imperfections). Ultimately, I took 125 scans of the negative in four different positions on the scanner bed, each with their own imperfections. After reading each image into Matlab, I corrected the contrast, gamma, and registered each image to align differences in position and rotation.

img = zeros(2692,5738,125);
for i = 1:125
    % read in scans and correct gamma and contrast limits
    img(:,:,i) = 1-mat2gray(double(imread(['Image_',num2str(i),'.bmp'])));
    img(:,:,i) = imadjust(img(:,:,i),[0.5, 1],[],10);

% register images, allowing for rotation and translation
mkdir Registered % directory for registered images
regImg = zeros(2692,5738,125);
[optimizer, metric] = imregconfig('monomodal'); % assume similar brightness and contrast
optimizer.MaximumIterations = 1000; % increase iterations
regImg(:,:,27)=img(:,:,27); % reference image

idx = [1:26, 28:125];
for i=1:length(idx)
    regImg(:,:,i) = imregister(img(:,:,i), img(:,:,27), 'rigid', optimizer, metric);

meanImg = mean(regImg,3);
figure; % display mean image

gamma and contrast corrected film scans
Noisy, gamma and contrast corrected scans.

The end result is pretty rough. However, we can remove the noise produced by the scanner by taking the average of each pixel across all scans.

averaged scan
The average of all 125 negative scans.

However, two sources of noise still remain. The various hairs and fingerprints present in the different scans have to be removed. Second, regular vertical lines are present throughout the negative. While these would not be visible in a proper scan, the increased image contrast has made them more present. While the fingerprints will be difficult to address, the lines are regular and can be removed my moving the image into the frequency domain.

axis image;

amplitude spectra
The amplitude spectra of the average image.

As seen above, a strong horizontal line is present throughout the frequency spectra of our image, representing the regular pattern of vertical lines. By removing these frequencies, we should also remove the lines without affecting the rest of our image. So, we will need a function to filter out frequencies ±2° from vertical. As the lines are closely spaced, we don't need to filter out low frequencies, only those above 200 cycles/image.

function thefilter = ofilter(oLow,oHigh,dim);
    % a bandpass filter of size dim that passes orientations between oLow and oHigh
    % horizontal and vertical are 0 and 90 degrees
    t1 = ((0:dim-1)-floor(dim/2))*(2/(dim));
    [t1,t2] = meshgrid(t1,t1);
    clear t2
    if oHigh>=oLow
        d=find((a=oHigh) & t1>0); % find out-of-band frequencies
        d=find((a=oHigh) & t1>0); % find out-of-band frequencies
    if ~isempty(d)

function thefilter = bpfilter(lowf,highf,dim);
    % a bandpass filter of size dim that passes frequencies between lowf and highf
    for yk=1:dim
        d=sqrt(dx.*dx + dy*dy);
        val1=find((d > lowf)&(d < highf));

function thefilter=filter2d(lowf,highf,orient1,orient2,dim);
    % combines ofilter and bpfilter to pass frequencies between lowf and highf
    % and orientations between orient1 and orient2

Given the function filter2d, we can apply it to each image to remove all frequencies above 200 cycles/image that are within 2° of vertical.

thefilter = filter2d(200,5738,88,92,5738);
thefilter = 1-thefilter(1524:4215,:);

for i=1:125
    amp = amp.*thefilter;
    ft = amp.*exp(sqrt(-1)*phase);

filtered image
Frequency filtered image.

Finally! We have removed the lines and are one step closer to having a clean picture. All that remains is to remove the various hairs and fingerprints present in each image. Fortunately, the position of these flaws varied between scans. I wiped down the negative, removing some prints and causing others to be deposited. Since presence or absence of these flaws are independent of the image, we can use independent component analysis (ICA) to isolate them and subtract them from our average image. The FastICA package is one of the best implementations of ICA.

filtImg=filtImg(59:2643,50:5678,:); % to avoid exceeding RAM capacity
[rows, columns, scans] = size(filtImg);
pcaArray = reshape(filtImg, rows * columns, scans)'; % reshape array for ICA
mu = mean(pcaArray,1); % the average image, ICA automatically subtracts this
[icasig, A] = fastica(pcaArray, 'lastEig', 8, 'numOfIC', 8); % take first 8 components
baseCase = mean(A(76:125,:),1); % scores for fingerprint-free images

% positive and negative parts of the first three components
ICAp1 = icasig(1,:); ICAp1(ICAp1 < 0) = 0;
ICAm1 = -icasig(1,:); ICAm1(ICAm1 < 0) = 0;
ICAp2 = icasig(2,:); ICAp2(ICAp2 < 0) = 0;
ICAm2 = -icasig(2,:); ICAm2(ICAm2 < 0) = 0;
ICAp3 = icasig(3,:); ICAp3(ICAp3 < 0) = 0;
ICAm3 = -icasig(3,:); ICAm3(ICAm3 < 0) = 0;

independent components
The first four independent components.

As seen above, ICA has isolated all the fingerprints and hairs that obscured our image within the first three independent components. Since all scans had defects, the absence of a flaw in one image implies the presence of another elsewhere. Therefore, each component captures both the presence (positive values & lighter pixels) and absence (negative values & darker pixels) of various features. As we want to make an image that contains none of these flaws, we separate the positve and negative values of each component.

All that remains is to find the ideal combination of component weights that removes all blemishes from the average image (mu). My approach was to generate random weights for two sets of components, add these to the average image, and indicate which looked better. It is a bit hacky, but a few iterations resulted in an image that looked half decent.

WeightArr = [rand(1,2).* 0.026 - 0.013; ...
    rand(1,2).* 0.016 - 0.008; ...
    rand(1,2).* 0.0608 - 0.0304];
for i = 1:100
    tmpArr = [rand(1,2).* 0.026 - 0.013; ...
        rand(1,2).* 0.016 - 0.008; ...
        rand(1,2).* 0.0608 - 0.0304];
    figure('units','normalized','outerposition',[0 0 1 1]);
    imshow(reshape((mu+ICAp1.*WeightArr(1,1)+ICAp2.*WeightArr(2,1)+ICAp3 ...
    .*WeightArr(3,1)+ICAm1.*WeightArr(1,2)+ICAm2.*WeightArr(2,2)+ICAm3 ...
    imshow(reshape((mu+ICAp1.*tmpArr(1,1)+ICAp2.*tmpArr(2,1)+ICAp3 ...
    .*tmpArr(3,1)+ICAm1.*tmpArr(1,2)+ICAm2.*tmpArr(2,2)+ICAm3 ...
    w = waitforbuttonpress;
    if w
        WeightArr = tmpArr;

imwrite(uint16(mat2gray(reshape((mu+ICAp1.*WeightArr(1,1)+ICAp2.*WeightArr(2,1) ... 
+ICAp3.*WeightArr(3,1)+ICAm1.*WeightArr(1,2)+ICAm2.*WeightArr(2,2)+ICAm3 ...
.*WeightArr(3,2) )', rows, columns) ).*65535),'FinalImg.tif');

the final image
The final image, after filtering, ICA, and averaging.

And here it is, the final image. While it still could use a lot of work, the end result is a far cry from what my scanner initially produced. However, the main lesson is that I should probably buy a better scanner.