Contents

SyntheticTerrain is a script created to create synthetic terrain

Allan et al.'s work https://doi.org/10.1109/AERO.2019.8741780 to create a s terrain generator for the lunar surface to create high resolution images from the currently available low resolution images from http://imbrium.mit.edu/

This script upscales the input image and outputs a synthetic higher resolution version of the input image:

1. separates the high and low frequencies of the input image 2. performs fractal synthesis by the diamond-square method on the high frequency portion of the image and uses the Catmull-Rom spline interpolation method to generate new data 3. performs bilinear interpolation on the low frequency component of the image 4. recombines the high and low frequencies of the input image and returns the final image

key input variables the user can change:

c_img - the input image (this should be a dem but the script was mostly tested on optical images. may need different frequency splitting method for low contrast dems avaialable.) the input image needs to be a square with even # of pixels. I recommend starting with a small 500x500 image.

scale - determines the upscaling of the input image. this needs to be a power of 2. i.e. scale = 2, 4, 8, 16 etc.

freq - determines the frequency split between high and low frequencies. make use of the figure showing the frequency split to determine what the frequency should be for the input image. This needs to be between [0,1].

interp - sets the interpolation method. by default this is 'spline' for Catmull-Rom Spline and 'linear', or linear interpolation, is only used on edge cases. However this can be changed to 'linear' for the entire image for comparison. This can only be 'spline' or 'linear'.

key outputs:

hmap - final heightmap. this is the upscaled input image. size will be scale*input image size

method - matrix with indices indicating if a pixel was interpolated with the square or diamond method and the iteration number. 1:iter_max are squares for each iteration. iter_max+1:iter_max+iter_max are diamonds for each iteration. useful for debugging or for understanding final image artifacts.

Online resources:

Diamond-square method https://stevelosh.com/blog/2016/06/diamond-square/

catmull rom spline https://en.wikipedia.org/wiki/Centripetal_Catmull%E2%80%93Rom_spline

%--------------------------------------------------------------------------

%--------------------------------------------------------------------------
%   Copyright (c) 2021 Princeton Satellite Systems, Inc.
%   All rights reserved.
%--------------------------------------------------------------------------
%   Since version 2021.1
%--------------------------------------------------------------------------

% Lunar image
file  = 'LunarSurfaceImage.png';
dem   = imread(file);
demt  = dem/max(dem(:));
c_img	= demt(1:512,1:512);

Fourier transform to separate high and low freq data

Taken from Maggie Kautz's homework for OPTI 512L from UArizona fourier transforms to split high and low frequencies by freq input variable

% Fourier transform original image
F = fftshift(fft2(fftshift(c_img)));

% Make a meshgrid the same size as input image
xi = linspace(-1,1,length(c_img));
eta = xi;
[XI, ETA] = meshgrid(xi,eta);

% Change this variable to change frequency split
freq = 0.01; % bounds [0,1]

% Use meshgrid to greate high and low pass filters
LPF = (abs(XI) < freq) & (abs(ETA) < freq); %low pass filter only lets in xi < abs(freq) frequencies
HPF = 1 - LPF; % high pass filter blocks xi < abs(freq) frequencies

% Use high and low pass filters on fourier transform of the original image
Img_LPF = F.*LPF;
Img_HPF = F.*HPF;

% Inverse fourier transform to return data we can use
img_lpf = fftshift(ifft2(fftshift(Img_LPF)));
img_hpf = fftshift(ifft2(fftshift(Img_HPF)));

% Take real portion of inverse fft for visualization purposes
lpf_data = real(img_lpf);
hpf_data = real(img_hpf);

% Make real portion of fft into rgb images for visualizations at the end of
% the code
lpf_img = imRGB(lpf_data,255);
hpf_img = imRGB(hpf_data,1);

% Plot frequency split to see if freq variable needs to be changed
NewFig('Frequency Split');
subplot(1,2,1)
imagesc(xi,eta, lpf_data)
colormap(gray)
axis square
title('Low Pass Filter')

subplot(1,2,2)
imagesc(xi,eta, hpf_data)
colormap(gray)
axis square
title('High Pass Filter')

sgtitle(sprintf('Frequency %0.2f',freq))

High freq: square-diamond interpolation

% perform square-diamond interpolation only on high frequency portion
image = hpf_data;

% get size of high frequency portion
[m,n] = size(image);

% set upscaling factor (needs to be power of 2: 2,4,8,16, etc.)
scale = 2;

% set interpolation method (either 'linear' or 'spline')
interp = 'spline';

% initilize heightmap and diamond-square method output
hmap = zeros(m*scale,n*scale);
method = zeros(m*scale,n*scale);

