load("qt","lapack","ginac");
C := [[[213,366],[0,0,0]],[[394,330],[30,0,0]],[[559,299],[60,0,0]],
     [[201,239],[0,30,0]],[[348,216],[30,30,0]],[[514,72],[60,30,30]],
     [[194,155],[0,60,0]],[[317,23],[30,60,30]],[[431,122],[60,60,0]]]
[[[213,366],[0,0,0]],[[394,330],[30,0,0]],[[559,299],[60,0,0]],[[201,239],[0,30,0]],[[348,216],[30,30,0]],[[514,72],[60,30,30]],[[194,155],[0,60,0]],[[317,23],[30,60,30]],[[431,122],[60,60,0]]]
A := matrix(18,11)
b := matrix(18,1)
for local i := 1 as i <= C:size() do begin
	local x := (C[i])[2]
	local p := (C[i])[1]
	A[2*i-1,1] := x[1]
	A[2*i-1,2] := x[2]
	A[2*i-1,3] := x[3]
	A[2*i-1,4] := 1
	A[2*i-1,9] := -p[1]*x[1]
	A[2*i-1,10] := -p[1]*x[2]
	A[2*i-1,11] := -p[1]*x[3]
	b[2*i-1] := p[1]
	A[2*i,5] := x[1]
	A[2*i,6] := x[2]
	A[2*i,7] := x[3]
	A[2*i,8] := 1.0
	A[2*i,9] := -p[2]*x[1]
	A[2*i,10] := -p[2]*x[2]
	A[2*i,11] := -p[2]*x[3]
	b[2*i] := p[2]
