3D (Surface) Reconstruction of Stereo Image

Hello,

i have been trying to reconstruct a 3D surface out of 2 given stereo images. Basically i need to create some kind of surface plot with the texture of a given image at the correct real world coordinates, so i cannot simply use the "texturemap" feature of matlab.

To make it a little bit more clear i try to do something like this:

Where i have the two pictures (from left and righ camera) and i try to create the corresponding surface in 3D-world coordinates. - i need this to visualle check my computations

Can someone please point me in the right direction how to do this ?

thanks

Answer 1

[Agnish Dutta]

Use these instructions as a starting point for what you are trying to accomplish:

The following shows the workflow for creating a 3D point cloud from two stereo images:

https://www.mathworks.com/help/vision/examples/depth-estimation-from-stereo-video.html

The instructions are for video frames, but the same functions apply for singular images as well.

To set the stereo parameters correctly, kindly go through the following:

https://www.mathworks.com/help/vision/ref/stereoparameters.html

You may have to rectify the stereo parameters, so here is the document with the necessary details:

https://www.mathworks.com/help/vision/ref/rectifystereoimages.html

I was able to find something for generating a 3-D surface from a points cloud on MATLAB central:

https://www.mathworks.com/matlabcentral/fileexchange/63731-surface-reconstruction-from-scattered-points-cloud-open-surfaces

[Christoph Betschoga]

Hello Agnish

I would like to thank you for your input. It helped me finde the right direction !

However as i do not have acces to the vision toolbox i had to find another way - i will post it below if someone with as little knowladge of this topic as me stumbles about the same problem ;)

Answer 2

[Christoph Betschoga]

With the help of Agnish i am currently working for a workaround - still in progress:

basic procedure:

1) Find camera matrix for cam0 and cam1 -> are known because stereo system is calibrated

2) i did a mapping workaround for a plane surface - in my case it is always like that:

2a) calculate real coordinates of corners of my rectangular surface -> can be done with known template position in both images when the cameras are calibrated: literature for example: "White Paper: Camera Calibration and Stereo Vision by Peter Hillmann)

2b) define a grid in a reference element and map the gridpoints to the real element (here "shape" functions as linear interpolation between corners are used for this - like in FINITE ELEMENTS, as i am a civil engineer ;) )

2c) then with the known perspective projection matrix of for instance cam0 define the location of the real world coordinates in the picture - one more mapping

2d) interpolate the colors at these positions (picture values)

2e) plot the surface

Example Code should work in case someone is interested, picture can be downloaded here: https://postimg.cc/qtGgsGVn.

I suppose this methode can somewho be generalized with some triangulation etc. but somehow the biggest problem in my case was the recognition of equivalemt templates in both picture, which is necessary to get the real world coordinates. I could not get it to work .... also tried rectification and so on but without any success.

greetings all

% info: distortion neglected
% CAM0 
% intrinsics
KK0 = [ 9.8490634765625000e+003     0.                          2.3789833984375000e+003 ;...
        0.                          9.8490634765625000e+003     1.7011253662109375e+003 ; ...
        0.                          0.                          1.                      ] ;

% extrinsics
rotvec0 = [ 0 0 0 ] ; 
transvec0 = [ 0 0 0 ].' ;

% camera matrix
PP0 = KK0 * [ RotMatrix(0, [1,0,0]), transvec0 ] ; %rodr. func. by Author: Jan Simon, Heidelberg
% no rotation and translation in CAM0 - hence (0,...

% CAM1 
% intrinsics
KK1 = [ 7.1750791015625000e+003     0.                          2.4779909667968750e+003 ; ...
        0.                          7.1750791015625000e+003     1.5908875732421875e+003 ; ...
        0.                          0.                          1.                      ] ;

% extrinsics
rotvec1 =   [ -1.4065706325327406e-002 1.3712001490606401e-001   2.2151575978542927e-003 ] ; 
transvec1 = [ -3.6780345263841093e+002 2.8410817957758999e-001  -5.4874311039726160e+002 ].' ;

% camera matrix
PP1 = KK1 * [ RotMatrix(norm(rotvec1), rotvec1), transvec1 ] ;


% for functional example
xkoor =     [  -392.555    394.961     394.704   -398.128 ] ;
ykoor =  -1*[  -043.410   -060.136     329.221    313.227 ] ;
zkoor =  -1*[ -2287.690  -2234.0191  -2223.241  -2268.321 ] ;

% grid points in x and y direction
dx =1904; dy =950;
xi = linspace(-1,1,dx) ; eta = linspace(-1,1,dy);

% define mapping from reference element (-1/-1 and -1/1) to real element 
% with xkoord, ykoor, zkoor
c_x(1) = sum(xkoor)/4;
c_x(2) = -(xkoor(1) - xkoor(2)- xkoor(3) + xkoor(4))/4;
c_x(3) = -(xkoor(1) + xkoor(2)- xkoor(3) - xkoor(4))/4;
c_x(4) = +(xkoor(1) - xkoor(2)+ xkoor(3) - xkoor(4))/4;

c_y(1) = sum(ykoor)/4;
c_y(2) = -(ykoor(1) - ykoor(2)- ykoor(3) + ykoor(4))/4;
c_y(3) = -(ykoor(1) + ykoor(2)- ykoor(3) - ykoor(4))/4;
c_y(4) = +(ykoor(1) - ykoor(2)+ ykoor(3) - ykoor(4))/4;

c_z(1) = sum(zkoor)/4;
c_z(2) = -(zkoor(1) - zkoor(2)- zkoor(3) + zkoor(4))/4;
c_z(3) = -(zkoor(1) + zkoor(2)- zkoor(3) - zkoor(4))/4;
c_z(4) = +(zkoor(1) - zkoor(2)+ zkoor(3) - zkoor(4))/4;

% grid on reference element
[X,Y] = meshgrid(xi,eta) ;
Z = X.*Y ;

% do the mapping 
c_mat = [c_x ; c_y ; c_z] ;
matmult = [ones(1,length(Y(:))); Y(:).'; X(:).'; Z(:).'] ;

% interpolated coordinates in world coordinate system
mapped = c_mat * matmult ;

% corresponding grid points to plot surface
surf_gridX = reshape(mapped(1,:),[dy dx]).'; 
surf_gridY = reshape(mapped(2,:),[dy dx]).'; 
surf_gridZ = reshape(mapped(3,:),[dy dx]).'; 

% map real coordinates with camera projection matrix to pixel coordinates
mapped_pic0 = PP0*[mapped ; ones(1,length(Y(:))) ];
mapped_pic0 = mapped_pic0(1:2,:) ./ mapped_pic0(3,:) ;

% read picture
[I] = imread('Frame-004635.png') ;

% interpolate for mapped pixel positions for texture
test=uint8(interp2(1:size(I,2),1:size(I,1),single(flipud(I(:,:))),mapped_pic0(1,:),mapped_pic0(2,:),'cubic',0)) ;
surf_color = reshape(test,[dy dx]).' ;

% plot figure
figure
ts = surf(surf_gridX,surf_gridY,surf_gridZ,surf_color,'EdgeColor','none','FaceColor','interp')
colormap(gray)
view(-18,71)