#!/usr/bin/perl -w
### An attempt to compress raw image data into the SPIFF format using the ISO standard JPEG algorithm
### Author: Gary Johnson
### 12.2.2001 - Project Begun

use strict;

die "Usage: jpeg.pl <input .bmp file> <output .jpg file>\n" if @ARGV != 2;

my @imgdata = ();
my @imgdata2 = ();
my ($rows, $cols);
$rows = 24;
$cols = 24;

####################################################################
### Read the data from INFILE into @imgdata

open (INFILE, "<$ARGV[0]") or die "Cannot read from $ARGV[0]: $!\n";
while (<INFILE>) { push(@imgdata, split(//, $_)); }
close (INFILE) or die "Cannot close $ARGV[0]: $!\n";

print "Bytes read in: ", scalar(@imgdata), "\n";

####################################################################
### Remove the BMP header (first 54 bytes) and save it in @bmheader

my @bmheader = ();

@bmheader[0..53] = @imgdata[0..53];
for (0..53) { shift @imgdata; }

print "Size of BitmapHeader: ", scalar(@bmheader), "\n";

####################################################################
### Parse the @imgdata into @R(ed)block, @G(reen)block, @B(lue)block

my @Rblock = ();
my @Gblock = ();
my @Bblock = ();

print "Splitting imgdata into Rblock, Gblock, Bblock...";
for (my $i=0; $i<$rows; $i++) {
	for (my $j=0; $j<$cols; $j++) {
		$Rblock[$i][$j] = $imgdata[3*($i*$cols+$j)];
		$Gblock[$i][$j] = $imgdata[3*($i*$cols+$j)+1];
		$Bblock[$i][$j] = $imgdata[3*($i*$cols+$j)+2];
	}
}
print "done\n";
print "Size of Rblock: ", scalar(@Rblock), "x", scalar(@{$Rblock[0]}), "\n";
printMatrix(\@Rblock, $rows, $cols);
print "Size of Gblock: ", scalar(@Gblock), "x", scalar(@{$Gblock[0]}), "\n";
printMatrix(\@Gblock, $rows, $cols);
print "Size of Bblock: ", scalar(@Bblock), "x", scalar(@{$Bblock[0]}), "\n";
printMatrix(\@Bblock, $rows, $cols);

####################################################################
### Split @Rblock into 8x8 pixel blocks

my @R8 = ();
my @G8 = ();
my @B8 = ();

print "Extracting an 8x8 pixel block from Rblock...";
for (my $i=0; $i<8; $i++) {
	for (my $j=0; $j<8; $j++) {
		$R8[$i][$j] = $Rblock[$i][$j];
	}
}
print "done\n";
print "Size of R8: ", scalar(@R8), "x", scalar(@{$R8[0]}), "\n";
printMatrix(\@R8, 8, 8);

####################################################################
### Create the Discrete Cosine Transform Matrix

my $pi = 4.0*atan2(1,1);  # typical way to get pi, using tan(pi/4)=1
my $piOver2N = $pi/(2*8);
my $invSqrt2 = 1/sqrt(2);
my $sqrt2overN = sqrt(2/8);
my $C = $invSqrt2;
my @dct = ();

print "Generating the DCT Matrix...";
for (my $i=0; $i<8; $i++) {
	for (my $j=0; $j<8; $j++) {
		$dct[$i][$j] = $sqrt2overN*$C*cos((2*$j+1)*$i*$piOver2N);
	}
	$C = 1;
}
print "done\n";
print "Size of DCT: ", scalar(@dct), "x", scalar(@{$dct[0]}), "\n";
printNumberMatrix(\@dct, 8, 8);

####################################################################
### Create the DCT Inverse Matrix

my @dctInv = ();

print "Generating the DCT Inverse Matrix...";
for (my $i=0; $i<8; $i++) {
	for (my $j=0; $j<8; $j++) {
		$dctInv[$i][$j] = $dct[$j][$i];
	}
}
print "done\n";
print "Size of DCT Inverse: ", scalar(@dct), "x", scalar(@{$dct[0]}), "\n";
printNumberMatrix(\@dctInv, 8, 8);

####################################################################
### Testing the matrixMultiply function

my @a = ([2,6],[3,5]);
my @b = ([8,2],[9,1]);
my @c = ();

print "Testing the matrixMultiply function...";
matrixMultiply(\@a, 2, 2, \@b, 2, 2, \@c);
print "done\n";
printNumberMatrix(\@a, 2, 2);
printNumberMatrix(\@b, 2, 2);
printNumberMatrix(\@c, 2, 2);

####################################################################
### Verify the matrices by checking that their product is the Identity
#
my @identity = ();

print "Multiplying dct and dctInv (This should be the Identity)...";
matrixMultiply(\@dct, 8, 8, \@dctInv, 8, 8, \@identity);
print "done\n";
printNumberMatrix(\@identity, 8, 8);

####################################################################
### DCT Transform - first pass
#
my @pixVec = ();
my @dctVec = ();
#
print "DCT Transform - first pass on 8x8 pixel block...";
for (my $i=0; $i<8; $i++) {
	for (my $j=0; $j<8; $j++) {
		$pixVec[0][$j] = ord($R8[$i][$j]);
	}
	matrixMultiply(\@dct, 8, 8, \@pixVec, 8, 1, \@dctVec);
	for (my $j=0; $j<8; $j++) {
		$R8[$i][$j] = $dctVec[0][$j];
	}
}
print "done\n";
print "Size of Transformed R8: ", scalar(@R8)*scalar(@{$R8[0]}), "\n";
printNumberMatrix(\@R8, 8, 8);

#
####################################################################
### DCT Transform - second pass
#
#
#
#
####################################################################
### Assemble the output data in @imgdata2

push(@imgdata2, @bmheader, @imgdata);

####################################################################
### Write the data from @imgdata2 to OUTFILE

open (OUTFILE, ">$ARGV[1]") or die "Cannot write to $ARGV[1]: $!\n";
print OUTFILE join('', @imgdata2);
close (OUTFILE) or die "Cannot close $ARGV[1]: $!\n";

print "Bytes written out: ", scalar(@imgdata2), "\n";

####################################################################
### Define Subroutines

sub matrixMultiply {
	my $m1 = shift;
	my $rows1 = shift;
	my $cols1 = shift;
	my $m2 = shift;
	my $rows2 = shift;
	my $cols2 = shift;
	my $product = shift;
	my $temp;

	die "Cannot multiply these two matrices: Sizes do not match.\n" unless $cols1 == $rows2;
	for (my $i=0; $i<$rows1; $i++) {
		for (my $j=0; $j<$cols2; $j++) {
			$temp = 0;
			for (my $k=0; $k<$rows2; $k++) {
				$temp += $$m1[$i][$k]*$$m2[$k][$j];
			}
			$$product[$i][$j] = $temp;
		}
	}
}

sub printMatrix {
	my $matrix = shift;
	my $rows = shift;
	my $cols = shift;
	for (my $i=0; $i<$rows; $i++) {
		print "[ ";
		for (my $j=0; $j<$cols; $j++) {
			print $$matrix[$i][$j], " ";
		}
		print "]\n";
	}
}

sub printNumberMatrix {
	my $matrix = shift;
	my $rows = shift;
	my $cols = shift;
	for (my $i=0; $i<$rows; $i++) {
		print "[ ";
		for (my $j=0; $j<$cols; $j++) {
			printf("% 7.5f ", $$matrix[$i][$j]);
		}
		print "]\n";
	}
}


####################################################################



