function max_pic = find_local_max(pic,h,true_max)
% Function by Jakob Lønborg Christensen s183985 F2020
%
% function that finds local maximum from strel h
%
% returns a logical matrice where 1=maxima and 0=not maxima
if nargin < 3
    true_max = 1;
end
if nargin == 1
    if size(pic,3) == 1
        h = ones(3,3);
    else
        h = ones(3,3,3);
    end
end
if true_max
    h((numel(h)+1)/2) = 0;
    max_pic = imdilate(pic,strel(h))<pic;
else
    max_pic = imdilate(pic,strel(h))==pic;
end
end




function max_pic = find_local_min(pic,h,true_max)
% Function by Jakob Lønborg Christensen s183985 F2020
%
% function that finds local minima from strel h
%
% returns a logical matrice where 1=maxima and 0=not maxima
pic = -pic;
if nargin < 3
    true_max = 1;
end
if nargin == 1
    if size(pic,3) == 1
        h = ones(3,3);
    else
        h = ones(3,3,3);
    end
end
if true_max
    h((numel(h)+1)/2) = 0;
    max_pic = imdilate(pic,strel(h))<pic;
else
    max_pic = imdilate(pic,strel(h))==pic;
end
end

function [centers,r_vec] = find_voxel_blobs(voxel_pic,t_vec,num_blobs,find_min,find_max)
if nargin == 3
    find_min = 1;
    find_max = 1;
end
if find_min && find_max
    voxel_pic_edit = abs(voxel_pic);
elseif find_min
    voxel_pic_edit = voxel_pic;
elseif find_max
    voxel_pic_edit = -voxel_pic;
else
    voxel_pic_edit = voxel_pic;
end
max_pic = find_local_max(voxel_pic_edit);
threshold = 0.4;
if num_blobs < 1
    num_max = sum(max_pic,'all');
    avg_val = sum(double(max_pic).*voxel_pic_edit,'all')/num_max;
    thresh_val = avg_val*threshold;
    max_pic = (double(max_pic).*voxel_pic_edit>thresh_val);
    num_blobs = sum(max_pic,'all');
end
centers = zeros(num_blobs,2);
voxel_pic_edit = voxel_pic_edit.*double(find_local_max(voxel_pic_edit));
voxel_pic_vec = reshape(voxel_pic_edit,numel(voxel_pic_edit),1);
[B,I] = maxk(voxel_pic_vec,num_blobs);
[centers(:,1),centers(:,2),rad_idx] = ind2sub(size(voxel_pic),I);
r_vec = sqrt(2*t_vec(rad_idx));
end

function g = get_gauss_filt(t,order_deriv,size)
% Function by Jakob Lønborg Christensen s183985 F2020
%
% Computes a 1d gaussian filter.
%
% Inputs:
%           t           = positive real numbers. The stan. dev. of the filter
%           order_deriv = {0,1,2}. Number of times the derivative of the
%                         gaussian is taken. default: 0
%           size        = Naturals. Makes the returned vector 2*size+1 long.
%                         default: around 3*x sd
% Outputs:
%           g           = vector of size 1xN. 1d gaussian filter kernel.
%
if nargin == 1
    order_deriv = 0;
end
if nargin < 3
    size = ceil(3*t);
end
x = -size:size;
g = 1/sqrt(2*pi*t^2)*exp(-x.^2/(2*t));
if order_deriv == 1
    g = g.*(-x/t);
elseif order_deriv == 2
    g = g.*(x.^2-t)/t^2;
end
end

