function fn = build(h, middle_index)
    len = length(h);
    max = len - middle_index;
    min = 1 - middle_index;

    function v = getter(in)
        if isnumeric(in)
            v = h(in + middle_index);
        else
            if isequal(in,'max')
                v = max;
            else
                v = min;
            end
        end
    end
    fn = @getter;
end

function z = correlate(f,g,j)
    z = 0;
    min = g('min')-j;
    max = g('max')-j;
    for k = min:max
        z = z + mirror(f,k)*g(k+j);
    end
end

function im2 = dwt(im, h, g, numOfIterations)
    n = size(im,1);
    i = 0;
    im2 = double(im);
    while i < numOfIterations
        w = 1:n/2^i;
        for r = w
            im2(r,w) = dwt1d(im2(r,w),h,g);
        end
        for c = w
            im2(w,c) = dwt1d(im2(w,c),h,g);
        end
        i = i + 1;
    end
end

function im2 = dwt1d(im, h, g)
    n = length(im);
    im2 = zeros(size(im));
    for j = 1:n/2
        im2(j) = correlate(im,h,-2*j);
        im2(n/2+j) = correlate(im,g,-2*j);
    end
end

function im = dwtinv( im, h, g, numOfIterations )
    n = size(im,1);
    i = numOfIterations - 1;
    while i >= 0
        w = 1:n/2^i;
        for r = w
            im(r,w) = dwtinv1d(im(r,w),h,g);
        end
        for c = w
            im(w,c) = dwtinv1d(im(w,c),h,g);
        end
        i = i - 1;
    end
    im = uint8(im);
end

function im = dwtinv1d(im, h, g)
    n = length(im);
    f = im(1:n/2);
    alpha = im(n/2+1:n);

    ming = g('min');
    maxg = g('max');
    minh = h('min');
    maxh = h('max');
    im(:) = 0;
    for j = 1:n
        for k = ceil((-maxh+j)/2) : floor((-minh+j)/2)
            im(j) = im(j) + mirror(f,k)*h(j-2*k);
        end
        for k = ceil((-maxg+j)/2) : floor((-ming+j)/2)
            im(j) = im(j) + mirror(alpha,k)*g(j-2*k);
        end
    end
end

function [h,ht,g,gt] = haar()
    t = 1/sqrt(2);
    h = build(t*[1,1], 1);
    ht = h;
    g = build(t*[1,-1], 1);
    gt = g;
end

function v = mirror(f, k)
    len = length(f);
    if k < 1
        k = 2 - k;
    elseif k > len
        k = len - (k - len);
    end
    v = f(k);
end

function im2 = runforward(im, numOfIterations)
    %im = rgb2gray(imread('lichtenstein.png'));
    [~,ht,~,gt] = haar;
    im2 = dwt( im, ht, gt, numOfIterations );
    %figure;
    %dwtshow(im2, numOfIterations);
end

function im3 = runbackward(im2, numOfIterations)
    [h,~,g,~] = haar;
    im3 = dwtinv(im2, h, g, numOfIterations);
    figure;
    imshow(im3);
end