#!/usr/bin/perl
##################
#
# Perform a "Discrete Cosine Transform" on an image,
# and output the result, scaled a bit for visibility.
# input:  tim.jpeg
# output: tim_DFT.jpeg
#
##################
use strict;
use warnings;
use GD;

my $input     = "tim.jpeg";
my $truecolor = 1;             # Yes; use 24-bits per pixel.
print "Reading in jpeg file '$input'... with 24 bits per pixel";
my $img    = GD::Image->newFromJpeg($input, $truecolor);
print "done.\n";

my ($width, $height) = $img->getBounds();
print " '$input' is $width x $height.\n";

my (@red, @green, @blue);
for my $w (0..$width-1){
  for my $h (0..$height-1){
    ($red[$w][$h], $green[$w][$h] , $blue[$w][$h]) = 
      $img->rgb( $img->getPixel($w,$h) );
  }
}

print "Doing DFT on rows...\n";
# Calculate the DFT matrix for the rows.
my @dft = ();
my $pi = 4*atan2(1,1);    # atan2(1,1) is pi/4, so...
for my $j (0..$width-1){
  for my $k (0..$width-1){
    $dft[$j][$k] = cos( $pi/$width * $j * ($k+0.5) );
  }
}
# Apply this DFT matrix to each row of the image, in each RGB channel.
for my $h (0..$height-1){
  my @newRed   = ();                # for each row,
  my @newGreen = ();
  my @newBlue  = ();
  for my $w (0..$width-1){
    $newRed[$w]  = 0;              # We're doing a sum,
    $newGreen[$w]= 0;              # so start off with zero.
    $newBlue[$w] = 0;
    for my $k (0..$width-1){
      $newRed[$w] += $red[$k][$h]     * $dft[$w][$k];
      $newGreen[$w] += $green[$k][$h] * $dft[$w][$k];
      $newBlue[$w] += $blue[$k][$h]   * $dft[$w][$k];
    }
  }
  for my $w (0..$width-1){         # Copy these transformed rows back.
    $red[$w][$h] = $newRed[$w];
    $green[$w][$h] = $newGreen[$w];
    $blue[$w][$h] = $newBlue[$w];
  }
}

print "Doing DFT on columns...\n";
# Calculate the DFT matrix for the columns.
@dft = ();
for my $j (0..$height-1){
  for my $k (0..$height-1){
    $dft[$j][$k] = cos( $pi/$height * $j * ($k+0.5) );
  }
}
# Apply same DFT matrix to each column of the image, in each RGB channel.
for my $w (0..$width-1){
  my @newRed   = ();                # for each row,
  my @newGreen = ();
  my @newBlue  = ();
  for my $h (0..$height-1){
    $newRed[$h]  = 0;              # We're doing a sum,
    $newGreen[$h]= 0;              # so start off with zero.
    $newBlue[$h] = 0;
    for my $k (0..$height-1){
      $newRed[$h] += $red[$w][$k]     * $dft[$h][$k];
      $newGreen[$h] += $green[$w][$k] * $dft[$h][$k];
      $newBlue[$h] += $blue[$w][$k]   * $dft[$h][$k];
    }
  }
  for my $h (0..$width-1){         # Copy these transformed columns back.
    $red[$w][$h] = $newRed[$h];
    $green[$w][$h] = $newGreen[$h];
    $blue[$w][$h] = $newBlue[$h];
  }
}


print "Normalizing...\n";
# Normalize to 0-255.
my $max = -1e6;
my $min =  1e6;
for my $w (0..$width-1){
  for my $h (0..$height-1){
    $max =  $red[$w][$h] if $red[$w][$h] > $max;
    $max =  $green[$w][$h] if $green[$w][$h] > $max;
    $max =  $blue[$w][$h] if $blue[$w][$h] > $max;
    $min =  $red[$w][$h] if $red[$w][$h] < $min;
    $min =  $green[$w][$h] if $green[$w][$h] < $min;
    $min =  $blue[$w][$h] if $blue[$w][$h] < $min;
  }
}
print "max = $max\n";
print "min = $min\n";
my $fudge = 500;   # saturate the highest pixels; push the rest.
my $convert = $fudge*255/($max-$min);   
for my $w (0..$width-1){
  for my $h (0..$height-1){
    $red[$w][$h] = $convert * ($red[$w][$h]-$min) % 255;
    $green[$w][$h] = $convert * ($green[$w][$h]-$min) % 255;
    $blue[$w][$h] = $convert * ($blue[$w][$h]-$min) % 255;
  }
}


# Let's see what we've got.
print "Forming new image...\n";
my $newImage = GD::Image->new($width,$height,$truecolor);
my $black = $newImage->colorAllocate(0,0,0);
my $white = $newImage->colorAllocate(255,255,255);
my $color;
for my $w (0..$width-1){
  for my $h (0..$height-1){
    $color=$newImage->colorAllocate($red[$w][$h],$green[$w][$h],$blue[$w][$h]);
    $newImage->setPixel($w,$h,$color);
  }
}
my $output = "tim_DFT.jpeg";
unlink $output;
print "Writing '$output'... ";
open    OUT, ">$output" or die "Oops couldn't create '$output'";
binmode OUT;
print   OUT $newImage->jpeg;
close   OUT;
print "done.\n";
