#!perl -l

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

use PDL;
use PDL::Matrix;

print my $red_in  = (mpdl [[58,48], [108,155], [186,80], [255,191], [331,48]])->transpose->append(ones(1));
print my $red_out = (mpdl [[471,15], [531,141], [603,90], [682,227], [747,107]])->transpose->append(ones(1));

print my $blue_in = (mpdl [[125,73], [197,158], [282,94]])->transpose->append(ones(1));
print my $blue_out_for_checking = (mpdl [[542,64], [622,175], [701,138]])->transpose->append(ones(1));

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

print my $G = mpdl
[
    [ $sxx, $sxy, $sx ],
    [ $sxy, $syy, $sy ],
    [ $sx,  $sy,  $s1 ],
];

print my $Q = mpdl
[
    [ $sxv, $sxw ],
    [ $syv, $syw ],
    [ $sv,  $sw  ],
];

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

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

print my $blue_out_calc = $abcdef001 x $blue_in;

print $blue_out_calc - $blue_out_for_checking; # no more than one pixel off