#!perl -l

use 5.012; # //, strict, say
use warnings;

use PDL 2.057;
#use PDL::Matrix;

print "perl version $]";
print "PDL version $PDL::VERSION";

print my $red_in  = (pdl [ # reference
    [419 , 17 ], #  p01
    [772 , 13 ], #  p02
    [348 , 102], #   p03
    [720 , 95 ], #  p04
    [58  , 238], #   p05
    [218 , 236], #   p06
    [55  , 333], #   p07
    [58  , 391], #   p08
    [60  , 498], #   p09
    [58  , 554], #   p10
    [1061, 553], #   p11
    [60  , 892], #   p12
])->append(ones(1))->transpose;
print my $red_out = (pdl [ # target_41
    [487, 466], # p01
    [531, 630], # p02
    [439, 447], # p03
    [486, 618], # p04
    [340, 328], # p05
    [358, 401], # p06
    [297, 343], # p07
    [274, 349], # p08
    [227, 366], # p09
    [205, 374], # p10
    [313, 809], # p11
    [77 , 420], # p12
])->append(ones(1))->transpose;

print my $blue_in = (pdl [ # reference
    [438,331],  # q01
    [432,392],  # q02
    [436,493],  # q03
    [436,552],  # q04
    [439,602],  # q05
])->append(ones(1))->transpose;
print my $blue_out_for_checking = (pdl [ # target_41
    [339,508],  # q01
    [314,518],  # q02
    [265,533],  # q03
    [245,537],  # q04
    [228,544],  # q05
])->append(ones(1))->transpose;

$\ = '';
print for
    "sx  = ", my $sx  = sum( $red_in->transpose->slice(0) ),                    "\n",
    "sy  = ", my $sy  = sum( $red_in->transpose->slice(1) ),                    "\n",
    "s1  = ", my $s1  = sum( $red_in->transpose->slice(2) ),                    "\n",
    "sxx = ", my $sxx = sum( $red_in->transpose->slice(0)**2 ),                 "\n",
    "syy = ", my $syy = sum( $red_in->transpose->slice(1)**2 ),                 "\n",
    "sxy = ", my $sxy = sum( $red_in->transpose->slice(0)*$red_in->transpose->slice(1) ),  "\n",
    "sv  = ", my $sv  = sum( $red_out->transpose->slice(0) ),                   "\n",
    "sw  = ", my $sw  = sum( $red_out->transpose->slice(1) ),                   "\n",
    "sxv = ", my $sxv = sum( $red_in->transpose->slice(0)*$red_out->transpose->slice(0) ), "\n",
    "sxw = ", my $sxw = sum( $red_in->transpose->slice(0)*$red_out->transpose->slice(1) ), "\n",
    "syv = ", my $syv = sum( $red_in->transpose->slice(1)*$red_out->transpose->slice(0) ), "\n",
    "syw = ", my $syw = sum( $red_in->transpose->slice(1)*$red_out->transpose->slice(1) ), "\n",
    ;

$\ = "\n";

print "G => ", my $G = pdl
[
    [ $sxx, $sxy, $sx ],
    [ $sxy, $syy, $sy ],
    [ $sx,  $sy,  $s1 ],
];

print "Q => ", my $Q = pdl
[
    [ $sxv, $sxw ],
    [ $syv, $syw ],
    [ $sv,  $sw  ],
];


print "P => ", my $P = inv($G) x $Q;

# print "P->transpose", $P->transpose;
# print "P->append->transpose", $P->append(zeros(1))->transpose;

print "abcdef001 => ", my $abcdef001 = $P->append(zeros(1))->transpose->set(2,2,1);

print "blue out calc => ", my $blue_out_calc = $abcdef001 x $blue_in;

print "error => ", $blue_out_calc - $blue_out_for_checking; # no more than one pixel off

__END__

perl version 5.032001
PDL version 2.080

[
 [ 419  772  348  720   58  218   55   58   60   58 1061   60]
 [  17   13  102   95  238  236  333  391  498  554  553  892]
 [   1    1    1    1    1    1    1    1    1    1    1    1]
]


[
 [487 531 439 486 340 358 297 274 227 205 313  77]
 [466 630 447 618 328 401 343 349 366 374 809 420]
 [  1   1   1   1   1   1   1   1   1   1   1   1]
]


[
 [438 432 436 436 439]
 [331 392 493 552 602]
 [  1   1   1   1   1]
]


[
 [339 314 265 245 228]
 [508 518 533 537 544]
 [  1   1   1   1   1]
]

sx  = 3887
sy  = 3922
s1  = 12
sxx = 2604611
syy = 2052390
sxy = 929565
sv  = 4034
sw  = 5551
sxv = 1608891
sxw = 2354880
syv = 945962
syw = 1755275
G =>
[
 [2604611  929565    3887]
 [ 929565 2052390    3922]
 [   3887    3922      12]
]

Q =>
[
 [1608891 2354880]
 [ 945962 1755275]
 [   4034    5551]
]

P =>
[
 [0.11504241 0.44420653]
 [-0.4325134 0.11994661]
 [ 440.26231  279.49489]
]

abcdef001 =>
[
 [0.11504241 -0.4325134  440.26231]
 [0.44420653 0.11994661  279.49489]
 [         0          0          1]
]

blue out calc =>
[
 [ 347.48895  320.41538  277.19169   251.6734  230.39286]
 [ 513.75967  518.41118  532.30261  539.37946  546.70941]
 [         1          1          1          1          1]
]

error =>
[
 [   8.488949   6.4153774   12.191694   6.6734037   2.3928612]
 [  5.7596723  0.41117608 -0.69739064   2.3794591    2.709409]
 [          0           0           0           0           0]
]