% parameters of projective (homographic) coordinate transformation
% between six 1-D ``pictures'' (five successive pairs sequence)
% just made up these numbers (e.g. from 2-D sciam sequence, took x-parameters
p=[...
   0.9985  -0.0135 -0.0008;...
   0.9996  -0.0158 -0.0005;...
   1.0033  -0.0162  0.0004;...
   1.0045  -0.0203  0.0015;...
   1.0055  -0.0256  0.0016...
  ];

% a little experiment: make a new p, based on a particular K:
% K=[2 -.5;0,1]
% % 10 degree pan...   do various pan angles:
% for thetaperdeg=-10:10
% pchirp(rand(1,10),thetaperdeg*o,K)
% end
% resulting in family of p values...
p=[...
0.9155588      -0.3588750       0.0844412;...
0.9155588      -0.3588750       0.0844412;...
1.0922291       0.3919737      -0.0922291;...
1.0824572       0.3504431      -0.0824572;...
1.0728293       0.3095245      -0.0728293;...
1.0633365       0.2691800      -0.0633365;...
1.0539702       0.2293735      -0.0539702;...
1.0447225       0.1900707      -0.0447225;...
1.0355855       0.1512384      -0.0355855;...
1.0265518       0.1128450      -0.0265518;...
1.0176142       0.0748602      -0.0176142;...
1.0087658       0.0372546      -0.0087658;...
1.0000000       0.0000000       0.0000000;...
0.9913104      -0.0369309       0.0086896;...
0.9826907      -0.0735644       0.0173093;...
0.9741350      -0.1099263       0.0258650;...
0.9656373      -0.1460414       0.0343627;...
0.9571920      -0.1819341       0.0428080;...
0.9487934      -0.2176281       0.0512066;...
0.9404361      -0.2531466       0.0595639;...
0.9321147      -0.2885123       0.0678853;...
0.9238241      -0.3237478       0.0761759;...
0.9155588      -0.3588750       0.0844412;...
  ];

X=[];   % initialize the ``X'' of Hartley, page 8
for m=1:5
  p11=p(m,1); p12=p(m,2);
  p21=p(m,3); p22=1;

  q=inv([p11,p12;p21,p22]);
  q=q/q(2,2);

  q11=q(1); q12=q(2);
  q21=q(3); q22=1;

  Xm=[q11-p11    q21-p12        0      ;...
      q12        q22-p11       -p12    ;...
...
         -p21    q11-p22    q21        ;...
      0          q12-p21    q22-p22    ;...
     ];
  X=[X;Xm];

end%for

[v d]=eig(X'*X);
d=diag(d);
[trash location]=min(d);  % choose least eigenvalue
answer=v(:,location);

C=[answer(1) answer(2);...
   answer(2) answer(3)...
  ];
K=chol(C);

disp('still need to see why K(2,2) not equal to 1')

K=K/K(2,2);

% final result:
%K =
%    2.0471   -0.4579
%         0    1.0000
% is close to what it should be

