use strict;
use Image::Magick;
use Math::Matrix;
use Math::Gradient qw(gradient);

my $rate        = 500;
my $tension	= 0.1;
my(@coords)	= (
                   [ 23, 24],     [123, 64],     [167,200],     [212,285],     [293,297],
                   [552,205],     [537, 71],     [358,122],     [262,130],     [238, 24]
                  );

my $image	= Image::Magick->new;           # Create new image
$image		->Set(size=>'720x480');         # Set size
$image		->ReadImage('xc:white');        # Make it all white
foreach my $ra_coord (@coords)                  # For every set of coords in the list
  {
  my($x,$y)     = @{$ra_coord}[0,1];		# Get the x and y
  $image	->Draw  (
      			 primitive	=> 'rectangle',
      			 points		=> (($x-3).','.($y-3).' '.($x+3).','.($y+3)),
      			 fill		=> 'red'
  			);			# Draw a small rectangle at each coord
  }
while (scalar(@coords) >= 4)                    # While there are atleast 4 entries in the list
  {
  for my $u (gradient(0,1,$rate))		# iterate from 0 to 1 in 500 steps
    {
    my($x,$y)	= &EvaluateCardinal2D(\@coords, $tension, $u);	# Hand paramaters to formula
    $image	->Set("pixel\[$x,$y\]"=>'red'); # Set that pixel red
    }
  shift(@coords);				# Remove the first entry of the list
  }
$image		->Write("Cardinal_Spline_Example.JPG"); # Save image

sub EvaluateCardinal2D
  {
  my($ra_coords,$T,$u)          = @_;
  my $s				= (1-$T)/2;
  my $u_matrix			= new Math::Matrix	# 4 x 1
  	(				# Matrix based off the point in the curve
  	 [($u ** 3),	($u ** 2),	($u),		(1)		]
  	);
  my $cardinal_matrix		= new Math::Matrix	# 4 x 4
  	(				# Guts of the Cardinal Spline formula
  	 [(-1 * $s),	(2 - $s),	($s - 2),	($s)		],
  	 [(2 * $s),	($s - 3),	(3-(2 * $s)),	(-1 * $s)	],
  	 [(-1 * $s),	(0),		($s),		(0)		],
  	 [(0),		(1),		(0),		(0)		],
  	);
  my $x_matrix			= new Math::Matrix	# 1 x 4
  	(				# X coords for point:
  	 [${${$ra_coords}[0]}[0]],	# 1
  	 [${${$ra_coords}[1]}[0]],	# 2
  	 [${${$ra_coords}[2]}[0]],	# 3
  	 [${${$ra_coords}[3]}[0]]	# 4
  	);
  my $y_matrix			= new Math::Matrix	# 1 x 4
  	(				# Y coords for point:
  	 [${${$ra_coords}[0]}[1]],	# 1
  	 [${${$ra_coords}[1]}[1]],	# 2
  	 [${${$ra_coords}[2]}[1]],	# 3
  	 [${${$ra_coords}[3]}[1]]	# 4
  	);
  my $xt	= int ($u_matrix * $cardinal_matrix * $x_matrix); # Compute for X
  my $yt	= int ($u_matrix * $cardinal_matrix * $y_matrix); # Compute for Y
  return($xt,$yt);
  }