% add existing data to heightmap
original_x = zeros(m*n,1);
original_y = zeros(m*n,1);
count = 1;
for i = 1:m
    for j = 1:n
        % position original data determined by scale*x+1,scale*y+1
        % starts at 1,1 when x=y=0
        hmap(scale*(i-1)+1,scale*(j-1)+1) = image(i,j);

        % save position of original data
        original_x(count) = scale*(i-1)+1;
        original_y(count) = scale*(j-1)+1;

        count = count + 1;
    end
end

% rename original data to square input variables
s_in_x = original_x;
s_in_y = original_y;

iteration = 1; % increases for each iteration of square+diamond
iter_max = log(scale)/log(2); % max # of iteration determined by scale
for i = 1:iter_max

[hmap, so_x, so_y, method] = square(hmap, interp, iteration, scale, m, n, s_in_x, s_in_y, method);
[hmap, s_in_x, s_in_y, method] = diamond(hmap, interp, iteration, scale, m, n, so_x, so_y, method);

iteration = iteration + 1;
end

Low freq: bilinear interpolation

% use matlab function to bilinear interpolate low frequency data
bi = imresize(lpf_data,scale,'bilinear');

Rejoin frequencies

recombine processed high and low frequencies

combo = hmap + bi;

% normalize
combo = combo/max(combo,[],'all');

% make rgb for visualization
combo_rgb = imRGB(combo,1);

Plotting

% Original terrain image
NewFig('Original Image');
imshow(c_img)
title('Original Lunar Terrain Image')

% subplots to see before and after for low and high frequencies
NewFig('Low Frequency');
subplot(1,2,1)
imshow(lpf_img)
title('Original Low Frequency')
subplot(1,2,2)
imshow(imRGB(bi,1))
title('Low Frequency Bilinear Interpolation')

NewFig('High Frequency');
subplot(1,2,1)
imshow(hpf_img)
title('Original High Frequency')
subplot(1,2,2)
imshow(imRGB(hmap,1))
title('High Frequency Catmull Rom Spline Interpolation')

% final output image
NewFig('Combined');
imshow(combo_rgb)
title(sprintf('Combined Low and High Frequency Simulated High Resolution Image. Scale %d Split Frequency %0.2f',scale, freq))

Figui

Square interpoloation

function [hmap, output_locations_x, output_locations_y, method] = square(hmap, interp, iteration, scale, im_m, im_n, input_location_x, input_location_y, method)
% function to perform square portion of square-diamond interpolation

fprintf('Performing square iteration %d out of %d\n',iteration,log(scale)/log(2))

% get size of input data
m = length(input_location_x);
n = length(input_location_y);

% initilize arrays to store square locations
square_locations_x = zeros(m,1);
square_locations_y = zeros(n,1);

count = 1;
for i=1:m % for size of input data

    % input data locations
    x = input_location_x(i);
    y = input_location_y(i);

    % if x or y are zero, print and break
    if (x == 0) | (y == 0)
        fprintf('Input data is zero. Miscalculated # of outputs of previous function\n breaking at %d of %d\n', i, m)
        break
    end

    % for each input pixel, find distance to next input pixel in x & y
    distance = scale/(2^(iteration-1));
    x_add = distance;
    y_add = distance;

    interp_2 = interp; % for all data use original interpolation method

    % check if pixel is on edge of image
    x_edge = (x+x_add)>=scale*im_m;
    y_edge = (y+y_add)>=scale*im_n;

    input_data=[]; % initialize input data
    % if input pixel on edge, only take existing data and do linear instead of spline
    % use linear on edge cases bc spline needs 4 data points
    if x_edge & ~y_edge
        input_data = [hmap(x,y),hmap(x,y+y_add)];
        interp_2 = 'linear';
    elseif y_edge & ~x_edge
        input_data = [hmap(x,y),hmap(x+x_add,y)];
        interp_2 = 'linear';
    elseif x_edge & y_edge
        input_data = [hmap(x,y)];
        interp_2 = 'linear';
    else
        % if input pixel not on edge, take surrounding four corner pixels
        try
            input_data = [hmap(x,y),hmap(x,y+y_add),hmap(x+x_add,y),hmap(x+x_add,y+y_add)];
        catch
            fprintf('caught error. usually is an out of bounds error. check input_data')
        end
    end

    % use distance to find placement of new square pixel
    sx = x + distance/2;
    sy = y + distance/2;

    % store square pixel locations
    square_locations_x(count) = sx;
    square_locations_y(count) = sy;

    % depending on if pixel is on edge or not, use linear or spline to
    % interpolate new square pixel value and store in hmap
    switch interp_2
        case 'spline'
            % take spline average
            hmap(sx,sy) = CatmullRomSpline(input_data(:));
        case 'linear'
            % take linear average
            hmap(sx,sy) = mean(input_data,'all');

    end
    count = count + 1;

    % store square method
    method(sx,sy) = iteration;

end

