%
%  Developed by Timothy D. Dorney
%               Rice University
%               December, 1999
%               tdorney@ieee.org
%
%  Coded using MATLAB 5.X.X with NO additional toolboxes.
%

% Original phantom
% 'img'
load img_orig_512.mat

% Reconstructed image with 800 angles, detector length = 1.5,
% 1000 receivers (projection lines), final image is 512x512 pixels
% 'rimage3'
load reconst_image_800_1.5_1000_512.mat

% Squared error for entire images
err1_img = (img - rimage3).^2;
err1 = sum(sum(err1_img));
disp('Squared error for entire image =')
disp(err1)

% imagesc(err1_img)

in = rimage3;
[r,c] = size(in);

% scale = 1/sqrt(2);
scale = 1;

% USER DEFINED threshold to determine where the edges are located.
tol = 0.1;

% Increase SPAN to 5 from 3 to test 2 adjacent pts. on either side
% of peak in HP instead of 1 adjacent pt.
% It does not improve the result!
SPAN = 3;

% ========================================
% HAAR WAVELET TRANSFORM => VERTICAL LINES
t1 = zeros(r,c);
t2 = zeros(r,c);
for k = 1:1:r
   for m = 1:1:(c-1)
      t1(k,m) = scale*(in(k,m) + in(k,(m+1)));
      t2(k,m) = scale*(-in(k,m) + in(k,(m+1)));
   end
   t1(k,c) = scale*(in(k,c));
   t2(k,c) = scale*(-in(k,c));
end

q2 = zeros(r,c+4);
q2(1:r,3:(c+2)) = abs(t2)./(max(max(abs(t2))));
t3 = zeros(r,c+4);
t3(1:r,3:(c+2)) = t2;
tz = zeros(r,c+4);
tz(1:r,3:(c+2)) = t1;

for x = 1:1:r
   for y = 3:1:(c+2)
      VFLAG = [1 1 1 1 1];
      if (q2(x,y) >= tol)
         if ((t3(x,y) > t3(x,y-1)) & (t3(x,y) > t3(x,y+1)))
            if (t3(x,y-1) > 0)
               t3(x,y) = t3(x,y) + t3(x,y-1);
               tz(x,y-1) = tz(x,y-1) - t3(x,y-1);
               tz(x,y) = tz(x,y) - t3(x,y-1);
               VFLAG(2) = 0;
            end
            if (t3(x,y+1) > 0)
               t3(x,y) = t3(x,y) + t3(x,y+1);
               tz(x,y+1) = tz(x,y+1) + t3(x,y+1);
               tz(x,y) = tz(x,y) + t3(x,y+1);
               VFLAG(4) = 0;
            end
            if ((t3(x,y-1) > t3(x,y-2)) & (t3(x,y-2) > 0) & (SPAN == 5))
               t3(x,y) = t3(x,y) + t3(x,y-2);
               VFLAG(1) = 0;
            end
            if ((t3(x,y+1) > t3(x,y+2)) & (t3(x,y+2) > 0) & (SPAN == 5))
               t3(x,y) = t3(x,y) + t3(x,y+2);
               VFLAG(5) = 0;
            end
            t3(x,y-2:y+2) = VFLAG.*t3(x,y-2:y+2);
         elseif ((t3(x,y) < t3(x,y-1)) & (t3(x,y) < t3(x,y+1)))
            if (t3(x,y-1) < 0)
               t3(x,y) = t3(x,y) + t3(x,y-1);
               tz(x,y-1) = tz(x,y-1) - t3(x,y-1);
               tz(x,y) = tz(x,y) - t3(x,y-1);
               VFLAG(2) = 0;
            end
            if (t3(x,y+1) < 0)
               t3(x,y) = t3(x,y) + t3(x,y+1);
               tz(x,y+1) = tz(x,y+1) + t3(x,y+1);
               tz(x,y) = tz(x,y) + t3(x,y+1);
               VFLAG(4) = 0;
            end
            if ((t3(x,y-1) < t3(x,y-2)) & (t3(x,y-2) < 0) & (SPAN == 5))
               t3(x,y) = t3(x,y) + t3(x,y-2);
               VFLAG(1) = 0;
            end
            if ((t3(x,y+1) < t3(x,y+2)) & (t3(x,y+2) < 0) & (SPAN == 5))
               t3(x,y) = t3(x,y) + t3(x,y+2);
               VFLAG(5) = 0;
            end
            t3(x,y-2:y+2) = VFLAG.*t3(x,y-2:y+2);
         end
      end
   end