end;
[A,b]
⎡⎛ 0   0   0   1  0   0   0    0     0       0       0    ⎞ ⎛ 213 ⎞⎤
⎢⎜                                                        ⎟ ⎜     ⎟⎥
⎢⎜ 0   0   0   0  0   0   0   1.0    0       0       0    ⎟ ⎜ 366 ⎟⎥
⎢⎜                                                        ⎟ ⎜     ⎟⎥
⎢⎜ 30  0   0   1  0   0   0    0   -11820    0       0    ⎟ ⎜ 394 ⎟⎥
⎢⎜                                                        ⎟ ⎜     ⎟⎥
⎢⎜ 0   0   0   0  30  0   0   1.0  -9900     0       0    ⎟ ⎜ 330 ⎟⎥
⎢⎜                                                        ⎟ ⎜     ⎟⎥
⎢⎜ 60  0   0   1  0   0   0    0   -33540    0       0    ⎟ ⎜ 559 ⎟⎥
⎢⎜                                                        ⎟ ⎜     ⎟⎥
⎢⎜ 0   0   0   0  60  0   0   1.0  -17940    0       0    ⎟ ⎜ 299 ⎟⎥
⎢⎜                                                        ⎟ ⎜     ⎟⎥
⎢⎜ 0   30  0   1  0   0   0    0     0     -6030     0    ⎟ ⎜ 201 ⎟⎥
⎢⎜                                                        ⎟ ⎜     ⎟⎥
⎢⎜ 0   0   0   0  0   30  0   1.0    0     -7170     0    ⎟ ⎜ 239 ⎟⎥
⎢⎜                                                        ⎟ ⎜     ⎟⎥
⎢⎜ 30  30  0   1  0   0   0    0   -10440  -10440    0    ⎟ ⎜ 348 ⎟⎥
⎢⎜                                                        ⎟,⎜     ⎟⎥
⎢⎜ 0   0   0   0  30  30  0   1.0  -6480   -6480     0    ⎟ ⎜ 216 ⎟⎥
⎢⎜                                                        ⎟ ⎜     ⎟⎥
⎢⎜ 60  30  30  1  0   0   0    0   -30840  -15420  -15420 ⎟ ⎜ 514 ⎟⎥
⎢⎜                                                        ⎟ ⎜     ⎟⎥
⎢⎜ 0   0   0   0  60  30  30  1.0  -4320   -2160   -2160  ⎟ ⎜ 72  ⎟⎥
⎢⎜                                                        ⎟ ⎜     ⎟⎥
⎢⎜ 0   60  0   1  0   0   0    0     0     -11640    0    ⎟ ⎜ 194 ⎟⎥
⎢⎜                                                        ⎟ ⎜     ⎟⎥
⎢⎜ 0   0   0   0  0   60  0   1.0    0     -9300     0    ⎟ ⎜ 155 ⎟⎥
⎢⎜                                                        ⎟ ⎜     ⎟⎥
⎢⎜ 30  60  30  1  0   0   0    0   -9510   -19020  -9510  ⎟ ⎜ 317 ⎟⎥
⎢⎜                                                        ⎟ ⎜     ⎟⎥
⎢⎜ 0   0   0   0  30  60  30  1.0   -690   -1380    -690  ⎟ ⎜ 23  ⎟⎥
⎢⎜                                                        ⎟ ⎜     ⎟⎥
⎢⎜ 60  60  0   1  0   0   0    0   -25860  -25860    0    ⎟ ⎜ 431 ⎟⎥
⎢⎜                                                        ⎟ ⎜     ⎟⎥
⎣⎝ 0   0   0   0  60  60  0   1.0  -7320   -7320     0    ⎠ ⎝ 122 ⎠⎦
v := (A^t*A)^(-1)*A^t*b
⎛   6.729327652246897333    ⎞
⎜                           ⎟
⎜   1.3066454617655438685   ⎟
⎜                           ⎟
⎜   -2.238960277276364933   ⎟
⎜                           ⎟
⎜   213.22004126984835398   ⎟
⎜                           ⎟
⎜  -0.5933225103918719457   ⎟
⎜                           ⎟
⎜  -2.2227873537201961302   ⎟
⎜                           ⎟
⎜  -6.0882261648313575936   ⎟
⎜                           ⎟
⎜   366.08058086446076315   ⎟
⎜                           ⎟
⎜ 0.0017555467142436832819  ⎟
⎜                           ⎟
⎜  0.008442792379581966193  ⎟
⎜                           ⎟
⎝ -0.0071105508747814286606 ⎠
K := matrix([
	[v[1],v[2],v[3],v[4]],
	[v[5],v[6],v[7],v[8]],
	[v[9],v[10],v[11],1.0]])
⎛   6.729327652246897333     1.3066454617655438685     -2.238960277276364933    213.22004126984835398 ⎞
⎜                                                                                                     ⎟
⎜  -0.5933225103918719457   -2.2227873537201961302    -6.0882261648313575936    366.08058086446076315 ⎟
⎜                                                                                                     ⎟
⎝ 0.0017555467142436832819  0.008442792379581966193  -0.0071105508747814286606           1.0          ⎠
u := K*vector([30,30,0,1])
⎛ 454.29923469022159002 ⎞
⎜                       ⎟
⎜ 281.59728494109872088 ⎟
⎜                       ⎟
⎝ 1.3059501728147694843 ⎠
u/u[3]
⎛ 347.86873507666168456 ⎞
⎜                       ⎟
⎜ 215.62636217135315096 ⎟
⎜                       ⎟
⎝          1.0          ⎠
M := lmatrix_t([
	[double(v[1]),double(v[2]),double(v[3])],
	[double(v[5]),double(v[6]),double(v[7])],
	[double(v[9]),double(v[10]),double(v[11])]])
