#   prog2p2.mpl -- Image Processing with Haar Wavelets
#   Written December 7, 2011 by Eric Olson for Math 761
#
#   This program was developed and tested using Maple 9.5
#   and should also work with other versions.

restart:
kernelopts(printbytes=false):
with(LinearAlgebra):
Digits:=15:
UseHardwareFloats:=true:

haar2d:=proc(X)
    local i,n,k,e,o,sqrt2,N;
    sqrt2:=sqrt(2.0);
    N:=RowDimension(X);
    e:=Vector(N/2);
    o:=Vector(N/2);
    for i from 1 to N do
        n:=N/2;
        while(n>=1) do
            for k from 1 to n do
                e[k]:=X[i,2*k-1];
                o[k]:=X[i,2*k]
            end;
            for k from 1 to n do
                X[i,k]:=(e[k]+o[k])/sqrt2;
                X[i,k+n]:=(e[k]-o[k])/sqrt2
            end;
            n:=n/2
        end
    end;
    for i from 1 to N do
        n:=N/2;
        while(n>=1) do
            for k from 1 to n do
                e[k]:=X[2*k-1,i];
                o[k]:=X[2*k,i]
            end;
            for k from 1 to n do
                X[k,i]:=(e[k]+o[k])/sqrt2;
                X[k+n,i]:=(e[k]-o[k])/sqrt2
            end;
            n:=n/2
        end
    end;
    X
end:

ihaar2d:=proc(X)
    local i,n,k,e,o,sqrt2,N;
    sqrt2:=sqrt(2.0);
    N:=RowDimension(X);
    e:=Vector(N/2);
    o:=Vector(N/2);
    for i from 1 to N do
        n:=1;
        while(n<=N/2) do
            for k from 1 to n do
                e[k]:=(X[i,k]+X[i,k+n])/sqrt2;
                o[k]:=(X[i,k]-X[i,k+n])/sqrt2
            end;
            for k from 1 to n do
                X[i,2*k-1]:=e[k];
                X[i,2*k]:=o[k]
            end;
            n:=n*2
        end
    end;
    for i from 1 to N do
        n:=1;
        while(n<=N/2) do
            for k from 1 to n do
                e[k]:=(X[k,i]+X[k+n,i])/sqrt2;
                o[k]:=(X[k,i]-X[k+n,i])/sqrt2
            end;
            for k from 1 to n do
                X[2*k-1,i]:=e[k];
                X[2*k,i]:=o[k]
            end;
            n:=n*2
        end
    end;
    X
end:

printf("prog2p2.mpl -- Image Processing with Haar Wavelets\n"
       "Written December 7, 2011 by Eric Olson for Math 761\n\n"):

N:=512:
epsilon:=0.1:

printf("Reading image1.dat...\n"):
fd:=fopen("image1.dat",READ):
buf:=readdata(fd,float,N):
fclose(fd):
y:=Matrix(buf,datatype=float):

printf("Processing...\n"):
y:=haar2d(y):
for n from 1 to N do
    for m from 1 to N do
        if abs(y[n,m])<epsilon
        then
            y[n,m]:=0.0
        end
    end
end:
y:=ihaar2d(y):

printf("Writing test.dat...\n"):
buf:=[seq([seq(y[i,j],j=1..N)],i=1..N)]:
fd:=fopen("test.dat",WRITE):
writedata(fd,buf,float,N):
fclose(fd):