function [output, n, i, image, baseline, scandata, laserdata] = calibrate(laserscans, tstep, tbase);

% No initial input necessary.
% Name baseline image "cbase.jpg", and the scan images
% "cscan_1.jpg",...,"cscan_n.jpg"

baseline = convert('cbase.jpg');
i = 0;
n = size(laserscans,2); 
camstats = auto_calib_cam2(baseline, gray);
tlaserdata = zeros(n,3);
laserdata = zeros(3,4);
ymax = size(baseline,1);

lines = zeros(n,2);
avgpoint = zeros(3,n);

while i<n,
	i = i+1;
	scandata = laserscans(:,i);
	scandata = [scandata (1:size(scandata,1))'];
	scandata = scandata(find(~isnan(sum(scandata,2))),:);
	points = pixel_2_3dpoint(scandata',camstats(6,1:2)', camstats(5,1:2)', camstats(7,1), rodrigues(camstats(1,:)), camstats(2,:)');
	avgpoint(:,i) = mean(points, 2);
	%find X as a function of Z
	lines(i,:) = polyfit(points(3,:), points(1,:), 1);
	
	%points(1:2,:) = points(1:2,:) ./ [points(3,:); points(3,:)];
	%tlaserdata(i,:) = polyfit(points(2,:), points(1,:), 2);
end;

%find (X,Z) intersection of lines
%find Z which minimizes the stdev of all the lines.
Z = fmin('stdevlines', -1000, 1000, [], lines)
%therefore X
X = meanlines(Z, lines)

for i = 1:n
	%find theta of planes
	theta(i) = atan2(avgpoint(1,i)-X,avgpoint(3,i)-Z)
end
%fit to polynomial
thetap = polyfit(linspace(tbase, tbase+tstep*(n-1), n),theta,3);


%tsteps = ((0:n-1)*tstep) + tbase;
%laserdata(1,:) = polyfit(tsteps', tlaserdata(:,1), 3);
%laserdata(2,:) = polyfit(tsteps', tlaserdata(:,2), 3);
%laserdata(3,:) = polyfit(tsteps', tlaserdata(:,3), 3);
output = zeros(10,4);
output(1,1:4) = thetap;
output(4:10,1:3) = camstats;

output(2,1) = X;
output(2,2) = Z;