#!/usr/bin/perl

##[[[[[[[[[[[[[[[[[[[[[[[
#Blur.pl
#takes a jpeg and, applies a cosine transformation, shows the transform,
# filters the transformation, and then displays the filtered image.
#NOTE: this is very slow and should only be used on small images. (128x192ish)
#
# # Modified by Jim M to run as a command line script.
# # Invoke from command line as "blur.pl test3.jpg  > output.jpg"
#
#
#by Ian Smith-Heisters ian@0x09.com
#March 8th, 2004
##]]]]]]]]]]]]]]]]]]]]]]]

use CGI qw(:standard);
use CGI::Carp qw(fatalsToBrowser);
use GD;
use Math::Matrix;
use Math::Trig;
use strict;
use warnings;

## Jim - get filename from command line
my $input_jpeg = shift;
die "Usage: ./blur.pl input.jpg > output.jpg" unless defined $input_jpeg;
print " input file is '$input_jpeg' \n";
exit;

#initialize global variable hash
my %g = ( filename => $input_jpeg );

#set truecolor
my $truecolor = 1;		# on

#Load any global variables from the URL
for my $key (keys %g){
  $g{$key} = param($key) if defined param($key);
}

#create an image from filename
my $image = GD::Image -> newFromJpeg($g{filename}, $truecolor) || 
  die "couldn't load original jpeg";

#get the width and height of the image
#  getBounds returns a two member list (width,height)
my ($width,$height) = $image -> getBounds();

### --------------------------------------------------
### Jim : Math::Matrix documentation doesn't discuss
###       creating empty matrices and then growing them -
###       so the following might work, but you're walking on thin ice.

# +++ returns a hash containing three arrays, each containing R, G, or B 
# information for each pixel in the image
sub get_channels {
  my $img = shift;
  my %temp;
  $temp{r} = new Math::Matrix ([]); # these all need to be matrices so we can
  $temp{g} = new Math::Matrix ([]); #  easily do matrix operations on them
  $temp{b} = new Math::Matrix ([]); #  using Math::Matrix.

  for (my $x = 0;$x<$height;$x++){ #start by looping over the columns
    for (my $y = 0;$y<$width;$y++){ #then loop over each row

      my $index = $img -> getPixel($x,$y); # get the index of the pixel
      my @rgb = $image->rgb($index); # coinvert the index to RGB format

      my $red = shift @rgb;	# put the red information in red
      my $green = shift @rgb;	# put the g info in green
      my $blue = shift @rgb;	# and (surprise!) put the b info in blue

      #throw red green and blue into the hash at their pixel coordinate
      $temp{r}[$x][$y] = $red;
      $temp{g}[$x][$y] = $green;
      $temp{b}[$x][$y] = $blue;
    }
  }
  return \%temp;		# return a ref to the hash
}

# +++ return an image given a ref to an rgb hash
sub get_image {
  my $rgb = shift;
  my $img = GD::Image->new($width, $height, $truecolor) || 
    die "Couldn't create the new jpeg."; # create a new image

  for (my $x=0; $x<$height; $x++) { # cycle through x pixels
    for (my $y=0; $y<$width; $y++){ # cycle through y pixels

      # create a list for the pixel we're working on from the rgb hash info
      my @rgbList = ($rgb->{r}[$x][$y], 
		     $rgb->{g}[$x][$y], 
		     $rgb->{b}[$x][$y]);

      # set the pixel to the data in rgbList
      $img -> setPixel($x, $y, $img->colorClosest(@rgbList));
    }
  }
  return $img;			# return the image.
}

# ------ Test all that. This shouldn't do anything more than display test.jpg
#print "Content-type: image/jpg\n\n";
#print get_image(get_channels($image))->jpeg;
# ------ end of tests

# SO! Now we can get an image, extract the rgb info into arrays, and put
#  that information back into a new image.

# Now we need to be able to take the rgb info and apply the cosine
#  transformation. So. First a function to make the cosine transformation
#  matrix:


# Generate a DCT matrix for dimensions (rows,columns). 
# This is defined in the Wikipedia as DCT-II.
sub get_dctm {
  my $n = shift;
  my $pi_div_n = pi/$n;  ## Jim: replaced pi/$n; renamed (j,k)->(row,col)
  my $dctm = new Math::Matrix ([]);
  for my $row (0..$n-1){
    for my $col (0..$n-1){
      $dctm->[$row][$col] = cos ($pi_div_n * $row * ($col + 0.5));
    }
  }
  return $dctm;
}

# Generate the inverse of the DCT. 
# This is defined in the Wikipedia as DCT-III multiplied by 2/n.
sub get_idctm {
  my $n = shift;
  my $pi_div_n = pi/$n;
  my $one_div_n = 1/$n;
  my $idctm = new Math::Matrix ([]);

  # Jim: Changed column index to start at 1.
  #      Replaced pi/$n with $pi_div_n
  for my $row (0..$n-1){
    for my $col (1..$n-1){
      $idctm->[$row][$col] = (2/$n) * cos ($pi_div_n * $col * ($row + 0.5));
    }
  }

  # Jim: The k=0 term is different; left column of matrix is just 0.5 * 2/$n
  for my $row (0..$n-1){
    $idctm->[$row][0] = $one_div_n;
  }
  
  return $idctm;
}



