Here is official page:
https://pypi.org/project/ip2geotools/
https://github.com/tomas-net/ip2geotools
This is my try
..
..
The results :
165.132.116.105
Seodaemun-gu
Seoul
KR
37.5790747
126.9367861
Thank you. ๐๐ป♂️
www.marearts.com
%//////////////////////////////////////////////////////////////////////////
%// Made by J.H.KIM, 2011 / feelmare@daum.net, feelmare@gmail.com //
%// blog : http://feelmare.blogspot.com //
%// Eight-Point Algorithm
%//////////////////////////////////////////////////////////////////////////
clc; clear all; close all;
% Corresponding points between two images
% sample #1 I11.jpg, I22.jpg
%{
load I11.txt; load I22.txt;
m1 = I11; m2 = I22;
%}
%sample #2 I1.jpg, I2.jpg
load I1.txt; load I2.txt;
m1 = I1; m2 = I2;
s = length(m1);
m1=[m1(:,1) m1(:,2) ones(s,1)];
m2=[m2(:,1) m2(:,2) ones(s,1)];
Width = 800; %image width
Height = 600; %image height
% Intrinsic Matrix
load intrinsic_matrix.txt
K = intrinsic_matrix;
% The matrix for normalization(Centroid)
N=[2/Width 0 -1;
0 2/Height -1;
0 0 1];
%%
% Data Centroid
x1=N*m1'; x2=N*m2';
x1=[x1(1,:)' x1(2,:)'];
x2=[x2(1,:)' x2(2,:)'];
% Af=0
A=[x1(:,1).*x2(:,1) x1(:,2).*x2(:,1) x2(:,1) x1(:,1).*x2(:,2) x1(:,2).*x2(:,2) x2(:,2) x1(:,1) x1(:,2), ones(s,1)];
% Get F matrix
[U D V] = svd(A);
F=reshape(V(:,9), 3, 3)';
% make rank 2
[U D V] = svd(F);
F=U*diag([D(1,1) D(2,2) 0])*V';
% Denormalize
F = N'*F*N;
%Verification
%L1=F*m1'; m2(1,:)*L1(:,1); m2(2,:)*L1(:,2); m2(3,:)*L1(:,3);
%%
%Get E
E=K'*F*K;
% Multiple View Geometry 259page
%Get 4 Possible P matrix
P4 = get4possibleP(E);
%Get Correct P matrix
inX = [m1(1,:)' m2(1,:)'];
P1 = [eye(3) zeros(3,1)];
P2 = getCorrectCameraMatrix(P4, K, K, inX)
%%
%Get 3D Data using Direct Linear Transform(Linear Triangular method)
Xw = Triangulation(m1',K*P1, m2',K*P2);
xx=Xw(1,:);
yy=Xw(2,:);
zz=Xw(3,:);
figure(1);
plot3(xx, yy, zz, 'r+');
%{
%This code is also run well instead of Triangulation Function.
nm1=inv(K)*m1';
nm2=inv(K)*m2';
% Direct Linear Transform
for i=1:s
A=[P1(3,:).*nm1(1,i) - P1(1,:);
P1(3,:).*nm1(2,i) - P1(2,:);
P2(3,:).*nm2(1,i) - P2(1,:);
P2(3,:).*nm2(2,i) - P2(2,:)];
A(1,:) = A(1,:)./norm(A(1,:));
A(2,:) = A(2,:)./norm(A(2,:));
A(3,:) = A(3,:)./norm(A(3,:));
A(4,:) = A(4,:)./norm(A(4,:));
[U D V] = svd(A);
X(:,i) = V(:,4)./V(4,4);
end
%}
%% clc; clear all; close all; %% image coordinate m1 = [ -88.6 1019.21739130435 1046.63913043478 -9643.52173913043 1488.59518796992; -216 -14137.7391304348 -629.04347826087 18775.3043478261 -464.421052631579]; m2 = [ 644.916058797683 8264.25802702142 2735.93716970091 -3264.00791696585 4601.62993487106; 237.341114262516 -16276.926674065 -591.505245744724 4076.20064528107 -313.261770012357]; % rotation matrix RealR = rodrigues( [-10 30 20]*pi/180 ); % translation matrix RealT = [20 30 200]'; RealA=rodrigues(RealR)*180/pi; % rotation matrix angle Test K = [1000 1 512;0 1000 384;0 0 1;]; %calibration matrix %% P1, P2 ๋ง๋ค๊ธฐ P1 = K*[eye(3) zeros(3,1)]; % left camera Projection matrix P2 = K*[RealR RealT]; % right camera projection matrix %% 3์ฐจ์ ์ ๋ง๋ค๊ธฐ %P1, P2๋ฅผ ์ด์ฉํ 3์ฐจ์ ์ ๋ณต์ % World 3D coordination reconstruction by Direct Linear Transformation % Multiple view geometry P.312 W=[]; for i=1:5 A=[ m1(1,i)*P1(3,:) - P1(1,:); m1(2,i)*P1(3,:) - P1(2,:); m2(1,i)*P2(3,:) - P2(1,:); m2(2,i)*P2(3,:) - P2(2,:)]; A(1,:) = A(1,:)/norm(A(1,:)); A(2,:) = A(2,:)/norm(A(2,:)); A(3,:) = A(3,:)/norm(A(3,:)); A(4,:) = A(4,:)/norm(A(4,:)); [u d v] = svd(A); W=[W v(:,4)/v(4,4)]; end %% ๋ค์ 3์ฐจ์ ์ ์์ ํฝ์ ์ ๋ง๋ค๊ธฐ % again calculate 2D image point from World coordinate, projection matrix reip1 = P1*W; reip1 = [reip1(1,:)./reip1(3,:); reip1(2,:)./reip1(3,:)] m1(:,1:5) reip2 = P2*W; reip2 = [reip2(1,:)./reip2(3,:); reip2(2,:)./reip2(3,:)] m2(:,1:5) % m1 and reip1 is same. Also, m2 and reip2 is same.
%% clc; clear all; close all; %% ๋ฐ์ดํฐ ๋ก๋ฉ %% Data Loading load rotation_matrices.txt load translation_vectors.txt load intrinsic_matrix.txt load distortion_coeffs.txt R = rotation_matrices; %R1์ Left์นด๋ฉ๋ผ์ ํจํดํ์ผ๋ก ๋ถํฐ์ Rotation Matrix %R1 is Rotation Matrix of Left Camera from Pattern board. Pattern Board has origin %coordinate (0,0,0) R1 = reshape(R(19,:),3,3); %R2 is Left Camera Rotation Matrix. R2 = reshape(R(20,:),3,3); %T1 is Left Camera Translation Matrix T1 = translation_vectors(19,:)'; %T1 is right Camera Translation Matrix T2 = translation_vectors(20,:)'; K = intrinsic_matrix; %Load Matched Coordinate load pattern19.txt; load pattern20.txt; m1 = pattern19'; m2 = pattern20'; %% Real R,T ๋ง๋ค๊ธฐ %% Make Real R,T that is relation between Left, Right Camera % R,T is made by R1,R2 and T1, T1 RealT = T2-T1; %T๋ ๊ทธ๋ฅ ๋นผ๋ฉด ๋๋ค. RealR = R2*inv(R1); RealA=rodrigues(RealR)*180/pi; %This is Angle %% P1, P2 ๋ง๋ค๊ธฐ % Make Projection matrix P1 = K*[eye(3) zeros(3,1)]; %P1 is reference Camera so P1=K[I:O] P2 = K*[RealR RealT]; %% 3์ฐจ์ ์ ๋ง๋ค๊ธฐ %P1, P2๋ฅผ ์ด์ฉํ 3์ฐจ์ ์ ๋ณต์ W=[]; %Make 3D coordinate using Direct Linear Transformation %W is wrold coordinate. for i=1:5 A=[ m1(1,i)*P1(3,:) - P1(1,:); m1(2,i)*P1(3,:) - P1(2,:); m2(1,i)*P2(3,:) - P2(1,:); m2(2,i)*P2(3,:) - P2(2,:)]; A(1,:) = A(1,:)/norm(A(1,:)); A(2,:) = A(2,:)/norm(A(2,:)); A(3,:) = A(3,:)/norm(A(3,:)); A(4,:) = A(4,:)/norm(A(4,:)); [u d v] = svd(A); W=[W v(:,4)/v(4,4)]; end %% ๋ค์ 3์ฐจ์ ์ ์์ ํฝ์ ์ ๋ง๋ค๊ธฐ % Now, make image coordiante using P1, P2 from W matrix. reip1 = P1*W; % reip1 is recovered image coordiante reip1 = [reip1(1,:)./reip1(3,:); reip1(2,:)./reip1(3,:)] %3์ฐจ์์์ ๋ณต์๋ ์ด๋ฏธ์ง ์ขํ % m1 is origin image coordinate m1(:,1:5) %์๋ ์ด๋ฏธ์ง ์ขํ reip2 = P2*W; %reip2 is recovered image coordiante reip2 = [reip2(1,:)./reip2(3,:); reip2(2,:)./reip2(3,:)] %3์ฐจ์์์ ๋ณต์๋ ์ด๋ฏธ์ง ์ขํ %m2 is origin image coordinate m2(:,1:5) %์๋ ์ด๋ฏธ์ง ์ขํ-------------------------------------------------------------------------
