Sunday, February 28, 2016

Lunar Orbiter photograph restoration, 1































Scan of an original "8x10" print, from the same Lunar orbiter 1 source as frame 1107. The photography was acquired with the 80 mm wide angle medium resolution (MR) lens, originally on 70 mm film, between August 18 to 29, 1966, and readout occurred through September 14, 1966.

Spacecraft position:
     altitude: 53.71 km
     latitude: 0.28
     longitude: 12.71
Illumination:
     sun azimuth: 88.5
     incident angle: 69.2
     phase angle: 68.8

Several spacecraft, as part of the Lunar Orbiter program, went to the moon, took detailed photographs and sent them back to earth. The recorded images have a primary artifact of mostly regular striping in one direction. It is not difficult to suppress most of the striping and other artifacts, resulting in images much closer to what we would see.

Detail, bottom-right corner of the above full scan.



The basic idea behind removing these vertical stripe artifacts is to adjust each pixel such that the average luminance of any vertical segment is nearly equal. To do this, the average luminance along each stripe, centered on each pixel, is calculated and the difference with respect to an arbitrary global average is subtracted from the pixel.

Why design restoration algorithms? Mostly because I really like looking at these images, wondering about what's there. Just below the resolution of these images is an unknown landscape, terra incognita. At various scales and image contrast, the moon looks like an interesting place.





Detail at full resolution, from near center-right of the above detail.
The moon is superficially stark and simple, with a limited range of terrain and feature types. On the other hand, every crater is a new place to explore, and walking between boulders and across slopes is easy to imagine.




There are thousands of Lunar Orbiter images available: see the Lunar Orbiter Photo Gallery. This library has much of the striping artifacts removed, but not aggressively: some banding is left intact to preserve subtle detail and large scale variation in apparent albedo.

Recently better photographs from a similar perspective are available from the Lunar Reconnaissance Orbiter. The lighting will generally be slightly different, but the lunar features haven't changed much in the last 50 years. A digital elevation model, at a resolution of about 400 m/pixel, is also available.

Today some of these striping removal and other corrections would have been applied directly, before distribution. At the time they were received, this wasn't possible or feasible, and to prevent loss of detail the artifacts were left intact.

I'm interested only as a test case; how well can reasonable corrections be made, without losing key information content? I want them to look somewhat like we'd see them visually, from orbit. The photography was carefully choreographed to have a low sun angle that highlights subtle surface features, ideal for both a photographic record and for visual navigation.

At the time there were some good guesses of what the surface looks like, but most of the details were new. Many of the features still aren't well understood. Any restoration can be compared to the newer, better images, as well as some newer, better reconstructions from a subset of the original recorded data (LOIRP).

This variety of simple restoration is now common, and the techniques are widely applicable, to imagery and other signals.


Lunar Orbiter background

The Lunar Orbiter program was a series of five unmanned lunar orbiter missions launched by the United States from 1966 through 1967. Intended to help select Apollo landing sites by mapping the Moon's surface,[1] they provided the first photographs from lunar orbit.

Lunar Orbiter Photographic Atlas of the Moon
[photo subsystem diagrams]


Lunar Orbiter Image Recovery Project (LOIRP)

Spaceflight Revolution, chapter on Lunar Orbiter Project history

The Lunar Orbiter spacecraft took photos, scanned them, transmitted them back to earth, and then were printed on paper.

















Lunar Orbiters shot their images on 70mm SO-243 photographic film. The film was developed aboard the spacecraft, scanned and sent back to Earth as analog data. One data stream demodulated the image and sent it to a high resolution kinescope (TV screen) where the image on the screen was photographed as a photographic positive. Those photographs were subsequently re-photographed to create a negative, then stitched together and photographed again to create photographic prints used for the requirements of the program.


Primary destriping algorithm


Most image information of interest is contained in variations at high spatial frequencies, but high frequency spatial differences between columns are partly artifact. A form of high pass filtering, or unsharp masking, can be used to reduce this artifact: subtract low spatial frequency components in the vertical direction, leaving high spatial frequency variations.

This could be done in the spatial frequency domain, and often is. Fourier analysis (DFT) gives an equivalent representation in terms of spatial frequency component amplitudes. The vertical low frequency components can be removed, by not including those components in a reconstruction.

Averaging, in the spatial domain, is similar to finding low spatial frequencies, with the length scale of averaging corresponding to the range of spatial frequencies. High spatial frequency components average to zero, so subtracting averages suppresses low spatial frequency features while preserving high spatial frequency features.

At each pixel, the mean luminance in a tall window in the direction of the striping, is found.

The vertical mean at each pixel is expected to be nearly independent of luminance values to either side. This mean change across the image over short ranges -- the average brightness -- is expected to be relatively constant.



Luminance profiles across two closely spaced vertical profiles. The top profile (red, x2) is on average brighter than the bottom profile (blue, x1).

