# So. We need to make two DCTMs. Each is a square, one with dimensions equal
#  to the number of rows in the image, and the other with dimensions equal to
#  the number of columns in the image. Then we need to step through each row
#  in the separated RGB Channels hash made by get_colors and multiply it by
#  the DCT matrix. Same for columns. Yay.

# First, a function that takes a channel hash and multiplies the rows and 
#  columns by the two hashes supplied. To invert a DCT, supply 1 as this sub's
#  4th argument.
sub dct_transform {
  my ($hashref, $rowdctm, $coldctm, $inverse) = @_;
  my ($w, $h) = $image->getBounds(); # get the height and width.

  # Return if the rows and cols in the matrices don't match the image size
  my @rowdim = $rowdctm->size;
  my @coldim = $coldctm->size;
  return undef unless ($rowdim[0] == $width && $coldim[0] == $height);
  my @hs = $hashref->{r}->size;
#  print "Image width $width, DCTM width $rowdim[0], Hash width $hs[1]\n";
#  print "Image height $height, DCTM height $coldim[0], Hash height $hs[0]\n"; 

#  print "Starting transform...\n" if !($inverse);
#  print "Starting inverse transformation...\n" if ($inverse);

  for my $key (keys %{$hashref}) { # we'll do this channel by channel
#    print "Starting rows.\n";
    for my $row (0..$h-1) {
      my $v = new Math::Matrix ($hashref->{$key}[$row]);
      $v = $v->multiply($rowdctm);
      # The inverse equation starts by adding (x0)/2. This should add that to
      #  the transformation
      if ($inverse) {		# if we're inverting do the following
	my ($x,$y) = $v->size;
	my $q = $v->[0][0];	# need to get this out before it changes
	for my $g (0..$y-1) {
	  $v->[0][$g] += $q/2;	# since $v is one dimensional use [0]
	}
      }
      $hashref->{$key}[$row] = $v->[0];	# put the transformed vector back
    }
#    print "Done with rows for $key channel.\n";

#    print "Starting columns\n";
    my $tm = $hashref->{$key}->transpose;
    for my $col (0..$w-1) { # since it's transposed, the columns work like rows
      my $v = new Math::Matrix ($tm->[$col]);
      $v = $v->multiply($coldctm);
      # The inverse equation starts by adding (x0)/2. This should add that to
      #  the transformation
      if ($inverse) {		# if we're inverting do the following
	my ($x,$y) = $v->size;
	my $q = $v->[0][0];
	for my $g (0..$y-1) {
	  $v->[0][$g] += $q/2; # add half of x0 to the vector
	}
      }

      $tm->[$col] = $v->[0];	# put the contents of the vector back
    }
    $hashref->{$key} = $tm->transpose; # flip it back and throw it in!
#    print "Done with columns for $key channel.\n"
  }
  return $hashref;
}


# If that works, we can apply some kind of radial cutoff to blur the image.
# This filter is destructive, doing it on the hash passed in, not a copy of the
#  hash.
sub hip_filter {
#  print "Applying Hi-pass filter to image.\n";
  my ($hashref, $cutoff) = @_;
  for my $h (0..$height-1) {
    for my $w (0..$width-1) {
      my $radius = sqrt($w**2 + $h**2);
      if ($radius > $cutoff) {
	for my $key (keys %{$hashref}) {
	  $hashref->{$key}[$h][$w] = 0;
	}
      }
    }
  }
}

# Try passing it the image, a row-dctm and a column-dctm, then the same with
#  inverted dctms to see if we get the same image back.

# get the rgb info
my $channels = get_channels($image);
#print "Got channels.\n";

# get the DCT Matrices for rows and columns.
my $row_dctm = get_dctm($width,$width);
my $col_dctm = get_dctm($height,$height);
#print "Got dctms.\n";

# get the inverse DCT Matrices for rows and columns
my $row_idctm = get_idctm($width,$width);
my $col_idctm = get_idctm($height,$height);
#print "Got idctms\n";

# transform the channel info with the DCTMs
my $transformed = dct_transform($channels,$row_dctm,$col_dctm);

#filter the transformed channel information
#hip_filter($transformed, $width/3);

# untransform the channel info with the inverse DCTMs
my $untransformed = dct_transform($transformed, $row_idctm, $col_idctm,1);

# get the image from the channel info.
my $new_image = get_image($untransformed);
#print %{(get_channels($new_image))};

# print "Content-type: image/jpg\n\n";   ## Jim - just output jpeg
print $new_image->jpeg;