end
t4 = [];
t4 = t3(1:r,3:c+2);
ty = [];
ty = tz(1:r,3:c+2);

% INVERSE HAAR WAVELET TRANSFORM

t5 = zeros(r,c);
for k = 1:1:r
   for m = 1:1:c
      t5(k,m) = (scale*(ty(k,m) - t4(k,m)))./2;
   end
end

% ==========================================
% HAAR WAVELET TRANSFORM => HORIZONTAL LINES
t6 = zeros(r,c);
t7 = zeros(r,c);
for k = 1:1:(r-1)
   for m = 1:1:c
      t6(k,m) = scale*(t5(k,m) + t5(k+1,m));
      t7(k,m) = scale*(-t5(k,m) + t5(k+1,m));
   end
end
for k = r
   for m = 1:1:c
      t6(k,m) = scale*(t5(k,m));
      t7(k,m) = scale*(-t5(k,m));
   end
end

q2 = zeros(r+4,c);
q2(3:(r+2),1:c) = abs(t7)./(max(max(abs(t7))));
t3 = zeros(r+4,c);
t3(3:(r+2),1:c) = t7;
tx = zeros(r+4,c);
tx(3:(r+2),1:c) = t6;

for x = 3:1:(r+2)
   for y = 1:1:c
      VFLAG = [1 1 1 1 1]';
      if (q2(x,y) >= tol)
         if ((t3(x,y) > t3(x-1,y)) & (t3(x,y) > t3(x+1,y)))
            if (t3(x-1,y) > 0)
               t3(x,y) = t3(x,y) + t3(x-1,y);
               tx(x-1,y) = tx(x-1,y) - t3(x-1,y);
               tx(x,y) = tx(x,y) - t3(x-1,y);
               VFLAG(2) = 0;
            end
            if (t3(x+1,y) > 0)
               t3(x,y) = t3(x,y) + t3(x+1,y);
               tx(x+1,y) = tx(x+1,y) + t3(x+1,y);
               tx(x,y) = tx(x,y) + t3(x+1,y);
               VFLAG(4) = 0;
            end
            if ((t3(x-1,y) > t3(x-2,y)) & (t3(x-2,y) > 0) & (SPAN == 5))
               t3(x,y) = t3(x,y) + t3(x-2,y);
               VFLAG(1) = 0;
            end
            if ((t3(x+1,y) > t3(x+2,y)) & (t3(x+2,y) > 0) & (SPAN == 5))
               t3(x,y) = t3(x,y) + t3(x+2,y);
               VFLAG(5) = 0;
            end
            t3(x-2:x+2,y) = VFLAG.*t3(x-2:x+2,y);
         elseif ((t3(x,y) < t3(x-1,y)) & (t3(x,y) < t3(x+1,y)))
            if (t3(x-1,y) < 0)
               t3(x,y) = t3(x,y) + t3(x-1,y);
               tx(x-1,y) = tx(x-1,y) - t3(x-1,y);
               tx(x,y) = tx(x,y) - t3(x-1,y);
               VFLAG(2) = 0;
            end
            if (t3(x+1,y) < 0)
               t3(x,y) = t3(x,y) + t3(x+1,y);
               tx(x+1,y) = tx(x+1,y) + t3(x+1,y);
               tx(x,y) = tx(x,y) + t3(x+1,y);
               VFLAG(4) = 0;
            end
            if ((t3(x-1,y) < t3(x-2,y)) & (t3(x-2,y) < 0) & (SPAN == 5))
               t3(x,y) = t3(x,y) + t3(x-2,y);
               VFLAG(1) = 0;
            end
            if ((t3(x+1,y) < t3(x+2,y)) & (t3(x+2,y) < 0) & (SPAN == 5))
               t3(x,y) = t3(x,y) + t3(x+2,y);
               VFLAG(5) = 0;
            end
            t3(x-2:x+2,y) = VFLAG.*t3(x-2:x+2,y);
         end
      end
   end
end
t8 = [];
t8 = t3(3:r+2,1:c);
tw = [];
tw = tx(3:r+2,1:c);

% INVERSE HAAR WAVELET TRANSFORM

t9 = zeros(r,c);
for k = 1:1:r
   for m = 1:1:c
      t9(k,m) = (scale*(tw(k,m) - t8(k,m)))./2;
   end
end

err2_img = (img - t9).^2;
err2 = sum(sum(err2_img));
fin = t9;

disp('The variable "fin" has the localized edge enhanced image.')
disp(' ')
disp('Squared error for entire image after edge enhancement =')
disp(err2)