Mean difference subtracted from x2, such that the mean luminance's are matched.

This correction for the mean at every pixel is applied to the full image, such that every luminance profile mean is the same. The range is stretched to take advantage of most of the available dynamic range of the 8-bit image.
































Spatial albedo differences across large scales are also suppressed. Broad mountain ranges, large regions in shadow and some plateaus and depressions have albedo differences at large scales, and they are artificially suppressed (little or no contrast) at these spatial scales. If the real variations could be distinguished from artifact, they could be added back.


Detail, bottom-right corner of the above full scan.


Detail at full resolution, from near center-right of the above detail.






























There are remaining artifacts:
  • vertical bands between 'orbit strips', which don't contain enough high spatial frequencies for reasonable reconstruction with this method.
  • variations in deep shadows, where there is little or no dynamic range available
  • dust, high spatial frequency bits and pieces
  • checkerboard texture, moire of scan grid and original sample spacing


Abbreviated program listing, Matlab

(Includes only the core algorithm, operating on a full image.)

% supress_vertical_stripes_v4( strImFilename )
%
% Subtract a mean difference of luminance values of vertical column segments
%  centered on each pixel, leaving high spatial frequency components intact.
 
function supress_vertical_stripes_v4_short( strImFilename )

    % Read image.
    nimIn  = imread( strImFilename ); % 8-bit single luminance channel image
    szimIn = size( nimIn );
 
    % Half-width of band window for mean calculation.
    hwBand = 300;

    % Find a mean of each vertical line segment centered on each pixel:
 
    % For each column:
    for ix = 1 : szimIn(2)

        ix % column index progress indicator
 
        % For each row:   
        for iy = 1 : szimIn(1)

            % Find mean luminance of a vertically alligned set of pixel values,
            % centered on the pixel at (ix,iy):
 
            yRange =   max( iy-hwBand + 1,         1 ) ...
                     : min( iy+hwBand,     szimIn(1) );

            vColSamples = double( nimIn( yRange, ix ) );

            % Find mean of only values nearest to the global mean, excluding
            % high and low values that may be saturated or otherwise less
            % linearly related to the striping artifact.
 
            lNearMean = vColSamples < 235 & vColSamples > 37;
            imMeanColSegment( iy, ix ) = mean( vColSamples( lNearMean ) );
        end
    end
 
 
    % Subtract mean and apply linear luminance transform:
 
    % Subtract difference of means, relative to a global mean, that puts the
    % black point near zero.
    imHorizontal_cor = double( nimIn ) - ( imMeanColSegment - 120 );

    % Stretch entire range, taking brightest pixels to near the white point.
    nimHorizontalHighPass = uint8( 1.40*imHorizontal_cor );

    % Write image.
    imwrite( nimHorizontalHighPass, [ 'destriped_' strImFilename ] );
    end    % end: function supress_vertical_stripes_v4



Full program listing, Matlab.

(Includes bits and pieces for debugging, testing and parameter estimation, and generation of illustrative graphics.)

%
% supress_vertical_stripes_v4( strImFilename )
%
% Subtract a mean difference of luminance values of vertical column segments
%  centered on each pixel, leaving high spatial frequency components intact.
%
% USAGE:
%           supress_vertical_stripes_v4( 'LunarOrbitor_crp_test.png' );
%           supress_vertical_stripes_v4( 'LunarOrbitor_crp.png' );       % full cropped scan
%
% ARGUMENTS:
%
%   strImFilename:  8-bit image full filename.
%                   Note: Must include single filename extension '.'.
%
%
% RETURN VALUES: (none)
%
% HARDCODED: (see below 'function' statement)
%
%       hwBand: % Half-width of column band window for mean calculation.
%
% CALLS: (none)
%
%
% Mark Dow, created       July  6, 2013 (horizontal_stripes_supress_v1,
%                                           derived from horizontal_stripes_slope_v1)
% Mark Dow, modified      July  7, 2013 (v2, simple mean in range +-.5 std)
% Mark Dow, modified  February 15, 2016 (v2, simplify and clarify)
% Mark Dow, modified  February 17, 2016 (v3, outlier range, example graphics)
% Mark Dow, modified  February 19, 2016 (supress_vertical_stripes_v4,
%                                           derived from horizontal_stripes_supress_v3)
% Mark Dow, modified  February 22, 2016 (v4, example profile graphics)
%
 
 
function supress_vertical_stripes_v4( strImFilename )

    % Read image.
    nimIn  = imread( strImFilename ); % 8-bit single luminance channel image
    szimIn = size( nimIn );

    % Show original full image:
    figure; imagesc( nimIn )
    colormap(gray)

    %%%%%%%%%%%%%%%%%%%%%%%%
    % hardcoded:
 
    bWrite = true;

    % Note: These may be overwritten for processing the full image.
    xCropMin =  711;  % limits of image region for 'escarpment' example
    xCropMax = 1160;  % in 'LunarOrbitor_crp_test.png'.
    yCropMin =  301;
    yCropMax =  950;
 %     xCropMin = 1;             % full image
