function model_lidar_and_validate_tform
%MODEL_LIDAR_AND_VALIDATE_TFORM encapsulates the overall process
% involved with modeling Lidar, developing a spatial transform
% to correct as-scanned image distortion, and validate the
% correction algorithm using a synthetic test image.
%
%USAGE: model_lidar_and_validate_tform
%Subfunctions: MAKE_LIDAR_TFORM
% MAKE_POL2CART_TFORM
% MAKE_JI2RT_TFORM
% Copyright 2004-2010 RBemis The MathWorks, Inc.
%% Step 1. Get test image (5x5 array of clouds) & define coordinates
im=imread('cloud_array.tif'); % read image from file
A=15*pi/180; % sweep amplitude (rad)
d=[10 20]; % distance range (m)
[M,N]=size(im); % scan image size (pixels)
% first display original image
x1=d(1)*cos(A); x2=d(2); dx=(x2-x1)/(N-1); xx=x1:dx:x2;
y2=d(2)*sin(A); y1=-y2; dy=(y2-y1)/(M-1); yy=y1:dy:y2;
figure, image(xx,yy,im), colormap gray
title('Original test image')
xlabel('Horizontal Distance (m)')
ylabel('Vertical Height (m)')
% then overlay nonuniform, polar scan coordinate system
m=11; n=11; % number scan angles & distances
i=[1:m]'; j=1:n; % angle(row) & distance(col) indices
[I,J]=meshgrid(i,j); % grid system representing scan points
U=[J(:) I(:)]; % scan points (pixel coordinates)
t1=make_ji2rt_tform(A,d,m,n); % map indices (I,J) <-> polar (R,TH)
V=tformfwd(U,t1); % scan points in polar coordinates
R=V(:,1); TH=V(:,2); % component values
t2=make_pol2cart_tform; % map polar (R,TH) <-> cartesian (X,Y)
W=tformfwd(V,t2); % scan points in cartesian coordinates
X=reshape(W(:,1),n,m); % X components
Y=reshape(W(:,2),n,m); % Y components
line_width=1.5; set(gca,'color','k')
line(X,Y,'color','g','linewidth',line_width)
line(X',Y','color','g','linewidth',line_width)
title('Original test image + scan grid')
axis equal, ax1=axis; shg
%f=getframe; imwrite(f.cdata,'Fig1.tif','tiff')
%% Step 2. Scan image = inverse transform of test image
%% (how scanner sees it)
T=make_lidar_tform(A,d,M,N);
im_inv=imtransform(im,fliptform(T),'xdata',[1 N],'ydata',[1 M]);
a=A*[-1:2/(m-1):1]';
x1=min(X(:)); x2=max(X(:)); % horizontal limits
dx=(x2-x1)/(n-1); % horizontal spacing
y1=min(Y(:)); y2=max(Y(:)); % vertical limits
dy=(y2-y1)/(m-1); % vertical spacing
x=x1:dx:x2; y=[y1:dy:y2]'; % horz/vert indices
figure, image(x,a,im_inv), colormap gray
title('Scan image (clouds as detected by lidar system)')
ylabel('Scan Angle (rad)'), xlabel('Radial Distance (m)')
set(gca,'color','k'), xlim(ax1(1:2))
[X2,Y2]=meshgrid(x,y); % similar grid
A3=repmat(a,[1 n]);
line(X2,A3,'color','g','linewidth',line_width)
line(X2',A3','color','g','linewidth',line_width)
%f=getframe; imwrite(f.cdata,'Fig2.tif','tiff')
%% Step 3. Correct image = forward transform of scan image
%% (illustrates distortion removal)
im_inv_fwd=imtransform(im_inv,T,'xdata',[1 N],'ydata',[1 M],'FillValues',0);
figure, image(x,y,im_inv_fwd), colormap gray, axis equal
title('Apply transform of correct scan distortion')
xlabel('Horizontal Distance (m)'), ylabel('Vertical Height (m)')
set(gca,'color','k')
line(X,Y,'color','g','linewidth',line_width)
line(X',Y','color','g','linewidth',line_width)
%f=getframe; imwrite(f.cdata,'Fig3.tif','tiff')
%% Step 4. Determine error between original & corrected images
%% (shows actual differences)
im_diff=imsubtract(im,im_inv_fwd);
figure, image(x,y,im_diff), colormap gray, axis equal
title('Difference between Original & Processed Images')
xlabel('Horizontal Distance (m)'), ylabel('Vertical Height (m)')
set(gca,'color','k')
%f=getframe; imwrite(f.cdata,'Fig4.tif','tiff')
%% NOTE: The clouds outside the scanner's view obviously get ignored.
%% Real differences are due to edge interpolation, which is
%% much less pronounced by comparison as expected.
%% CONCLUSION: The correction worked very well!
figure(gcf-3)
return
% End of main function (model_lidar_and_validate_tform)
% ----------------------------------------------------------------------
% NOTES TO PRESENTER
%
% Originally this was a plain script. As a function however it
% runs in its own workspace. This has advantages when presenting
% the Lidar case study to an audience (ie, 2002 AIP seminars).
%
% Modeling Lidar to develop a TFORM was only 1 intermediate step
% in the overall story. This diversion into model development
% and spatial transform validation is a tangential story which
% could be skipped to shorten the presentation. Or it could be
% its own topic talk (application of spatial transform to correct
% nonuniform scan distortion and to convert acquired images from
% polar to cartesian coordinates).
%
% If time permits, step through the code to explain what MATLAB
% actually did in fine detail; however, if time is tight, just
% run this (wham) and talk from the figures using the following
% explanations as a guide.
% ----------------------------------------------------------------------
% FIGURE EXPLANATION
%
% 1. The original test image is a regular pattern. The overlaid
% grid shows how the scanner collects data (nonuniformly in
% polar coordinates).
%
% 2. How the scanner sees the test pattern. Notice the non-
% uniform and curvelinear relationship or mapping
% between the two coordinate systems. (flip back & forth)
%
% 3. The spatial transform corrects scan images for distortion.
% Notice it crops off scenery outside its field of view.
%
% 4. How well did it really work? Compare actual differences
% between the original and corrected images. The obvious
% difference is cropped off clouds outside the working
% envelope of the scanner. More subtle edge effects around
% each cloud are due to interpolation. These cartoon clouds
% had particularly strong contrast between white and black
% around their edges. Real vapor clouds tend to have softer
% edges due to diffusion, so the interpolation artifact is
% only a minor concern.
% ----------------------------------------------------------------------
function [T,x,y]=make_lidar_tform(A,d,m,n)
%MAKE_LIDAR_TFORM creates TFORM object suitable for use in
% mapping between lidar scan images and cartesian space.
%
%USAGE: [T,x,y]=make_lidar_tform(A,d,m,n)
%
%INPUTS: A = maximum amplitude of oscillation (rad)
% d = [d1 d2] (1x2) vector if near & far distances (eg, m)
% m = number rows (input & output images same size)
% n = number cols (ditto)
%
%OUTPUTS: T = TFORM object (composite)
% x = column index vector (distance units)
% y = row index vector (distance units)
if nargin~=4
error('USAGE: [T,x,y]=make_lidar_tform(A,d,m,n)')
end
% key points for original image (pixel coordinates)
i=[1:m]'; j=1:n; [I,J]=meshgrid(i,j);
% key points for original image (spatial coordinates)
ctr=mean(i); period=2*(m-1); th=A*sin(2*pi*(i-ctr)/period);
r1=d(1); dr=diff(d)/(n-1); r=r1+(j-1)*dr; [R,TH]=meshgrid(r,th);
% key points for warped image (spatial coordinates)
[X,Y]=pol2cart(TH,R);
% key points for resampled image (spatial coordinates)
x1=min(X(:)); x2=max(X(:)); dx=(x2-x1)/(n-1);
y1=min(Y(:)); y2=max(Y(:)); dy=(y2-y1)/(m-1);
x=x1:dx:x2; y=[y1:dy:y2]'; [X2,Y2]=meshgrid(x,y);
t1=make_ji2rt_tform(A,d,m,n);
t2=make_pol2cart_tform;
t3=maketform('box',[X2([1 end]') Y2([1 end]')],[J([1 end]') I([1 end]')]);
T=maketform('composite',t3,t2,t1);
% ----------------------------------------------------------------------
function t = make_ji2rt_tform(A,d,m,n)
%MAKE_JI2RT_TFORM creates TFORM object suitable for use with
% lidar scan images (see MAKE_LIDAR_TFORM). This TFORM maps
% between J-I (image) indices and R-TH (scan) coordinates.
%
%USAGE: T=make_ji2rt_tform(A,d,m,n)
%
%INPUTS: A = maximum amplitude of oscillation (rad)
% d = [d1 d2] (1x2) vector if near & far distances (eg, m)
% m = number rows (input & output images same size)
% n = number cols (ditto)
%
%OUTPUT: T = TFORM object (custom)
i=[1:m]'; j=1:n; % row/col indices
ctr=mean(i); period=2*(m-1); % vertical constants
r1=d(1); dr=diff(d)/(n-1); % horizontal terms
% bundle scan parameters into structure (tdata)
Params.A=A; % amplitude of oscillation (rad)
Params.ctr=ctr; % middle row = X axis
Params.period=period; % period = twice image/frame height
Params.r1=r1; % nearest scan distance (ie, m)
Params.dr=dr; % radial distance range (ie, m)
% define custom transform
t=maketform('custom',2,2,@ji2rt,@rt2ji,Params);
% ----------------------------------------------------------------------
function V=ji2rt(U,t)
% FORWARD sine transform implementation
Params=t.tdata; % scan system parameters
A=Params.A; % amplitude of oscillation (rad)
ctr=Params.ctr; % middle row = X axis
period=Params.period; % period = twice image/frame height
r1=Params.r1; % nearest scan distance (ie, m)
dr=Params.dr; % radial distance range (ie, m)
j=U(:,1); i=U(:,2); % inputs (spatial coordinates)
th=A*sin(2*pi*(i-ctr)/period); % scan angles (deg)
r=r1+(j-1)*dr; % radial distances (ie, m)
V=[r th]; % outputs (spatial coordinates)
% ----------------------------------------------------------------------
function U=rt2ji(V,t)
% INVERSE sine transform implementation
Params=t.tdata; % scan system parameters
A=Params.A; % amplitude of oscillation (rad)
ctr=Params.ctr; % middle row = X axis
period=Params.period; % period = twice image/frame height
r1=Params.r1; % nearest scan distance (ie, m)
dr=Params.dr; % radial distance range (ie, m)
r=V(:,1); th=V(:,2); % inputs (spatial coordinates)
j=(r-r1)/dr+1; % col
i=asin(th/A)*period/2/pi+ctr; % row
%NOTE: The ASIN function is only valid for inputs in the
% interval [-1,1]. Invalid inputs are detected and
% replaced with NaN (not a number) by the following line.
i(abs(th)>A)=NaN; % ignore invalid row values
U=[j i]; % outputs (spatial coordinates)
% ----------------------------------------------------------------------
function T=make_pol2cart_tform
%MAKE_POL2CART_TFORM creates TFORM object suitable for use with
% lidar scan images (see MAKE_LIDAR_TFORM). This TFORM maps
% between R-TH (polar) and X-Y (cartesian) coordinates.
%USAGE: T=make_pol2cart_tform
%
%INPUTS: none
%
%OUTPUT: T = TFORM object (custom)
T=maketform('custom',2,2,@pol2cart_fwd,@pol2cart_inv,[]);
% ----------------------------------------------------------------------
function V=pol2cart_fwd(U,t)
% FORWARD polar-cartesian transform implementation
r=U(:,1); % input #1 = radius (eg, m)
th=U(:,2); % input #2 = angle (rad)
[x,y]=pol2cart(th,r); % polar => cartesian conversion
V=[x y]; % outputs (spatial coordinates)
% ----------------------------------------------------------------------
function U=pol2cart_inv(V,t)
% INVERSE polar-cartesian transform implementation
x=V(:,1); y=V(:,2); % inputs (spatial coordinates)
[th,r]=cart2pol(x,y); % cartesian => polar conversion
U=[r th]; % output #1 = radius (ie, m)
% output #2 = angle (rad)
