#!/usr/bin/perl
#
# perspective_coords [-fx]   x1,y1 .. x4,y4    u1,v1 .. u4,v4
#
# Given two sets of four coordinate pairs, work out 8 constants needed
# to map the first set of coordinates to the second set of coordinates,
# using a perspective transform.
#
# The sixteen numbers given may be comma or space seperated, as one single
# argument or over a number of arguments.
#
# To use the constants generated use the following formulas..
#
#           c1*x + c2*y + c3
#    u  =  ------------------
#           c7*x + c8*y + 1
#
#           c4*x + c5*y + c6
#    v  =  ------------------
#           c7*x + c8*y + 1
#
# By default only the 8 constants needed for the above equation are returned
# (in the order c1 to c8).
#
# However you can also specify an '-fx' option which will output the above
# equations in the form of an ImageMagick "-fx" expression, to lookup colors
# the first 'u' image.
#
####
#
# The above equations can be re-written (by division) into the following form.
#
#   c7*x*u + c8*y*u + u  =  c1*x + c2*y + c3
#   c7*x*v + c8*y*v + v  =  c4*x + c5*y + c6
#
# And simplified (by subtraction) to this form.
#
#    u  =  c1*x + c2*y + c3 - c7*x*u - c8*y*u
#    v  =  c4*x + c5*y + c6 - c7*x*v - c8*y*v
#
# which allows us to generate the matrix of simultaneous equations,
# and thus work out the 8 constants needed.
#
####
#
# Special thanks goes to   Hugemann <Auto@hugemann.de>  for the discussions
# in the IM users mailing list which allowed the formulas to be worked out.
#
# Script by Anthony Thyssen  (December 2006)
#
# Requires the "Math::MatrixReal" perl module to be installed. No compilation
# is needed to build this module.  for the latest version of this module see
# http://leto.net/code/Math-MatrixReal/
#
use strict;
use Math::MatrixReal;
use FindBin;
my $prog = $FindBin::Script;

# Output the program comments as the programs manual
sub Usage {
  print STDERR "$prog: ", @_, "\n";
  @ARGV = ( "$FindBin::Bin/$prog" );
  while( <> ) {
    next if 1 .. 2;
    last if /^###/;
    s/^#$//; s/^# //;
    print STDERR "Usage: " if 3 .. 3;
    print STDERR;
  }
  exit 10;
}

@ARGV = map( /([^\s,]+)/g, @ARGV);   # Spilt arguments by spaces and commas

my $FX = 0;     # output a ImageMagick -fx perspective transformation formula
my $TEST = 0;   # Apply the constants to the input x,y points to test

if ( $ARGV[0] eq '-fx' ) {
  $FX = 1;  shift;
}

if ( @ARGV != 16 ) {
  Usage "Incorrect number of arguments, ", scalar(@ARGV),
        " given\n\t\t\t16 numbers (2 sets of 4 coordinate pairs) are needed.";
}

# Assign the coordinates of the two triangles (remove commas)
my ( $x1, $y1, $x2, $y2, $x3, $y3, $x4, $y4,
     $u1, $v1, $u2, $v2, $u3, $v3, $u4, $v4 ) = @ARGV;

# Convert coordinates into a matrix of simultanious equation
# Such that ...    A * c = r

my $A = Math::MatrixReal->new_from_rows(
  [  [ $x1, $y1, 1.0, 0.0, 0.0, 0.0, -$u1*$x1, -$u1*$y1 ],
     [ 0.0, 0.0, 0.0, $x1, $y1, 1.0, -$v1*$x1, -$v1*$y1 ],
     [ $x2, $y2, 1.0, 0.0, 0.0, 0.0, -$u2*$x2, -$u2*$y2 ],
     [ 0.0, 0.0, 0.0, $x2, $y2, 1.0, -$v2*$x2, -$v2*$y2 ],
     [ $x3, $y3, 1.0, 0.0, 0.0, 0.0, -$u3*$x3, -$u3*$y3 ],
     [ 0.0, 0.0, 0.0, $x3, $y3, 1.0, -$v3*$x3, -$v3*$y3 ],
     [ $x4, $y4, 1.0, 0.0, 0.0, 0.0, -$u4*$x4, -$u4*$y4 ],
     [ 0.0, 0.0, 0.0, $x4, $y4, 1.0, -$v4*$x4, -$v4*$y4 ]
  ] );

my $r = Math::MatrixReal->new_from_cols(
  [ [ $u1, $v1, $u2, $v2, $u3, $v3, $u4, $v4 ] ] );

# Debuging:  output an inverse matrix
#print $A->inverse();

# Solve the matrix for the constants
my ($dim, $c, $base) = $A->decompose_LR()->solve_LR($r);
if ( $dim ) {
  print STDERR
    "$prog: Coordinates given do not form solvable quadrangles.";
  exit 1;
}

# Extract resulting column vector as an array.
my @c = map { $c->element($_,1) }  ( 1 .. 8 );

# Test results with original coordinates
if ( 0 ) {
  print "Test Results against Original coordinates\n";
  printf "%5.1f,%-5.1f -> %5.1f,%-5.1f (should be %5.1f,%-5.1f )\n",
    $x1, $y1,
    ($c[0]*$x1 + $c[1]*$y1 + $c[2]) / ($c[6]*$x1 + $c[7]*$y1 + 1),
    ($c[3]*$x1 + $c[4]*$y1 + $c[5]) / ($c[6]*$x1 + $c[7]*$y1 + 1),
    $u1, $v1;
  printf "%5.1f,%-5.1f -> %5.1f,%-5.1f (should be %5.1f,%-5.1f )\n",
    $x2, $y2,
    ($c[0]*$x2 + $c[1]*$y2 + $c[2]) / ($c[6]*$x2 + $c[7]*$y2 + 1),
    ($c[3]*$x2 + $c[4]*$y2 + $c[5]) / ($c[6]*$x2 + $c[7]*$y2 + 1),
    $u2, $v2;
  printf "%5.1f,%-5.1f -> %5.1f,%-5.1f (should be %5.1f,%-5.1f )\n",
    $x3, $y3,
    ($c[0]*$x3 + $c[1]*$y3 + $c[2]) / ($c[6]*$x3 + $c[7]*$y3 + 1),
    ($c[3]*$x3 + $c[4]*$y3 + $c[5]) / ($c[6]*$x3 + $c[7]*$y3 + 1),
    $u3, $v3;
  printf "%5.1f,%-5.1f -> %5.1f,%-5.1f (should be %5.1f,%-5.1f )\n",
    $x4, $y4,
    ($c[0]*$x4 + $c[1]*$y4 + $c[2]) / ($c[6]*$x4 + $c[7]*$y4 + 1),
    ($c[3]*$x4 + $c[4]*$y4 + $c[5]) / ($c[6]*$x4 + $c[7]*$y4 + 1),
    $u4, $v4;
  print "\n";
}

if ( ! $FX ) {
  # Output the constants
  map { printf "%9.6f\n", $_ }  @c;
}
else {
  # Output a IM -fx expression
  printf "det=%+.6f*i %+.6f*j +1;\n".
         "xx=(%+.6f*i %+.6f*j %+.6f)/det;\n".
         "yy=(%+.6f*i %+.6f*j %+.6f)/det;\n",
         @c[6,7,0..5];
}