A Linear Method to Derive 3D Projective Invariants from 4 Uncalibrated Images

Algorithm 1

X = RandomReal[{−1000, 1000}, {6, 4}];

M = RandomReal[{−1, 1}, {4, 3, 4}];

T = RandomReal[{0, 1}, {4, 6}];

x = RandomReal[{0, 1}, {4, 6, 3}];

a = RandomReal[{0, 1}, {4, 6, 4}];

b = RandomReal[{0, 1}, {4, 6, 4}];

u = RandomReal[{0, 1}, {4, 6}];

v = RandomReal[{0, 1}, {4, 6}];

X[[1,4]] = 1;

X[[2,4]] = 1;

X[[3,4]] = 1;

X[[4,4]] = 1;

X[[5,4]] = 1;

X[[6,4]] = 1;

XT = Transpose[{X[[1]], X[[2]], X[[3]], X[[4]]}];

A = LinearSolve[XT, X[[5]]];

B = LinearSolve[XT, X[[6]]];

Inv1 = (A[[1]] B[[2]])/(A[[2]] B[[1]]);

Inv2 = (A[[1]] B[[3]])/(A[[3]] B[[1]]);

Inv3 = (A[[1]] B[[4]])/(A[[4]] B[[1]]);

Print[“The three invariants computed from 3D point

locations: ”, Inv1, “ ”, Inv2, “ ”, Inv3];

For[i = 1, i <= 4, i++, For[j = 1, j <= 6, j++,

x[[i, j]] = M[[i]] ⋅ X[[j]];

u[[i, j]] = x[[i, j, 1]]/x[[i, j, 3]];

v[[i, j]] = x[[i, j, 2]]/x[[i, j, 3]];

]];

For[i = 1, i <= 4, i++, For[j = 5, j <= 6, j++,

a[[i, j, 1]] = u[[i, 1]] − u[[i, j]];

a[[i, j, 2]] = u[[i, 2]] − u[[i, j]];

a[[i, j, 3]] = u[[i, 3]] − u[[i, j]];

a[[i, j, 4]] = u[[i, 4]] − u[[i, j]];

b[[i, j, 1]] = v[[i, 1]] − v[[i, j]];

b[[i, j, 2]] = v[[i, 2]] − v[[i, j]];

b[[i, j, 3]] = v[[i, 3]] − v[[i, j]];

b[[i, j, 4]] = v[[i, 4]] − v[[i, j]];

]];

For[i = 1, i <= 4, i++,

T[[i, 1]] = (a[[i, 5, 3]] b[[i, 5, 4]] − a[[i, 5, 4]] b[[i, 5, 3]])

   (a[[i, 6, 1]] b[[i, 6, 2]] − a[[i, 6, 2]] b[[i, 6, 1]]);

T[[i, 2]] = (a[[i, 5, 4]] b[[i, 5, 2]] − a[[i, 5, 2]] b[[i, 5, 4]])

   (a[[i, 6, 1]] b[[i, 6, 3]] − a[[i, 6, 3]] b[[i, 6, 1]]);

T[[i, 3]] = (a[[i, 5, 2]] b[[i, 5, 3]] − a[[i, 5, 3]] b[[i, 5, 2]])

   (a[[i, 6, 1]] b[[i, 6, 4]] − a[[i, 6, 4]] b[[i, 6, 1]]);

T[[i, 4]] = (a[[i, 5, 1]] b[[i, 5, 4]] − a[[i, 5, 4]] b[[i, 5, 1]])

   (a[[i, 6, 2]] b[[i, 6, 3]] − a[[i, 6, 3]] b[[i, 6, 2]]);

T[[i, 5]] = (a[[i, 5, 3]] b[[i, 5, 1]] − a[[i, 5, 1]] b[[i, 5, 3]])

   (a[[i, 6, 2]] b[[i, 6, 4]] − a[[i, 6, 4]] b[[i, 6, 2]]);

T[[i, 6]] = (a[[i, 5, 1]] b[[i, 5, 2]] − a[[i, 5, 2]] b[[i, 5, 1]])

   (a[[i, 6, 3]] b[[i, 6, 4]] − a[[i, 6, 4]] b[[i, 6, 3]]);

];

I1 = −Det[{

{T[[1,1]], T[[1,4]] + T[[1,5]], T[[1,3]], T[[1,6]]},

{T[[2,1]], T[[2,4]] + T[[2,5]], T[[2,3]], T[[2,6]]},

{T[[3,1]], T[[3,4]] + T[[3,5]], T[[3,3]], T[[3,6]]},

{T[[4,1]], T[[4,4]] + T[[4,5]], T[[4,3]], T[[4,6]]}

}]/Det[{

{T[[1,1]], T[[1,4]], T[[1,5]], T[[1,6]]},

{T[[2,1]], T[[2,4]], T[[2,5]], T[[2,6]]},

{T[[3,1]], T[[3,4]], T[[3,5]], T[[3,6]]},

{T[[4,1]], T[[4,4]], T[[4,5]], T[[4,6]]}

}];

I2 = −Det[{

{T[[1,4]] + T[[1,6]], T[[1,2]], T[[1,3]], T[[1,5]]},

{T[[2,4]] + T[[2,6]], T[[2,2]], T[[2,3]], T[[2,5]]},

{T[[3,4]] + T[[3,6]], T[[3,2]], T[[3,3]], T[[3,5]]},

{T[[4,4]] + T[[4,6]], T[[4,2]], T[[4,3]], T[[4,5]]}

}]/Det[{

{T[[1,4]], T[[1,2]], T[[1,6]], T[[1,5]]},

{T[[2,4]], T[[2,2]], T[[2,6]], T[[2,5]]},

{T[[3,4]], T[[3,2]], T[[3,6]], T[[3,5]]},

{T[[4,4]], T[[4,2]], T[[4,6]], T[[4,5]]}

}];

I3 = −Det[{

{T[[1,5]] + T[[1,6]], T[[1,2]], T[[1,3]], T[[1,4]]},

{T[[2,5]] + T[[2,6]], T[[2,2]], T[[2,3]], T[[2,4]]},

{T[[3,5]] + T[[3,6]], T[[3,2]], T[[3,3]], T[[3,4]]},

{T[[4,5]] + T[[4,6]], T[[4,2]], T[[4,3]], T[[4,4]]}

}]/Det[{

{T[[1,5]], T[[1,6]], T[[1,3]], T[[1,4]]},

{T[[2,5]], T[[2,6]], T[[2,3]], T[[2,4]]},

{T[[3,5]], T[[3,6]], T[[3,3]], T[[3,4]]},

{T[[4,5]], T[[4,6]], T[[4,3]], T[[4,4]]}

}];

Print[“The three invariants computed from 2D projections: ”,

I1, “ ”, I2, “ ”, I3];