[scilab-Users] Matrix transformations

[Samuel Gougeon]

Yes Denis, i agree that A\B is meaningless, since the 4x4 result can be 
applied
to a set of four and only four points and that it has no geometrical 
meaning.

To stick with geometry, a rotation o translation o scaling  transform in 
2D is linear and affine
and then can be described as P2 = A.P1 + B
where A is a 2x2 constant matrix , B is a 2x1 constant column,
P1 is the column of coordinates of a given point
P2 is the column of coordinates of the transformed point
If one wants to transform a whole set of m points,
P1 & P2 should be 2xm matrices (each column is a point in P1 and its 
transformed in P2).

Then, A have 4 coefficients and B have 2 ; so 6 coefficients have to be 
determined.
To do that, exactly 3 independant points (=6 coeffs) and their transform 
give 6
equations and are needed and sufficient.

Now, Capitao Obvio is giving 4 points. With all of these, the problem is 
overdetermined (8 equations).
This leads to a linear* fitting* problem, that may be processed using 
*reglin()*:

Renaming his A => P1  and his B => P2, the best 2x2 A matrix and 2x1 B 
column
such that P2 = A.P1 + B -- according to a least square criterium -- are 
given by
[ A, B ] = reglin(P1,P2);
That gives:
P1 = [
52.363965, 4.892435 ;
52.384611, 4.898272 ;
52.379267, 4.881105 ;
52.369678, 4.911489
]
P2 = [
0, 0 ;
2131, 1534 ;
2131, 0 ;
0, 1534
]

-->[A, B] = reglin(P1.',P2.')

B  =
   - 5747613.8
   - 3211459.3
  A  =
     112982.37  - 34458.388
     55002.261    67714.214

-->A*P1.' + B*ones(P1(:,1).')  // Checking:
  ans  =

     5.5829359    2137.0833    2124.8527  - 5.5189172
   - 35.430415    1495.3941    39.012139    1569.0241

-->P2.'    // To be compared to the previous fit
  ans  =
     0.    2131.    2131.    0.
     0.    1534.    0.       1534.

-->// For an exact fit, only 3 points must be processed, for example the 
3 first ones:
-->P1 = P1(1:3,:), P2 = P2(1:3,:)
  P1  =
     52.363965    4.892435
     52.384611    4.898272
     52.379267    4.881105
  P2  =
     0.       0.
     2131.    1534.
     2131.    0.

-->[A, B] = reglin(P1(1:3,:).',P2.')
  B  =
   - 5754011.4
   - 3170858.6
  A  =
     113176.65  - 35231.318
     53769.285    72619.382

-->A*P1.' + B*ones(P1(:,1).')  // Checking:
  ans  =
     9.313D-10    2131.    2131.
     4.657D-10    1534.    4.657D-10

-->P2.'
  ans  =
     0.    2131.    2131.
     0.    1534.    0.

This might help more than A\B ;)
Regards

Samuel


Le 14/10/2011 10:33, CRETE Denis a écrit :
>
> Hello,
>
> I do not agree with the approach given below, as the transform acts on 
> the plane (RxR) to another plane, i.e. it acts from R² to R². Hence, 
> the transform, if linear, would be described by a 2x2 matrix (and not 
> 4x4).
>
> It is very plausible that the initial difficulty (finding a linear 
> transform such as rotation o translation o scaling) comes from a 
> distortion of the "picture" on which a "path" should be mapped (these 
> terms are generic). So, the transform is not linear; it can be 
> represented by the system:
>
> X=g_x(x,y)
>
> Y=g_y(x,y)
>
> where x and y are given in the ma trix A and X,Y in B.
>
> The function g_x and g_y can be expanded to the second order (for the 
> simplest description of non-linearity)
>
> X=a_x+b_x*x+c_x*y+d_x*x²+e_x*x*y+f_x*y²
>
> Y=a_y+b_y*x+c_y*y+d_y*x²+e_y*x*y+f_y*y²
>
> The 12 coefficients a_x, a_y, b_x, b_y... f_x and f_y are determined 
> with the points contained in A and B. However, the number of points 
> must be increased to 6 in order to get 12 equations (giving 
> [X_1;Y_1;X_2;Y_2...X_6 ;Y_6] as a product of a 12x12 matrix C and the 
> vector [a_x; a_y;b_x;...f_y]) which can be inverted to get the vector 
> [a_x, a_y,b_x,...f_y] as a product of the 12x12 matrix inv(C) and the 
> vector [X_1;Y_1;X_2;Y_2...X_6;Y_6].
>
> But with usual optical distortion, it is likely that the second order 
> terms are zero (when the origin is chosen at the center of the image) 
> and the 3rd order terms should be considered. One can try the following
>
> X=a_x + b_x*x + c_x*y + d_x*x*y^2
>
> Y=a_y + b_y*x + c_y*y + d_y*x^2*y
>
> with proper choice of the origin (center of the image). Only the 
> matrix C needs to be changed (each pair of row is 
> [1,x_i,y_i,x_i*y_i^2; 1,x_i,y_i,x_i^2*y_i] for i=1...4) , otherwise 
> the procedure is the same as for the 2nd order case presented above. I 
> have already used this kind of transform to correct image distortion. 
> In this case, careful selection of the points is necessary: I fear 
> that choosing the 4 corners of the image is not appropriate. Choosing 
> the center, one corner and 2 "middle-side" points  of the image might 
> give better results.
>
> HTH
>
> Denis
>

Le 12/10/2011 18:25, Capitao Obvio a écrit :
> Hi!
>
> I'm starting to work with Scilab, and I'm trying to figure out how to 
> work with matrices.
>
> My goal is to find a way to convert any point from one system into 
> another.
>
> My (failed) approach so far was to estimate the necessary rotation, 
> translation and scale; but it didn't work quite well, because the 
> second map is slightly distorted (probably due perspective).
>
> Being a bit more specific:
>
> I have the following points:
>
> A = [
> [52.363965, 4.892435],
> [52.384611, 4.898272],
> [52.379267, 4.881105],
> [52.369678, 4.911489],
> ]
>
>
> They correspond to these points:
>
> B = [
> [0, 0],
> [2131, 1534],
> [2131, 0],
> [0, 1534],
> ]
>
> Do you have any suggestions on how to calculate tje transform from A to B?
>
> Thanks in advance,
>
> Cap.
>

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://lists.scilab.org/pipermail/users/attachments/20111017/f757067c/attachment.htm>