%     xCropMax = szimIn(2);
%     yCropMin = 1;
%     yCropMax = szimIn(1);
 
    % Half-width of column band window for mean calculation.
    hwBand = 300;

    % Example pixel coordinates.
    xEx1 = 282;
    yEx1 = 325;
    xEx2 = 286;
    yEx2 = 325;

    % end: hardcoded
    %%%%%%%%%%%%%%%%%%%%%%%%
 
    iExt = strfind( strImFilename, '.' );

    strCrop = '';

    % If cropped:
    if  xCropMin ~= 1 || xCropMax ~= szimIn(2) ...
      || yCropMin ~= 1 || yCropMax ~= szimIn(1)

        % Show crop region image:
        figure; imagesc( nimIn( yCropMin : yCropMax, xCropMin : xCropMax ) )
        colormap(gray)

        strCrop = [ '_' num2str(xCropMin) '_' num2str(xCropMax) ...
                    '_' num2str(yCropMin) '_' num2str(yCropMax) ];
    end
 
 
%     % Example profile graphics:
%    
%     % Note: Somewhat specific 'LunarOrbitor_crp_test.png', with the crop boundaries:
%     %     xCropMin =  711;  % limits of image region for 'escarpment' example
%     %     xCropMax = 1160;  % in 'LunarOrbitor_crp_test.png'.
%     %     yCropMin =  401;
%     %     yCropMax =  850;
%    
%     nimInEx(1:yCropMax-yCropMin+1,1:xCropMax-xCropMin+1,1) = 0;
%     nimInEx(:,:,1) = nimIn( yCropMin : yCropMax, xCropMin : xCropMax );
%     nimInEx(:,:,2) = nimIn( yCropMin : yCropMax, xCropMin : xCropMax );
%     nimInEx(:,:,3) = nimIn( yCropMin : yCropMax, xCropMin : xCropMax );
%    
%     % Column of luminance values centered on example pixel, rotated to horizontal.
%     vEx1Col = double( nimInEx( yEx1-hwBand + 1 : yEx1+hwBand, xEx1 ) );
%     figure; imagesc( vEx1Col' );
%     colormap(gray)
%    
%     vEx2Col = double( nimInEx( yEx2-hwBand + 1 : yEx2+hwBand, xEx2 ) );
%     figure; imagesc( vEx2Col' );
%     colormap(gray)
%    
%    
%     meanEx1Col = mean( vEx1Col )
%     stdEx1Col  = std(  vEx1Col )
%            
%     meanEx2Col = mean( vEx2Col )
%     stdEx2Col  = std(  vEx2Col )
%    
%     lNearMean = vEx1Col < 235 & vEx1Col >  37;
%     meanEx1Col = mean( vEx1Col( lNearMean ) )
%     lNearMean = vEx2Col < 235 & vEx2Col >  37;
%     meanEx2Col = mean( vEx2Col( lNearMean ) )
%    
%     % Raw column value plots.
%     figure;
%     plot( vEx1Col ); hold on
%     plot( vEx2Col, 'Color', 'r' );
%     plot( ( vEx2Col - vEx1Col ), 'Color', 'm' );
%     line( [1 2*hwBand+1], [0 0], 'Color', [.5 .5 .5], 'LineStyle', ':' )
%     set( gca, 'FontSize', 14 )
%     ylim( [-30 210 ] )
%     ylabel( 'luminance' )
%     xlabel( 'y' )
%     legend( 'x1', 'x2', 'difference' )
%    
%     % After mean shift and feathered reshifted low values.
%    
%     vEx2Col_cor = vEx2Col + (meanEx1Col-meanEx2Col);
%    
%     % Feather reshifted dark saturated region.
%     % Linear shift range below 40 such that dark point at ~30 sees no shift.
%     iLow = find( vEx2Col_cor < 30 );
%     vEx2Col_cor(iLow) = 30;
%    
%     figure;
%     plot( vEx1Col ); hold on
%     plot( vEx2Col + (meanEx1Col-meanEx2Col), 'Color', 'r', 'LineStyle', ':' );
%     plot( vEx2Col_cor, 'Color', 'r' );
%     plot( ( vEx2Col_cor - vEx1Col ), 'Color', 'm' );
%     % line( [1 2*hwBand+1], [(meanEx1Col-meanEx2Col) (meanEx1Col-meanEx2Col)], 'Color', [.5 .5 .5], 'LineStyle', ':' )
%     line( [1 2*hwBand+1], [0 0], 'Color', [.5 .5 .5], 'LineStyle', ':' )
%    
%     set( gca, 'FontSize', 14 )
%     ylim( [-30 210 ] )
%     ylabel( 'luminance' )
%     xlabel( 'y' )
%     legend( 'x1', 'x2, mean shifted', 'x2, corrected', 'difference' )
%    
%     % Draw example pixel profile tracks.
%     nimInEx = double(nimInEx)/255;
%     nimInEx( yEx1-hwBand + 1 : yEx1+hwBand, xEx1, 1 ) = 0;
%     nimInEx( yEx1-hwBand + 1 : yEx1+hwBand, xEx1, 2 ) = 0;
%     nimInEx( yEx1-hwBand + 1 : yEx1+hwBand, xEx1, 3 ) = 1;
%     nimInEx( yEx1, xEx1-2:xEx1+2, 1 ) = 0;
%     nimInEx( yEx1, xEx1-2:xEx1+2, 2 ) = 0;
%     nimInEx( yEx1, xEx1-2:xEx1+2, 3 ) = 1;
%    
%     nimInEx( yEx2-hwBand + 1 : yEx2+hwBand, xEx2, 1 ) = 1;
%     nimInEx( yEx2-hwBand + 1 : yEx2+hwBand, xEx2, 2 ) = 0;
%     nimInEx( yEx2-hwBand + 1 : yEx2+hwBand, xEx2, 3 ) = 0;
%     nimInEx( yEx2, xEx2-2:xEx2+2, 1 ) = 1;
%     nimInEx( yEx2, xEx2-2:xEx2+2, 2 ) = 0;
%     nimInEx( yEx2, xEx2-2:xEx2+2, 3 ) = 0;
%     figure; imagesc( nimInEx );
%    
%     % imwrite( nimInEx, [ strImFilename(1:iExt-1) '_crop_expix' strImFilename(iExt:end) ] );
 
 
    % Find a mean of each vertical line segment centered on each pixel:
 
    % For each column:
    for ix = xCropMin : xCropMax

        ix % column index progress indicator
 
        % For each row:   
        for iy = yCropMin : yCropMax

            % Find mean luminance of a vertically alligned set of pixel values,
            % centered on the pixel at (ix,iy):
 
            yRange =   max( iy-hwBand + 1,         1 ) ...
                     : min( iy+hwBand,     szimIn(1) );

            vColSamples = double( nimIn( yRange, ix ) );

            % Find mean of only values nearest to the global mean, excluding
            % high and low values that may be saturated or otherwise less
            % linearly related to the striping artifact.
 
            lNearMean = vColSamples < 235 & vColSamples >  37;
            imMeanColSegment( iy-yCropMin+1, ix-xCropMin+1 ) ...
                      = mean( vColSamples( lNearMean ) );
        end
    end
 
 
%     max( imMeanColSegment(:) )    % highest average luminance
%     min( imMeanColSegment(:) )    % lowest average luminance
 
    % Subtract mean and apply linear luminance transform:
 
    % Subtract difference of means, relative to a mean that puts the black point near zero.
    imHorizontal_cor =   double( nimIn( yCropMin : yCropMax, xCropMin : xCropMax ) ) ...
                       - ( imMeanColSegment - 120 );%     figure; hist( double( imMeanColSegment(:) ), 200 )
 
    nimHorizontalHighPass = uint8( 1.40*imHorizontal_cor );
 %     % Luminance histogram of original.
%     figure; hist( double( nimIn(:) ), 0:255 )
%     xlim( [-.5 255.5] )
%    
%     % Transformed luminance histogram.
%     figure; hist( double( imHorizontal_cor(:) ), 0:255 )
%     xlim( [-.5 255.5] )
%     figure; hist( double( nimHorizontalHighPass(:) ), 0:255 )
%     xlim( [-.5 255.5] )
 
    % Show cleaned up image.
    figure; imagesc( nimHorizontalHighPass )
    colormap(gray)

    figure; imagesc( imMeanColSegment )
    colormap(gray)
 %     figure; plot( sum(imMeanColSegment,2) )
 
    if bWrite

        % Write cropped, cleaned and artifact (vertical low frequency) images:
 
        imwrite( nimIn( yCropMin : yCropMax, xCropMin : xCropMax ), ...
                [ strImFilename(1:iExt-1) '_crop'  strCrop strImFilename(iExt:end) ] );
        imwrite( nimHorizontalHighPass, ...
                [ strImFilename(1:iExt-1) '_clean' strCrop strImFilename(iExt:end) ] );
        imwrite( uint8( imMeanColSegment ), ...
                [ strImFilename(1:iExt-1) '_VLF'  strCrop strImFilename(iExt:end) ] );
    end
 
end    % end: function supress_vertical_stripes_v4
Posted by at

No comments:

Post a Comment

Subscribe to: Post Comments (Atom)