findOrtho

Find an orthographic projection relating a set of 3-dimensional points to a corresponding set of noise-corrupted 2-dimensional points.

[P,Normal] = findOrtho(Y,X)

Finds the best-fitting orthographic projection relating a set of 2 dimensional points to corresponding 3 dimensional points where either may be corrupted by noise. It uses constrained optimization to obtain the least-squares solution with an orthogonal projection matrix, up to global rescaling and translation. The output is an affine transformation matrix, P, approximating Y = [X,ones(X,1)]*P.

Inputs:

Y - 2D points as Nx2 matrix

X - 3D points as Nx3 matrixk corresponding to Y

Outputs:

Pout - 3D to 2D affine transformation, approximating

Y = [X,ones(N,1)]*P

where P is an orthographic projection up to rescaling and translation.

Normal - The normal vector to the plane of the projection.

Example:

%Create a random orthographic projection and random data

[Ptrue,~] = svds(randn(3),2);

rsc = randn; %Add random rescaling

Ptrue = Ptrue*rsc; %True projection

Ptrue(4,3) = 1;

X = randn(10,3) ;

X(:,4) = 1;

Y = X*Ptrue(:,1:2) + .1*rsc*randn(size(X,1),2); % Projection with noise

% Find the projection

[P,Normal] = findOrtho(Y(:,1:2),X(:,1:3));

% Plot results

figure

subplot(1,2,1) % Data in Y

scatter(Y(:,1),Y(:,2),50,1:size(X,1),'markerfacecolor','flat')

axis equal

grid on

ax = axis;

title('2D data with noise')

cax = caxis;

subplot(2,2,2) % Data in X

scatter3(X(:,1),X(:,2),X(:,3),50,1:size(X,1),'markerfacecolor','flat')

axis equal vis3d

%Add the true and recovered projection planes

[mgx,mgy] = meshgrid(1:2,1:2);

mgxy = [mgx(:),mgy(:)+2];

axl = ax(mgxy); axl(:,end+1) =1;

axproj = axl*pinv(P);

axproj2 = axl*pinv(Ptrue);

hold on, h1=surf(reshape(axproj(:,1),2,2),reshape(axproj(:,2),2,2),reshape(axproj(:,3),2,2),'facecolor','none'); h2=surf(reshape(axproj2(:,1),2,2),reshape(axproj2(:,2),2,2),reshape(axproj2(:,3),2,2),'facecolor','none','edgecolor','r');

title('3D data')

caxis(cax)

% Same plot but with the camera oriented to the projection plane

subplot(2,2,4)

scatter3(X(:,1),X(:,2),X(:,3),50,1:size(X,1),'markerfacecolor','flat')

hold on

surf(reshape(axproj(:,1),2,2),reshape(axproj(:,2),2,2),reshape(axproj(:,3),2,2),'facecolor','none');

surf(reshape(axproj2(:,1),2,2),reshape(axproj2(:,2),2,2),reshape(axproj2(:,3),2,2),'facecolor','none','edgecolor','r');

axis equal vis3d

caxis(cax)

set(gca,'CameraPosition',Normal*40,'cameraUpVector',P(1:3,2),'cameraTarget',-Normal)

title('Camera oriented to plane')

legend([h1 h2],{'Inferred plane','True plane'},'position',[0.83,0.9,0.1,0.044])