The projection equation from world (x,y,z) to image (u,v) is
We can estimate A from the control points:
% Control points in image coords and world coords
u = [1333 1265 1227 1481 937 842 667];
v = [1066 910 822 499 558 539 532];
x = [ 498899 498895 498892 498854 498734 498682 498584];
y = [4915388 4915389 4915391 4915430 4915441 4915442 4915425];
z = [6.268 6.270 6.264 11.512 4.764 6.698 6.705];
mx=mean(x); my=mean(y); % mean x,y values
xzm=x-mx; yzm=y-my; % x,y with mean removed (zm=zero mean)
% Camera projection matrix
imagePoints=[u',v'];
worldPointszm=[xzm',yzm',z'];
A=estimateCameraProjection(imagePoints,worldPointszm); % camProjection
I removed the mean values from x and y before computing A, to improve numeric accuracy. We can add the mean x,y values back later.
Each red point in the image is at sea level, plus or minus the tide. Let us assume that we know the z value (z=zc) for an unknown world point, and we want to find the x,y values for the unknown world point. Then the equation above becomes
Rearrange equation above:
The last equation above is what we need to convert (u,v) to (x,y), when z=zc is known. If the vector obtained by doing the left hand side above does not have the 3rd element=1, then divide the vector by the third element to get a normalized vector.
I estimated pixel values (u,v) for three points along the red shoreline, and I assume zc=0 along the shoreline. I used the formula above to compute the world coordinates (x,y) of the image points. Script attached. I plotted the green points and the red shoreline points in world coordinates in 3D. The screenshot below shows a view of the points from directly overhead.
