% limpa
clear

% carrega e salva a imagem em uma variável
imagem = imread('test3.png');
imshow(imagem);

% recebe dados do usuário
[x1,y1] = ginput(2); % L1: primeira reta do primeiro par de retas paralelas
line(x1, y1, 'color', [.8, .0, .0], 'linewidth', 3); % desenha linha
[x2,y2] = ginput(2); % L2: segunda reta do primeiro par de retas paralelas
line(x2, y2, 'color', [.8, .0, .0], 'linewidth', 3); % desenha linha
[x3,y3] = ginput(2); % L3: primeira reta do segundo par de retas paralelas
line(x3, y3, 'color', [.0, .8, .0], 'linewidth', 3); % desenha linha
[x4,y4] = ginput(2); % L4: segunda reta do segundo par de retas paralelas
line(x4, y4, 'color', [.0, .8, .0], 'linewidth', 3); % desenha linha
[x5,y5] = ginput(2);  % L5: primeira reta do primeiro par de retas ortogonais
line(x5, y5, 'color', [.0, .0, .8], 'linewidth', 3); % desenha linha
[x6,y6] = ginput(2);  % L6: segunda reta do primeiro par de retas ortogonais
line(x6, y6, 'color', [.0, .0, .8], 'linewidth', 3); % desenha linha
[x7,y7] = ginput(2);  % L7: primeira reta do segundo par de retas ortogonais
line(x7, y7, 'color', [.0, .8, .8], 'linewidth', 3); % desenha linha
[x8,y8] = ginput(2);  % L8: segunda reta do segundo par de retas ortogonais
line(x8, y8, 'color', [.0, .8, .8], 'linewidth', 3); % desenha linha


% acha equações das linhas
A = [x1(1) y1(1) 1];
B = [x1(2) y1(2) 1];
L1 = cross(A,B);
L1 = L1./norm(L1);

A = [x2(1) y2(1) 1];
B = [x2(2) y2(2) 1];
L2 = cross(A,B);
L2 = L2./norm(L2);

A = [x3(1) y3(1) 1];
B = [x3(2) y3(2) 1];
L3 = cross(A,B);
L3 = L3./norm(L3);


A = [x4(1) y4(1) 1];
B = [x4(2) y4(2) 1];
L4 = cross(A,B);
L4 = L4./norm(L4);

A = [x5(1) y5(1) 1];
B = [x5(2) y5(2) 1];
L5 = cross(A,B);
L5 = L5./norm(L5);

A = [x6(1) y6(1) 1];
B = [x6(2) y6(2) 1];
L6 = cross(A,B);
L6 = L6./norm(L6);

A = [x7(1) y7(1) 1];
B = [x7(2) y7(2) 1];
L7 = cross(A,B);
L7 = L7./norm(L7);

A = [x8(1) y8(1) 1];
B = [x8(2) y8(2) 1];
L8 = cross(A,B);
L8 = L8./norm(L8);

% Mapeia os pontos de intersecção das retas paralelas
X1 = cross(L1, L2);
X1 = X1./norm(X1);
X2 = cross(L3, L4);
X2 = X2./norm(X2);

% Encontra a linha do infinito
LInf = cross(X1, X2);

% Faz a retificação afim através da linha do infinito
H = [1 0 0; 0 1 0; LInf];
transformação = projective2d(H');
out = imwarp(imagem, transformação);
imshow(out);

% Mapeia as linhas ortogonais para as respectivas novas linhas na imagem retificada
L5eq = H'\L5';
L5 = L5eq./norm(L5eq);
L6eq = H'\L6';
L6 = L6eq./norm(L6eq);
L7eq = H'\L7';
L7 = L7eq./norm(L7eq);
L8eq = H'\L8';
L8 = L8eq./norm(L8eq);