function [voxel_pic] = make_blob_voxel(pic,t_vec)
voxel_pic = zeros([size(pic),length(t_vec)]);
for i=1:length(t_vec)
    g2 = get_gauss_filt(t_vec(i),2);
    g0 = get_gauss_filt(t_vec(i),0);
    lap_1 = imfilter(pic,g2,'replicate');
    lap_1 = imfilter(lap_1,g0','replicate');
    lap_2 = imfilter(pic,g2','replicate');
    lap_2 = imfilter(lap_2,g0,'replicate');
    lap_full = lap_1 + lap_2;
    voxel_pic(:,:,i) = lap_full*t_vec(i)^2;
end
end

function pic_highlight = plot_blob_maxima(pic,t,num_blobs,lambda)
if nargin == 3
    lambda = 0.5;
end
r = sqrt(2*t);
centers = zeros(num_blobs,2);
Idx = zeros(1,2);

pic_edit = abs(pic);
pic_circle = (pic-min(pic,[],'all'))*((max(pic,[],'all')-min(pic,[],'all'))^-1);

h = fspecial('disk',r)/max(fspecial('disk',r),[],'all');
for i=1:num_blobs

[~,idx] = max(pic_edit,[],'all','linear');
[Idx(1),Idx(2)] = ind2sub(size(pic),idx);

pic_black = zeros(size(pic));
pic_black(Idx(1),Idx(2)) = 1;

pic_black = 1-imfilter(pic_black,h);
pic_circle = pic_circle.*pic_black + (1-pic_black)*(1-lambda).*pic_circle+...
    (1-pic_black)*lambda.*reshape([0,1,0.3],1,1,3);
centers(i,:) = Idx;
pic_edit = pic_edit.*pic_black;
end

imshow(pic_circle)
end

function pic_plot = plot_circles_on_pic(pic,centers,r_vec)
lambda = 0.7;
r_vec = round(r_vec);
n = size(centers,1);
if length(r_vec) == 1
    r_vec = repmat(r_vec,n,1);
end
color = reshape([0,1,0.3],1,1,3);
pic_plot = pic;
for i=1:n
    pic_black = zeros(size(pic));
    h_big = fspecial('disk',r_vec(i))/max(fspecial('disk',r_vec(i)),[],'all');
    h_small = fspecial('disk',r_vec(i)*0.9-1)/max(fspecial('disk',r_vec(i)*0.9-1),[],'all');
    pic_black(centers(i,1),centers(i,2)) = 1;
    pic_black = imfilter(pic_black,h_big) - imfilter(pic_black,h_small);
    pic_plot = pic_plot.*(1-pic_black) + pic_plot.*pic_black*(1-lambda) + pic_black*lambda.*color;
end
end




%% WEEK 2 advanced billed analyse
clc,clear
DirW2 = 'C:\Users\jakob\Desktop\DTU\Videregående Billedanalyse\Week2\EX_2_data\EX_2_data\';
%% ex 2.1.1 (1)
t = 10;
g0 = get_gauss_filt(t,0);
g1 = get_gauss_filt(t,1);
g2 = get_gauss_filt(t,2);
plot([g0;g1;g2]')
%% (2)

blob_unif = im2double(imread([DirW2,'test_blob_uniform.png']));
imshow(blob_unif)
%% (2) (3)
t_vec = 2.^(1:8);
figure
for i=1:8
    g2 = get_gauss_filt(t_vec(i),2);
    g0 = get_gauss_filt(t_vec(i),0);
    blob_unif_filt_1 = imfilter(blob_unif,g2);
    blob_unif_filt_1 = imfilter(blob_unif_filt_1,g0');
    blob_unif_filt_2 = imfilter(blob_unif,g2');
    blob_unif_filt_2 = imfilter(blob_unif_filt_2,g0);
    blob_unif_filt = blob_unif_filt_1 + blob_unif_filt_2;
    subplot(2,4,i), imagesc(blob_unif_filt), title(['t=',num2str(t_vec(i))])
    %max(blob_unif_filt,[],'all')
end
%% ex 2.1.2 (1) Detecting blobs on one scale


t_vec = 2.^((1:5)+6);
figure
for i=1:5
    g2 = get_gauss_filt(t_vec(i),2);
    g0 = get_gauss_filt(t_vec(i),0);
    blob_unif_filt_1 = imfilter(blob_unif,g2,'replicate');
    blob_unif_filt_1 = imfilter(blob_unif_filt_1,g0','replicate');
    blob_unif_filt_2 = imfilter(blob_unif,g2','replicate');
    blob_unif_filt_2 = imfilter(blob_unif_filt_2,g0,'replicate');
    blob_unif_filt = blob_unif_filt_1 + blob_unif_filt_2;
    subplot(2,5,i),imshow(blob_unif) , title(['t=',num2str(t_vec(i))])
    subplot(2,5,5+i), plot_blob_maxima(blob_unif_filt,t_vec(i),14)
end
%% ex 2.1.3 Detecting blobs on multiple scales
pic1 = blob_unif;
pic2 = im2double(imread([DirW2,'test_blob_varying.png']));
pic3 = im2double(imread([DirW2,'CT_lab_high_res.png']));
pic4 = pic3(300:end-300,300:end-300);
g0 = get_gauss_filt(2,0);
pic5 = imfilter(pic4,g0);
pic5 = imfilter(pic5,g0');
pic = pic2;
t_vec = (8:2:64).^2/2;
[voxel_pic] = make_blob_voxel(pic,t_vec);
[centers,r_vec] = find_voxel_blobs(voxel_pic,t_vec,16,1,1);
plot_circles_on_pic(pic,centers,r_vec)
%% ex 2.1.4 Detecting blobs in real data
pic = im2double(imread([DirW2,'CT_lab_high_res.png']));
pic = pic(300:end-300,300:end-300);
g0 = get_gauss_filt(2,0);
pic = imfilter(imfilter(pic,g0),g0');
%imshow(pic)

t_vec = (7:0.5:10).^2/2;
[voxel_pic] = make_blob_voxel(pic,t_vec);
[centers,r_vec] = find_voxel_blobs(voxel_pic,t_vec,200,0,1);
plot_circles_on_pic(pic,centers,r_vec)
%% ex 2.1.5 Localize blobs. show centers with gaussian filter
pic = im2double(imread([DirW2,'CT_lab_high_res.png']));
%pic2 = pic(300:(end-300),300:(end-300));
g0 = get_gauss_filt(12,0);
pic3 = imfilter(imfilter(pic,g0),g0');
col = reshape([0,0.3,1],1,1,3);
imshow(find_local_max(pic3).*col+pic*0.3*max(pic,[],'all')^-1)
%% show resulting circles with guassian filter
t_vec = (5:0.3:11).^2/2;
[voxel_pic] = make_blob_voxel(pic,t_vec);

[centers,r_vec] = find_voxel_blobs(voxel_pic.*find_local_max(pic3),t_vec,-1,0,1);

imshow(pic)
viscircles(centers(:,[2,1]), r_vec,'Color','g','LineWidth',0.3);
%% compare centers in gaussian filter (red) and imerode+gaussian filter (green)

g0 = get_gauss_filt(12,0);
pic4 = imerode(pic,fspecial('disk',3));
pic5 = imfilter(imfilter(pic4,g0),g0');
pic_test = pic*0.3.*reshape(ones(3,1),1,1,3);
pic_test(:,:,2) = pic_test(:,:,2) +  find_local_max(pic5);
pic_test(:,:,1) = pic_test(:,:,1) +  find_local_max(pic3);
imshow(pic_test)
%% show resulting circles with imerode+gaussian filter
t_vec = (5:0.3:11).^2/2;
[voxel_pic] = make_blob_voxel(pic,t_vec);

[centers,r_vec] = find_voxel_blobs(voxel_pic.*find_local_max(pic5),t_vec,-1,0,1);

imshow(pic)
viscircles(centers(:,[2,1]), r_vec,'Color','g','LineWidth',0.3);
%% load 4 best fiber pics
pic_names = {'CT_lab_med_res.png','CT_synchrotron.png','Optical.png','SEM.png'};
for i=1:length(pic_names)
    pic_list{i} = im2double(imread([DirW2,pic_names{i}]));
end
pic_list{4} = pic_list{4}(1:692,:);
pic_list{3}(1:43,1:171) = 0;
montage(pic_list)

%%
erode_vec = [2,3,3,3];
gauss_vec = [3,8,8,8];
r_vals = [3,5.5;
          3.3,7.5;
          3.5,8;
          5,9];
for i=4
name = pic_names{i};
picture = pic_list{i};
g0 = get_gauss_filt(gauss_vec(i),0);
picture4 = imerode(picture,fspecial('disk',erode_vec(i)));
picture5 = imfilter(imfilter(picture4,g0),g0');
t_vec = linspace(r_vals(i,1),r_vals(i,2),12).^2/2;
[voxel_pic] = make_blob_voxel(picture,t_vec);

[centers,r_vec] = find_voxel_blobs(voxel_pic.*find_local_max(picture5),t_vec,-1,0,1);

imshow(picture)
viscircles(centers(:,[2,1]), r_vec,'Color','g','LineWidth',0.3);
centers_list{i} = centers;
r_vec_list{i} = r_vec;
end
%%
subplot(2,2,1), imshow(picture4)
subplot(2,2,2), imshow(picture*0.3+find_local_max(picture5))
subplot(2,2,3), hist(r_vec,12)
subplot(2,2,4), imshow(picture5*max(picture5,[],'all')^-1)