% combine input locations and square locations to pass to diamond function
output_locations_x = [input_location_x; square_locations_x];
output_locations_y = [input_location_y; square_locations_y];
end

Performing square iteration 1 out of 1

Diamond interpolation

function [hmap, output_locations_x, output_locations_y, method] = diamond(hmap, interp, iteration, scale, im_m, im_n, input_location_x, input_location_y, method);
% function to perform diamond portion of square-diamond interpolation

fprintf('Performing diamond iteration %d out of %d\n',iteration,log(scale)/log(2))

% get size of input data
m = length(input_location_x);
n = length(input_location_y);

% initilize arrays to store square locations
diamond_locations_x = zeros(2*m - im_m*iteration - im_n*iteration,1);
diamond_locations_y = zeros(2*n - im_m*iteration - im_n*iteration,1);

count = 1;
for i=1:m % for size of input data

    % input data locations
    x = input_location_x(i);
    y = input_location_y(i);

    % if x or y are zero, print and break
    if (x == 0) || (y == 0)
        fprintf('input data is zero. miscalculated # of outputs of previous function\n breaking at %d of %d\n', i, m)
        break
    end

    % distance to next input pixel in x & y
    distance = scale/(2^(iteration-1));

    % if input location is on the bottom or right edge, don't create a new
    % diamond past the edge of the frame
    start = 1;
    stop = 2;
    if (x + distance/2) >= scale*im_m
        % only create second diamond
        start = 2;
        stop = 2;
    end
    if (y + distance/2) >= scale*im_n
        % only create first diamond
        start = 1;
        stop = 1;
    end
    if ((x + distance/2) >= scale*im_m) & ((y + distance/2) >= scale*im_n)
        % create no diamonds and continue to next input location
        fprintf('Found corner. Continue happening at i = %d out of %d.\n',i,m)
        continue
    end

    % create two diamonds for each input location
    for k=start:stop
        if k == 1
            % first diamond location
            xd = x + distance/2;
            yd = y;
        end
        if k == 2
            % second diamond location
            xd = x;
            yd = y + distance/2;
        end

        % store new diamond location for output later
        diamond_locations_x(count) = xd;
        diamond_locations_y(count) = yd;

        % this is the x,y to add to diamond location to find our data to
        % average
        x_add = distance/2;
        y_add = distance/2;

        interp_2 = interp; % for all data use original interpolation method

        data = []; % initialize data variable
        % find four surrounding pixels for diamond average calculation
        % make sure the x,y we are adding is not past the image frame
        x_edge = (xd + x_add) > scale*im_m;
        y_edge = (yd + y_add) > scale*im_n;

        % edge effects for x
        if xd == 1
            % take only pixel below
            down = hmap(xd + x_add,yd);
            data = [data, down];
            interp_2 = 'linear';
        end
        if x_edge
            % take only pixel above
            up = hmap(xd - x_add,yd);
            data = [data, up];
            interp_2 = 'linear';
        end
        if ~x_edge & xd~=1 % if not on edge
            % take both pixels above and below
            up = hmap(xd - x_add,yd);
            down = hmap(xd + x_add,yd);
            data = [data, up, down];
        end

        %edge effects for y
        if yd == 1
            % take only pixel to the right
            right = hmap(xd,yd + y_add);
            data = [data, right];
            interp_2 = 'linear';
        end
        if y_edge
            % take only pixel to the left
            left = hmap(xd,yd - y_add);
            data = [data, left];
            interp_2 = 'linear';
        end
        if yd~=1 && ~y_edge % if not on edge
            % take both pixels to the right and left
            left = hmap(xd,yd - y_add);
            right = hmap(xd,yd + y_add);
            data = [data, left, right];
        end

        % depending on if pixel is on edge or not, use linear or spline to
        % interpolate new square pixel value and store in hmap
        switch interp_2
            case 'linear' % linear interpolation
                hmap(xd,yd) = mean(data,'all');
            case 'spline' % Catmull Rom spline
                hmap(xd,yd) = CatmullRomSpline(data);
        end
        count = count + 1;

        % store diamond method
        method(xd,yd) = iteration + log(scale)/log(2);
    end
end

% combine input locations and square locations to pass to square function
output_locations_x = [input_location_x; diamond_locations_x];
output_locations_y = [input_location_y; diamond_locations_y];
end

Performing diamond iteration 1 out of 1
Found corner. Continue happening at i = 262144 out of 524288.
Found corner. Continue happening at i = 524288 out of 524288.

Generate an RGB image

function to make input an rgb image. can convert uint 8 to double by using scale = 255.

function image_out = imRGB(image_in,scale)

image_out(:,:,1) = image_in/scale;
image_out(:,:,2) = image_in/scale;
image_out(:,:,3) = image_in/scale;

end


%--------------------------------------

Published with MATLAB® R2020a