#!/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 separated, 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 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]; }