#!/usr/bin/perl
##################
#
# Blur and image by doing a DFT, a filter, and an inverse DFT.
# input: tim.jpeg
# outpu: tim_blur.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";
#print " and has " . $img->colorsTotal . " colors allocated.\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) );
  }
}

# Used by all of what follows.
my @dft;
my $pi = 4*atan2(1,1);    # atan2(1,1) is pi/4, so...



# ---- DFT rows ----------------------------------

print "Doing DFT on rows...\n";
# Calculate the DFT matrix for the rows.
@dft = ();
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];
      $newBlue[$w] += $blue[$k][$h]   * $dft[$w][$k];
      $newGreen[$w] += $green[$k][$h] * $dft[$w][$k];
    }
  }
  for my $w (0..$width-1){         # Copy these transformed rows back.
    $red[$w][$h] = $newRed[$w];
    $blue[$w][$h] = $newBlue[$w];
    $green[$w][$h] = $newGreen[$w];
  }
}

# ----- DFT columns ------------------------------------------------

# Apply same DFT matrix to each column of the image, in each RGB channel.
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) );
  }
}
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];
      $newBlue[$h] += $blue[$w][$k]   * $dft[$h][$k];
      $newGreen[$h] += $green[$w][$k] * $dft[$h][$k];
    }
  }
  for my $h (0..$height-1){         # Copy these transformed columns back.
    $red[$w][$h] = $newRed[$h];
    $green[$w][$h] = $newGreen[$h];
    $blue[$w][$h] = $newBlue[$h];
  }
}

# === blur ==========================================================

my $cutoff = $width/3;
print "Appling radial DFT chop at radius=$cutoff .. \n";
for my $w (0..$width-1){
  for my $h (0..$height-1){
    my $radius = sqrt($w**2 + $h**2);
    if ($radius > $cutoff){
      $red[$w][$h] =   0;
      $green[$w][$h] = 0;
      $blue[$w][$h] =  0;
    }
  }
}

# ----- inverse DFT rows -------------------------------------------

print "Doing inverse DFT on rows...\n";
# Calculate the DFT matrix for the rows.
@dft = ();
for my $j (0..$width-1){
  for my $k (0..$width-1){
    $dft[$j][$k] = cos( $pi/$width * ($j+0.5) * $k );
  }
}
# 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]  = $red[0][$h]/2;              # We're doing a sum;
    $newGreen[$w]= $green[0][$h]/2;            # inverse starts here.
    $newBlue[$w] = $blue[0][$h]/2;
    for my $k (1..$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];
    }
  }
  my $scale = 2/$width;
  for my $w (0..$width-1){         # Copy these transformed rows back.
    $red[$w][$h] = $scale*$newRed[$w];
    $green[$w][$h] = $scale*$newGreen[$w];
    $blue[$w][$h] = $scale*$newBlue[$w];
  }
}

# --- inverse DFT columns --------------------------------------

print "Doing inverse 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+0.5) * $k );
  }
}
# 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]  = $red[$w][0]/2;          # Inverse transform
    $newGreen[$h]= $green[$w][0]/2;       # starts here rather than at zero.
    $newBlue[$h] = $blue[$w][0]/2;
    for my $k (1..$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];
    }
  }
  my $scale = 2/$height;
  for my $h (0..$height-1){         # Copy these transformed columns back.
    $red[$w][$h] = $scale*$newRed[$h];
    $green[$w][$h] = $scale*$newGreen[$h];
    $blue[$w][$h] = $scale*$newBlue[$h];
  }
}

# Let's see what we've got.
print "Forming new image...\n";
my $newImage = GD::Image->newTrueColor($width,$height);
#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_blur.jpeg";
unlink $output;
my ($width2, $height2) = $newImage->getBounds();
print " Output image is $width2 x $height2 .\n";
#print " and has " . $newImage->colorsTotal . " colors allocated.\n";
print "Writing '$output'... ";
open    OUT, ">$output" or die "Oops couldn't create '$output'";
binmode OUT;
print   OUT $newImage->jpeg;
close   OUT;
print "done.\n";
