close all;
clear;
clc;

original_img = imread('102_6.tif');
original_img = im2double(original_img);
original_img = original_img(:, 1:end-2);
img = imgaussfilt(original_img).*0.1; % 0.4 % 0.1
ramp = linspace(0, 0.85, size(img, 2)); % 0.55 % 0.85
img = img + repmat(ramp, size(img, 1), 1);

%{
figure(1);
subplot(1, 2, 1);
imshow(original_img);
title('original image');
subplot(1, 2, 2);
imshow(img);
title('noisy image, this should be processed')
%}

% ---
% please do not edit the code above, write your solution under this line,
% please handle original_img as unknown
% ---

%parameters
sigma = 3.8;
gamma = 1.0;
lambda = 8.5;

figure(2);

img_gabor = zeros(334, 346);

%application of the Gabor kernel
for i = 0:7
    theta = pi*i/8;
    [gb_r, gb_i] = get_gabor_kernel(theta, sigma, gamma, lambda);
    out_img = apply_gabor_kernel(img, gb_r, gb_i);
    
    subplot(4, 4, 2*i + 1);
    imagesc(gb_r);
    colormap('gray');
    axis equal;
    title('real part of Gabor, theta = ', num2str(theta));

    subplot(4, 4, 2*i + 2)
    imagesc(out_img);
    colormap('gray');
    axis equal;
    title('result of the specific filtering');
    
    %summing
    img_gabor = img_gabor + out_img;
    
end

%normalization
img_gabor = img_gabor/max(img_gabor(:));

figure(3);

subplot(1, 2, 1);
imshow(img);
title('noisy input');
subplot(1, 2, 2);
imagesc(img_gabor);
colormap('gray');
axis('equal');
title('Gabor filtered, normilized');

%calculating the frequency spectrum
FT = fftshift(fft2(img_gabor));

%calculating the magnitude part
magn = log(abs(FT)+1);
M_FT = (magn - min(magn(:))) / max(magn(:)) - min(magn(:));

%zeroing out + inverse transformation
zeroed_out = FT;
zeroed_out(160:180, 160:180) = 0;
I = abs(ifft2(ifftshift(zeroed_out)));

%thresholding + zeroing out
loc = I>=0.1;
I_m = I;
I_m(loc) = 255;
I_m(~loc) = 0;
I_m(:, 1:7)= 0;
I_m(:, end-7:end)=0;

%dilation with a structuring element
structuring_el = strel('disk', 11, 8);
I_m_2 = imdilate(I_m, structuring_el);

%masked version of the Gabor filtered image
img_gabor_masked =  img_gabor + ~I_m_2;

%bradley local thresholding
bradley_thresholded = bradley(img_gabor_masked, [15, 15], 10);

%erosion
I_m_2_eroded = imerode(I_m_2, structuring_el);
bradley_eroded = bradley_thresholded + ~I_m_2_eroded;
bradley_eroded = bradley(bradley_eroded, [15, 15], 10);

figure(4)

subplot(2,3,1);
imshow(M_FT);
title('magnitude of FFT');

subplot(2,3,2);
imshow(I_m);
title('baseline mask');

subplot(2,3,3);
imshow(I_m_2);
title('dilated mask');

subplot(2,3,4);
imagesc(img_gabor_masked);
title('masked version of the Gabor-filtered');
colormap('gray');

subplot(2,3,5);
imagesc(bradley_thresholded);
colormap('gray');
title('after Bradley');

subplot(2,3,6);
imagesc(bradley_eroded);
colormap('gray');
title('boundaries of Bradley corrected');

%filling in the holes
bw = bradley_eroded;
filled = imfill(~bw, 'hole');
holes = filled & bw;
bigholes = bwareaopen(holes, 20);
smallholes = holes & ~bigholes;
small_holes_removed = bw & ~smallholes;

figure(5)
subplot(2,3,1);
imshow(bw);
title('boundaries of Bradley corrected');

subplot(2,3,2);
imshow(filled);
title('filled holes');

subplot(2,3,3);
imshow(holes);
title('holes');

subplot(2,3,4);
imshow(bigholes);
title('big holes');

subplot(2,3,5);
imshow(smallholes);
title('small holes');

subplot(2,3,6);
imshow(small_holes_removed);
title('small holes removed');

%skeletonization
skelton = bwmorph(~small_holes_removed, 'skel', 'Inf');
skeltonized = ~skelton;

%deleting small segments
deleted = regionprops(~skeltonized, 'Area', 'PixelIdxList');
cleaned = ~skeltonized;
for k = 1:size(deleted, 1)
   if(deleted(k).Area < 20)
      cleaned(deleted(k).PixelIdxList) = 0;
   end
end

%detection of endpoints and branchpoints
endpoints = bwmorph(cleaned, 'endpoints');
[ep_x, ep_y] = ind2sub([334 346], find(endpoints > 0));

branchpoints = bwmorph(cleaned, 'branchpoints');
[bp_x, bp_y] = ind2sub([334 346], find(branchpoints > 0));

%finding the real endpoints and branchpoints
endpoints = padarray(endpoints,[3 3],0,'both');
for idx1 = 1:size(endpoints, 1)
    for idx2 = 1:size(endpoints, 2)
     if(endpoints(idx1, idx2) == 1)
         for i = -3:3
             for j = -3:3
                if(endpoints(idx1+i, idx2+j)==1 && i ~= 0 && j ~= 0)
                    endpoints(idx1,idx2) = 0;
                end
             end            
         end         
     end
    end
end
[real_ep_x, real_ep_y] = ind2sub([334 346], find(endpoints > 0));

branchpoints = padarray(branchpoints,[3 3],0,'both');
for idx1 = 1:size(branchpoints, 1)
    for idx2 = 1:size(branchpoints, 2)
     if(branchpoints(idx1, idx2) == 1)
         for i = -3:3
             for j = -3:3
                if(branchpoints(idx1+i, idx2+j)==1 && i ~= 0 && j ~= 0)
                    branchpoints(idx1,idx2) = 0;
                end
             end            
         end         
     end
    end
end
[real_bp_x, real_bp_y] = ind2sub([334 346], find(branchpoints > 0));

figure(6);
subplot(1,3,1);
imshow(skeltonized);
title('skeletonized image');

subplot(1,3,2);
imshow(cleaned);
hold on;
plot(ep_y, ep_x, 'ro', bp_y, bp_x, 'gs');
title('cleaned from small independent line segments');
legend('end-points', 'branch-points');

subplot(1,3,3);
imshow(cleaned);
hold on;
plot(real_ep_y, real_ep_x, 'ro', real_bp_y, real_bp_x, 'gs');
title('only the real end- and br-points')
legend('real end-points', 'real branch-points');





