test_invperspmapsetup.m

逆透視轉換Matlab程式

% Script to test the inverse perspective mapping setup and processing
% functions, pixelsToWorld, getInterpMap, and getWorldImage.
%
% Author:
%----------
% Eric Johnson
% University of Utah
% CS 5320/6320 - Computer Vision 
% April 14, 2007

clear; clc;

%------------------
%% pixelsToWorld.m
%------------------
% Put in the parameters for the DARPA web cam images (since that's still
% all we have to work with).  The values are those which were developed in
% the invPerspCalibrate.m script.
camera.name = 'Web Cam';
camera.m = 480; % Rows (height)
camera.n = 640; % Columns (width)
camera.h = 1.5; % Distance camera was above the ground (meters)
camera.alphaTot = 30*pi/180; % Total HALF viewing angle (corner to corner - radians)
camera.theta0 = -0.01373524283527; % Camera tilt angle below horizontal (radians)
% That's correct, a negative sign.  The camera was pointed slightly up
% relative to horizontal.

% Start with the pixel to world mapping function.
% Keep track of how long this takes to run.  Sonce it turns out to be
% really fast, average over 10 times.
tic;
for i = 1:10
    [xMap, yMap] = pixelsToWorld(camera);
end
t = toc;
t = t/10;
disp('Time required to set up the pixel to world map:');
disp(sprintf('\t%.3f seconds', t));

% Plot the xMap versus yMap to make sure that the points fan out as
% expected.  Downsample them a bit so the grid will be sparse enough to
% tell what is going on.
hF1 = figure('Name', 'xMap, yMap', 'Units', 'pixels', 'Position', ...
    [50 50 640 480]);
sampby = 50;
xPlot = xMap(:, 1:sampby:end);
yPlot = yMap(:, 1:sampby:end);
% Plot all of the columns as lines, then all of the rows.
rPlot = size(xPlot,1);
cPlot = size(xPlot,2);
for r = 1:rPlot
    line(xPlot(r, :), yPlot(r, :), 'LineStyle', '-', 'Color', 'b', ...
        'Marker', '.', 'MarkerSize', 5);
end
for c = 1:cPlot
    line(xPlot(:,c), yPlot(:,c), 'LineStyle', '-', 'Color', 'b', ...
        'Marker', '.', 'MarkerSize', 5);
end
title('Image Pixels Mapped to World Frame', 'FontSize', 12);
xlabel('x (meters)', 'FontSize', 12);
ylabel('y (meters)', 'FontSize', 12);
axis tight;
axis equal;
drawnow;

%-----------------
%% getInterpMap.m
%-----------------

% Now set up the desired xRange, yRange and step based on our observations
% about the usable area of the IPM image from the calibration images.
params.xRange = [min(xMap(:,1)), 25];
params.yRange = [-8 8];
params.mIPM = 300;
params.step = (params.yRange(2) - params.yRange(1))/(params.mIPM-1);
% Keep track of how long this takes to run. (It's slow, so only do it once.)
tic;
[xGrid, yGrid, interpMap] = getInterpMap(xMap, yMap, camera, params);
t = toc;
disp('Time required to set up the interpolation map:');
disp(sprintf('\t%.3f seconds', t));
data.interpMap = interpMap;
data.xGrid = xGrid;
data.yGrid = yGrid;

% Load up an image to test this on and run it through the getWorldImage
% function.
I = imread('roadCalPic.ppm');

% Again, we want an idea of how fast it runs, so keep track of the elapsed
% time.  It turns out to be pretty quick, so do it 10 times in a row and
% average.
tic;
for i = 1:10
    [Wc, Wg] = getWorldImage(I, data);
end
t = toc;
t = t/10;
disp('Time required to apply the interpolation map for an RGB image:');
disp(sprintf('\t%.3f seconds', t));

% Display the resulting IPM images and set the axes to show the world
% coordinates instead of the rows and columns.  Also show the original
% image for comparison.
hF2 = figure('Name', 'IPM results', 'Units', 'inches', 'Position', ...
    [1 1 6.5 6]);

gap = 0.075;
fs.title = 10;
fs.label = 9;
fs.axes = 8;
width = (1-4*gap)/2;
height = (1-4*gap)/2;
top = (1-4*gap)/2+3*gap;
bot = gap;
left = gap;
right = (1-4*gap)/2+3*gap;

hA = axes('Units', 'normalized', 'Position', ...
    [left top width height], ...
    'FontSize', fs.axes);
image(I);
set(hA, 'DataAspectRatio', [1 1 1]);
title('Original Image', 'FontSize', fs.title);
xlabel('Column', 'FontSize', fs.label);
ylabel('Row', 'FontSize', fs.label);

hA = axes('Units', 'normalized', 'Position', ...
    [right top width height], ...
    'FontSize', fs.axes);
image(Wc, 'XData', params.xRange, 'YData', fliplr(params.yRange));
set(hA, 'YDir', 'normal', 'DataAspectRatio', [1 1 1]);
title('Inverse-Perspective-Mapped Image', 'FontSize', fs.title);
xlabel('x = Forward Distance (meters)', 'FontSize', fs.label);
ylabel('y = Side-to-Side Distance (meters)', 'FontSize', fs.label);
drawnow;

% To get the grayscale ones to display properly on the same figure, we need
% to create equivalent RGB versions with the intensities repeated equally
% in each channel.
hA = axes('Units', 'normalized', 'Position', ...
    [left bot width height], ...
    'FontSize', fs.axes);
Ig = mean(double(I),3)/255;
image(repmat(Ig, [1, 1, 3]));
set(hA, 'DataAspectRatio', [1 1 1]);
title('Original Image', 'FontSize', fs.title);
xlabel('Column', 'FontSize', fs.label);
ylabel('Row', 'FontSize', fs.label);

hA = axes('Units', 'normalized', 'Position', ...
    [right bot width height], ...
    'FontSize', fs.axes);
image(repmat(Wg, [1, 1, 3]), 'XData', params.xRange, ...
    'YData', fliplr(params.yRange));
set(hA, 'YDir', 'normal', 'DataAspectRatio', [1 1 1]);
title('Inverse-Perspective-Mapped Image', 'FontSize', fs.title);
xlabel('x = Forward Distance (meters)', 'FontSize', fs.label);
ylabel('y = Side-to-Side Distance (meters)', 'FontSize', fs.label);
drawnow;