⎛  6.72933     1.30665     -2.23896   ⎞
⎜ -0.593323    -2.22279    -6.08823   ⎟
⎝ 0.00175555  0.00844279  -0.00711055 ⎠
D := qrdecomposition(M^(-1))
⎡⎛ 0.976146   -0.149891  0.157069  ⎞ ⎛ 0.15818  0.00230436  -49.5192 ⎞⎤
⎢⎜ -0.216517  -0.618483   0.75538  ⎟,⎜    0      0.162339   -30.5158 ⎟⎥
⎣⎝ -0.01608   -0.77137   -0.636184 ⎠ ⎝    0         0       89.4704  ⎠⎦
M^(-1)
⎛  0.154407    -0.0220837  -29.7109 ⎞
⎜ -0.0342487   -0.100903   97.1794  ⎟
⎝ -0.00254353   -0.12526   -32.5844 ⎠
D[1]*D[2]
⎛  0.154407    -0.0220837  -29.7109 ⎞
⎜ -0.0342487   -0.100903   97.1794  ⎟
⎝ -0.00254353   -0.12526   -32.5844 ⎠
W := D[2]^(-1)
⎛ 6.3219  -0.089738   3.46838  ⎞
⎜   0      6.15996    2.10099  ⎟
⎝   0         0      0.0111769 ⎠
T := lmatrix_t([[1,0,0],[0,1,0],[0,0,1]])
⎛ 1  0  0 ⎞
⎜ 0  1  0 ⎟
⎝ 0  0  1 ⎠
lambda := 1.0/W[3,3]
89.47041585008822095
P := W*T*lambda
⎛ 565.623  -8.0289  310.317 ⎞
⎜    0     551.135  187.976 ⎟
⎝    0        0        1    ⎠
R := T*D[1]^(-1)
⎛ 0.976146   -0.216517  -0.01608  ⎞
⎜ -0.149891  -0.618483  -0.77137  ⎟
⎝ 0.157069    0.75538   -0.636184 ⎠
P*R/lambda
⎛  6.72933     1.30665     -2.23896   ⎞
⎜ -0.593323    -2.22279    -6.08823   ⎟
⎝ 0.00175555  0.00844279  -0.00711055 ⎠
tt := -(P*R/lambda)^(-1)*lvector_t([double(v[4]),double(v[8]),1.0])
⎛ 4.87261  ⎞
⎜ -52.9384 ⎟
⎝  78.982  ⎠
yy := (P*R)^(-1)*lvector_t([320,240,1])
⎛ 0.160939  ⎞
⎜ 0.693003  ⎟
⎝ -0.709293 ⎠
beta := -tt[3,1]/yy[3,1]
111.353152369761545515
yy*beta+tt
⎛ 22.7937 ⎞
⎜ 24.2296 ⎟
⎝    0    ⎠
image := QImage("platte.jpg")
 
result := QImage(200,200)
for i := 0 as i < 200 do
	for j := 0 as j < 200 do begin
		local b := K*vector([0.3*i,0.3*j,0,1])
		local x := b[1]/b[3]
		local y := b[2]/b[3]
		result[i,j] := image[x,y]
  end
[26,21,17]
result
 
result:save("rectified.jpg")
K
⎛   6.729327652246897333     1.3066454617655438685     -2.238960277276364933    213.22004126984835398 ⎞
⎜                                                                                                     ⎟
⎜  -0.5933225103918719457   -2.2227873537201961302    -6.0882261648313575936    366.08058086446076315 ⎟
⎜                                                                                                     ⎟
⎝ 0.0017555467142436832819  0.008442792379581966193  -0.0071105508747814286606           1.0          ⎠
v := vector([1,2,3,1])
⎛ 1 ⎞
⎜   ⎟
⎜ 2 ⎟
⎜   ⎟
⎜ 3 ⎟
⎜   ⎟
⎝ 1 ⎠
w := K*v
⎛ 215.84577901379724425 ⎞
⎜                       ⎟
⎜ 342.77700515213442617 ⎟
⎜                       ⎟
⎝ 0.9973094788490633297 ⎠
w/w[3]
⎛ 216.42808334970631369  ⎞
⎜                        ⎟
⎜ 343.70174195848751444  ⎟
⎜                        ⎟
⎝ 0.99999999999999999995 ⎠