Two view Demo: an example

http://vision.ucla.edu//MASKS/twoview_reconstruction/TwoViewDemo.m
On the link above some sample code in matlab for fulfilling some triangulation operations.
% two view example of calibrated structure from motion
% recovery
close all; clear;

im0 = rgb2gray(imread('boxes1.bmp','bmp'));
im1 = rgb2gray(imread('boxes8.bmp', 'bmp'));
[ydim, xdim] = size(im0);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% intrinsic parameter matrix
Khat = [653.91 0 320.52; 0 649.72 220.40; 0 0 1];

%%%%%%%%%%%%%%%%%%%
% undo distortion
im0 = distortion(im0,0.03);
im1 = distortion(im1,0.03);

%%%%%%%%%%%%%%%%%%%%%
load boxpoints  % 13 points seen in two frames frames
NPOINTS =  size(xim(:,:,1),2);  

figure(1); imagesc(im0); colormap gray; hold on; axis equal; axis off;
plot(xim(1,:,1), xim(2,:,1),'w.');
figure(2); imagesc(im1); colormap gray; hold on; axis equal; axis off;
plot(xim(1,:,2), xim(2,:,2),'w.');


FRAMES = 2;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% compute retinal coordinates
for i = 1:FRAMES
 xn(:,:,i) = inv(Khat)*xim(:,:,i);
end

%%%%%%%%%%%%%%%%%%%%
% initialize motion
[T0, R0] = dessential(xn(:,:,1), xn(:,:,2));

lambda = zeros(1,NPOINTS);
XE = zeros(3,NPOINTS,FRAMES);
[X,lambda]=compute3DStructure(xn(:,:,1),xn(:,:,2), R0, T0);
figure(FRAMES + 1);
XE(:,:,1) = X(:,:,1);
XE(:,:,2) = X(:,:,2);
plot3(XE(1,:,1),XE(2,:,1),XE(3,:,1),'r.'); hold on;

% pause
[l3d, inc]  = boxseqlines(XE(1:3,:,1));
lnum = size(l3d,2);
for i=1:lnum
  line([l3d(1,i) l3d(4,i)], [l3d(2,i) l3d(5,i)], [l3d(3,i) l3d(6,i)]);
end;
xlabel('x'); ylabel('y'); zlabel('z'); % axis equal;
title('Euclidean reconstruction'); view(220,20); axis equal; box on;
hold on;
pause;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% compute residuals and reproject them
xres0=Khat*[X(1,:,1)./X(3,:,1); X(2,:,1)./X(3,:,1); ones(1,NPOINTS)];
xres1=Khat*[X(1,:,2)./X(3,:,2); X(2,:,2)./X(3,:,2); ones(1,NPOINTS)];
figure(1); hold on; plot(xres0(1,:), xres0(2,:),'y.'); 
figure(FRAMES); hold on; plot(xres1(1,:), xres1(2,:),'y.');
pause;

xd0 = xim(1:2,:,1) - xres0(1:2,:);
xd1 = xim(1:2,:,2) - xres1(1:2,:);
res = [xd0'; xd1'];
f = [norm(res'*res)/2];
% pause
% two view example of calibrated structure from motion
% recovery
close all; clear;

im0 = rgb2gray(imread('boxes1.bmp','bmp'));
im1 = rgb2gray(imread('boxes8.bmp', 'bmp'));
[ydim, xdim] = size(im0);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% intrinsic parameter matrix
Khat = [653.91 0 320.52; 0 649.72 220.40; 0 0 1];

%%%%%%%%%%%%%%%%%%%
% undo distortion
im0 = distortion(im0,0.03);
im1 = distortion(im1,0.03);

%%%%%%%%%%%%%%%%%%%%%
load boxpoints  % 13 points seen in two frames frames
NPOINTS =  size(xim(:,:,1),2);  

figure(1); imagesc(im0); colormap gray; hold on; axis equal; axis off;
plot(xim(1,:,1), xim(2,:,1),'w.');
figure(2); imagesc(im1); colormap gray; hold on; axis equal; axis off;
plot(xim(1,:,2), xim(2,:,2),'w.');


FRAMES = 2;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% compute retinal coordinates
for i = 1:FRAMES
 xn(:,:,i) = inv(Khat)*xim(:,:,i);
end

%%%%%%%%%%%%%%%%%%%%
% initialize motion
[T0, R0] = dessential(xn(:,:,1), xn(:,:,2));

lambda = zeros(1,NPOINTS);
XE = zeros(3,NPOINTS,FRAMES);
[X,lambda]=compute3DStructure(xn(:,:,1),xn(:,:,2), R0, T0);
figure(FRAMES + 1);
XE(:,:,1) = X(:,:,1);
XE(:,:,2) = X(:,:,2);
plot3(XE(1,:,1),XE(2,:,1),XE(3,:,1),'r.'); hold on;

% pause
[l3d, inc]  = boxseqlines(XE(1:3,:,1));
lnum = size(l3d,2);
for i=1:lnum
  line([l3d(1,i) l3d(4,i)], [l3d(2,i) l3d(5,i)], [l3d(3,i) l3d(6,i)]);
end;
xlabel('x'); ylabel('y'); zlabel('z'); % axis equal;
title('Euclidean reconstruction'); view(220,20); axis equal; box on;
hold on;
pause;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% compute residuals and reproject them
xres0=Khat*[X(1,:,1)./X(3,:,1); X(2,:,1)./X(3,:,1); ones(1,NPOINTS)];
xres1=Khat*[X(1,:,2)./X(3,:,2); X(2,:,2)./X(3,:,2); ones(1,NPOINTS)];
figure(1); hold on; plot(xres0(1,:), xres0(2,:),'y.'); 
figure(FRAMES); hold on; plot(xres1(1,:), xres1(2,:),'y.');
pause;

xd0 = xim(1:2,:,1) - xres0(1:2,:);
xd1 = xim(1:2,:,2) - xres1(1:2,:);
res = [xd0'; xd1'];
f = [norm(res'*res)/2];
% pause