#!/usr/bin/perl

$LOG2 =	0.693147180559945309417232121458176568075500134360255254120680009;
$PI =	3.141592653589793238462643383279502884197169399375105820974944592;

$InFile =	$ARGV[0].".raw";
$OutFile =	$ARGV[0].".wavelet.raw";
$NewFile =	$ARGV[0].".recovered.raw";

print ("encoding: $InFile -> $OutFile \n");

$Char;
open(IN, $InFile);
binmode(IN);
$Size =		sqrt(-s IN);
print $Size;
for my $A (0 .. $Size - 1) {
	for my $B (0 .. $Size - 1) {
		read(IN, $Char, 1);
		$Data->[$A][$B] =	ord($Char);
	}
}
close(IN);

#$Data =	[[64,128,64,128],[128,192,64,128],[128,64,64,128],[192,128,64,64]];
#$Data =		[[4,4,14,14],[4,4,14,14],[14,14,24,24],[14,14,24,24]];
#$Data =	[[7,12],[34,21]];

print "\nDone reading the data...\n";

$Power =	log($#$Data + 1) / $LOG2;
$CurrentS = $Data;
$Out;
for my $G (0 .. $Power - 1) {
	my $N =	2**($Power - $G) - 1;
	my $M =	2**($Power - $G - 1) - 1;
	$LastS =	$CurrentS;
	for my $A (0 .. $M) {
		for my $B (0 .. $M) {
#			print "$A\t$B\n";
			my $NowH =	($LastS->[2*$A+1][2*$B+1] + $LastS->[2*$A][2*$B+1] - $LastS->[2*$A+1][2*$B] - $LastS->[2*$A][2*$B]) / 4;
			my $NowV =	($LastS->[2*$A+1][2*$B+1] - $LastS->[2*$A][2*$B+1] + $LastS->[2*$A+1][2*$B] - $LastS->[2*$A][2*$B]) / 4;
#			my $NowD =	($LastS->[2*$A+1][2*$B+1] - $LastS->[2*$A][2*$B+1] - $LastS->[2*$A+1][2*$B] + $LastS->[2*$A][2*$B]) / 4;
			my $NowD =	(-$LastS->[2*$A+1][2*$B+1] + $LastS->[2*$A][2*$B+1] + $LastS->[2*$A+1][2*$B] - $LastS->[2*$A][2*$B]) / 4;
			my $NowS =	($LastS->[2*$A+1][2*$B+1] + $LastS->[2*$A][2*$B+1] + $LastS->[2*$A+1][2*$B] + $LastS->[2*$A][2*$B]) / 4;
			$CurrentS->[$A][$B] =	$NowS;
			$Out->[$A][($N+1)/2 + $B] =				$NowH;
			$Out->[($N+1)/2 + $A][$B] =				$NowV;
			$Out->[($N+1)/2 + $A][($N+1)/2 + $B] =	$NowD;
			if($G == $Power - 1) {
				$Out->[0][0] =						$LastS->[0][0];
			} else {
				$Out->[$A][$B] =					$NowS;
			}
		}
		print "\r$A\t\t";
	}
}

print "\nDone calculating the wavelet coëfficients...";

#Displayer($Out);

open(OUT, '>'.$OutFile);
binmode(OUT);
for my $A (0 .. $#$Out) {
	for my $B (0 .. $#$Out) {
		if($A == 0 && $B == 0) {
			print OUT chr(Chop($Out->[$A][$B]));
		} else {
## n^3 ##
			my $Cube =	$Out->[$A][$B];
			if($Cube < 0) {
				$Cube =	-128 * ((-$Cube / 128) ** (1/5));
			} else {
				$Cube =	128 * (($Cube / 128) ** (1/5));
			}
			print OUT chr(Chop($Cube + 128));
## SQRT ## 
#			my $Sqrtized =	$Out->[$A][$B];
#			if($Sqrtized < 0) {
#				$Sqrtized =	-128 * sqrt(-$Sqrtized / 128);
#			} else {
#				$Sqrtized =	128 * sqrt($Sqrtized / 128);
#			}
#			print OUT chr(Chop($Sqrtized + 128));
## Linear ##
#			print OUT chr(Chop($Out->[$A][$B] + 128));
		}
	}
}
close(OUT);

#undef $Out;
undef $CurrentS;
undef $LastS;
undef $Data;

print "\nDone with the forward transform...\n";

##########
## 5UB5 ##
##########

sub Chop {
	if($_[0] > 255) {
		return 255;
	} elsif($_[0] < 0) {
		return 0;
	} else {
		return int($_[0]);
	}
}