hair removal

clear all; close all; clc;
%% Input images information
%addpath('C:\Users\TBL\Documents\MATLAB\Wrinkle Quantification\Skin_Demo\Single\Wrinkle_Detection\')
%img = imread('bb _333.jpg');

%% Negative degree
%addpath('C:\Users\TBL\Documents\MATLAB\1708_AMORE\02_YJ\AF_backhand')
%addpath('C:\Users\TBL\Documents\MATLAB\1708_AMORE\01_KB\AF_forehead')
%addpath('C:\Users\TBL\Documents\MATLAB\1708_AMORE\02_YJ\AF_forehead')
%addpath('C:\Users\TBL\Documents\MATLAB\1708_AMORE\02_YJ\AF_forehead')
%% zero degree
%addpath('C:\Users\TBL\Documents\MATLAB\1708_AMORE\01_KB\AF_backhand')
%addpath('C:\Users\TBL\Documents\MATLAB\1708_AMORE\01_KB\AF_backhand')
%addpath('C:\Users\TBL\Documents\MATLAB\1708_AMORE\10_MR\AF_forearm')
%addpath('C:\Users\TBL\Documents\MATLAB\1708_AMORE\08_MA\   forearm')
%% Positive degree
%addpath('C:\Users\TBL\Documents\MATLAB\1708_AMORE\06_HG\AF_forehead') %Net
%volosa
%addpath('C:\Users\TBL\Documents\MATLAB\1708_AMORE\04_HG\AF_forehead')
%addpath('C:\Users\TBL\Documents\MATLAB\1708_AMORE\03_PU\AF_forehead')
addpath('C:\Users\TBL\Documents\MATLAB\1708_AMORE\05_HR\AF_forehead')  %Volosy
%addpath('C:\Users\TBL\Documents\MATLAB\1708_AMORE\06_HG\AF_forehead')


%addpath('C:\Users\TBL\Documents\MATLAB\1708_AMORE\06_HG\forehead');


%addpath('C:\Users\TBL\Documents\MATLAB\1708_AMORE\07_EJ\AF_forearm')
%addpath('C:\Users\TBL\Documents\MATLAB\1708_AMORE\10_MR\forehead')


img = imread('Image0499.jpg');
%img = imread('Image0499.jpg');

%img = rgb2gray(img);
%imshow(img);
[H,W] = size(img);

%% Parameters setup
[lm,li] = max(img(5:H,51:W));
minDepth = min(li)
maxDepth = max(li)
%%

diff = 20; %it was 2
l_crop = 50;
r_crop = 1;
t_crop = 14+minDepth-0.75.*minDepth;
b_crop = 0.75*(H-maxDepth);
%b_crop = 2*minDepth;
%%
up = H-2;%443
gsigma = 4;
down = 222;
% %% Preprocessing
[img_crop,H_crop,W_crop] = crop(img,H,W,l_crop,r_crop,t_crop,b_crop);
% imshow(img_crop)
% imshow(img_crop)
%%
figure()
imshow(img_crop);
%%

%%
[lm,li] = max(img_crop(:,1))
[rm,ri] = max(img_crop(:,W_crop))
shif = ri-li;
rotDeg = atand(shif/W_crop);
if abs(rotDeg) > 1
    img_crop = imrotate(img_crop,rotDeg);
else
    rotDeg = 0;
end
[H,W] = size(img_crop);
%%
figure()
imshow(img_crop);
%%
[row,col] = find(img_crop);
maxRow = max(row);
minRow = min(row);

if rotDeg > 0
    % Case this is rotated to left side
    colLeft = find(img_crop(maxRow,:));
    cLeft = max(colLeft);

    rowLeft = find(img_crop(:,cLeft));
    rLeft = min(rowLeft);

    %[rLeftM,rLeftI] = max(img_crop(:,cLeft));%-30

    rLeftI = rLeft-(shif-160); %-60
    %rLeftI = rLeft-min(20,(shif-100)); %-60
    [img_crop,H_crop,W_crop] = crop(img_crop,H,W,cLeft,cLeft,1,rLeftI);
elseif rotDeg < 0
    colRight = find(img_crop(maxRow,:));
    cRight = min(colRight);
    cRight = W-cRight;
    rowRight = find(img_crop(:,cRight));
    rRight = min(rowRight);
    %[rRightM,rRightI] = max(img_crop(:,cRight));%-30
    rRight = rRight+abs(shif); %+40
    [img_crop,H_crop,W_crop] = crop(img_crop,H,W,cRight,cRight,rRight,rRight);
else
end
%rLeft = 144;


% cRight = cLeft;
% rTop = rLeft;
% rBottom = rTop;


%cRight = W-cLeft
%%

%rRight = H-rLeft

%cLeft = 152;
%rLeft = 126;


 %left right top bot
%colRight = find(img_crop(minRow,:));
%%


img_filt = imgaussfilt(img_crop,gsigma);
img_filt = imgaussfilt(img_crop,gsigma);

%%

%%
diff = 1;
img_diff = filt_diff(img_filt,H_crop,W_crop,diff);
up = H_crop-2;
[img_line_surf,ind_line_surf,avg_index,ref_line,line1] = line_surf3(img_diff,H_crop,W_crop,diff,up,down);


for i = 1:W_crop
  img_crop(1:line1(i),i) = 0;
end

diff = 1;
img_diff = filt_diff(img_filt,H_crop,W_crop,diff);
up = H_crop-2;
[img_line_surf,ind_line_surf,avg_index,ref_line,line1] = line_surf3(img_diff,H_crop,W_crop,diff,up,down);



% img_med_surf = medfilt1(line1);
% figure()
% imshow(img_crop);
% hold on
% plot(img_med_surf)
% plot(line1)
% hold off


%imshow(img_med_surf,[])


img_line_surf_med = medfilt2(img_crop,[4 4]);

img_diff = filt_diff(img_line_surf_med,H_crop,W_crop,diff);
[img_line_surf,ind_line_surf,avg_index,ref_line,line1] = line_surf3(img_diff,H_crop,W_crop,diff,up,down);

for i = 1:W_crop
  img_line_surf(1:line1(i)-1,i) = 255;
end

img_line_surf = medfilt2(img_line_surf,[5 11]);
%img_line_surf = medfilt2(img_line_surf);

figure()
imshow(img_line_surf);






%{




%%
% generate horizontal edge emphasis kernel
%h = fspecial('sobel',10);
% invert kernel to detect vertical edges
%h = h';
%h = [-1,0,0,0,1;-1,0,0,0,1;-1,0,0,0,1;-1,0,0,0,1;-1,0,0,0,1;-1,0,0,0,1]
%-2,0,2;-1,0,1];
%img_crop = imfilter(img_crop,h);
%%





% Project (sum) image vertically to get the horizontal profile.
horizontalProfile = sum(img_crop, 1);
horizontalProfile = horizontalProfile.^3;
plot(horizontalProfile);

threshold1 = mean(horizontalProfile)+2.5*std(horizontalProfile);
%threshold = 30e+11;

for i = 1:(W_crop-1)
    lineis = (horizontalProfile(i) > threshold1);
    if horizontalProfile(i) > threshold1
        line1(i) = line1(i-5);
    end
end


figure()
imshow(img_crop,[])
hold on
plot(line1)
hold off


%sum(horizontalProfile)
%{
figure()
imshow(img_crop,[])
hold on
plot(line1)
hold off
%}

%}


% [Gx, Gy] = imgradientxy(img_crop,'prewitt');
% figure()
% imshowpair(Gx, Gy, 'montage');
%  title('Directional Gradients: x-direction, Gx (left), y-direction, Gy (right), using Prewitt method')
% axis off;
%
% meanGx = mean(mean(abs(Gy)))


%
% for i = 1:H_crop
%     for j=1:W_crop
%         if Gx(i,j) > mean(Gy)
%             img_crop(i,j) = 0;
%         end
%     end
% end



%}








%{


%%
%img_crop = double(img_crop);
img_filt = imgaussfilt(img_crop,gsigma);
%%

%%
diff = 1;
img_diff = filt_diff(img_filt,H_crop,W_crop,diff);
up = H_crop-2;
[img_line_surf,ind_line_surf,avg_index,ref_line,line1] = line_surf3(img_diff,H_crop,W_crop,diff,up,down);



for i = 1:W_crop
  img_crop(1:line1(i),i) = 0;
end

img_filt = imgaussfilt(img_crop,gsigma);


[H,W] = size(img_crop);
diff = 2;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

med = 80;
%img_med = imgaussfilt(img_crop,gsigma);
img_med = padarray(img_filt',med,'replicate','both');
img_med = padarray(img_med',med,'replicate','both');
img_med = medfilt2(img_med,[med*2 med*2]);
img_med = img_med([med+1:H+med],[med+1:W+med]);
img_filt1 = zeros(H+diff,W);
img_filt2 = zeros(H+diff,W);

img_filt1([diff+1:H+diff],:) = img_med;
img_filt2([1:H],[1:W]) = img_med;

%%
img_med = img_filt1 - img_filt2;
img_med = img_med([diff+1:H+diff],:);
%%

%%
% figure();
% imshow(img_med,[]);
[H,W] = size(img_crop);
img_line2 = zeros(H,W);
%diff = 10;
for i = 1:W
    [local_max,local_index] = min(img_med(:,i));
    img_line2(local_index,i) = 255; %+2 %-10    ubrat -diff
    line2(i) = local_index; %brat -diff
end
%%
% figure()
% imshow(img_line2,[])
%%

figure(1)
imshow(img_crop,[])
hold on
plot(line1);
plot(line2);
hold off


%%
figure(2)
imshow(img_crop,[])
hold on
plot(line2);
hold off



%% Roughness

[peak_height,valley_depth] = peak_valley(line1,line2,W);
[avg_rough,rms_rough,max_rough,total_rough] = roughness(peak_height,valley_depth);


%% Visualization


%}















































%someIm = img_crop(maxRow,:);
%[row,col] = find(img_crop());
% ind1 = img_crop(maxRow,1) != 0;
% %ind1 = A(maxRow,:) == 1;
% A1 = img_crop(ind1,:);
%}


%%
%{
figure()
imshow(img_crop);
figure()
imshow(rotImg);
%{
%img_crop = double(img_crop);
img_filt = imgaussfilt(img_crop,8);
%%
%imshow(img_filt)

%%
% %% Wrinkle surface
img_diff = filt_diff(img_filt,H_crop,W_crop,diff);
[img_line_surf,ind_line_surf,avg_index,ref_line,line1] = line_surf3(img_diff,H_crop,W_crop,diff,up,down);

%%
[lm,li] = max(img_crop(:,1))
[rm,ri] = max(img_crop(:,W_crop))

rotDeg = atand((ri-li)/W_crop);
rotImg = imrotate(img_crop,rotDeg);
%%

%%
%{
figure()
imshow(img_crop);
figure()
imshow(rotImg);
%}
%%

l_crop2 = 20;
r_crop2 = 20;
t_crop2 = 200;
b_crop2 = 10;

% %% Preprocessing
[img_crop2,H_2,W_2] = crop(rotImg,H_crop,W_crop,l_crop2,r_crop2,t_crop2,b_crop2);

%%
figure()
imshow(img_crop2);

%%


%}





























































































%%

%figure()
%imshow(img_filt,[]);
%figure();
%imshow(im,[]);


% figure();
% imshow(img_filt);
% figure()
% plot(I)
% ylim([0 450])



%
% % figure()
% % imshow(img_filt)
%
%
% % for i = 1:W
% %     [M,I] = max(img_filt(:,i));
% % end
% [M,I] = max(img_filt);
% minRow = min(I);
%
% for i = 1:W
%     diffShift(i) = I(i)-minRow;
%     lowLine1(i) = I(i)+100-diffShift(i);
% end
%
%
%
% for i = 1:W
%     img_filt_new([minRow:lowLine1(i)],i) = 255;
%     %new_img_filt(:,:) = img_filt((I(i):I(i)+200),:);
% end
%
% % imshow(img_filt_new);
%
%
% im=zeros(443,970);
% a=200;
% b=100;
% c = 200;
% for i=1:size(im,2)
%     j=443-I(i);
%     im(a+1-c:(a+j-9),i)=img_filt(I(i)-c:(I(i)-10+j),i);
% end
%
% figure();
% imshow(im,[]);
%
% figure();
% imshow(img_filt);
% figure()
% plot(I)
% ylim([0 450])











% im=zeros(1000,970);
% a=100;
% b=100;
% for i=1:size(im,2)
%     j=443-I(i);
%     im(a+1:(a+j-9),i)=img_filt(I(i)-b:(I(i)-10+j-b),i);
% end
% figure();
% imshow(im);
%
% figure();
% imshow(img_filt);










% %% Wrinkle surface
% img_diff = filt_diff(img_filt,H_crop,W_crop,diff);
% [img_line_surf,ind_line_surf,avg_index,ref_line,line1] = line_surf3(img_diff,H_crop,W_crop,diff,up,down);
%
% %%
%
% med = 80;
% img_med = imgaussfilt(img_crop,gsigma);
% img_med = padarray(img_med',med,'replicate','both');
% img_med = padarray(img_med',med,'replicate','both');
% img_med = medfilt2(img_med,[med*2 med*2]);
% img_med = img_med([med+1:H+med],[med+1:W+med]);
% img_filt1 = zeros(H+diff,W);
% img_filt2 = zeros(H+diff,W);
%
% img_filt1([diff+1:H+diff],:) = img_med;
% img_filt2([1:H],[1:W]) = img_med;
%
% img_med = img_filt1 - img_filt2;
% img_med = img_med([diff+1:H+diff],:);
% %%
%
% % figure();
% % imshow(img_med,[]);
%
% img_line2 = zeros(H,W);
% %diff = 10;
% for i = 1:W
%     [local_max,local_index] = max(img_med(:,i));
%     img_line2(local_index-diff,i) = 255; %+2 %-10
%     line2(i) = local_index-diff;
% end
% %%
%
%
% for i = 1: W
%     %if(line2(i)-line1(i) < 0)
%         %img_wrinkle([line2(i):line1(i)],i,1) = 255;
%         %img_wrinkle_only([line2(i):line1(i)],i,1) = 255;
%     depth(i) = diff+line1(i) - line2(i);
%     %end
% end
%
% maxDepth = max(depth);
% colormap jet
% cm = colormap % returns the current color map
% %cm(1, :) = [1 1 1];
% %colormap(cm)
% %%
% for i = 1: W
%     normDepth(i) = depth(i)/maxDepth;
%     colorID = max(1, sum(normDepth(i) > [0:1/length(cm(:,1)):1]));
%     myColor = cm(colorID, :); % returns your color
%     depthColor(i,:,:,:) = myColor;
% end
% img_wrinkle(:,:,1) = zeros(size(img_crop));
% img_wrinkle(:,:,2) = zeros(size(img_crop));
% img_wrinkle(:,:,3) = zeros(size(img_crop));
%
% %img_wrinkle_only = zeros(size(img_crop));
% %%
% for i = 1: W
%     %if(depth(i) > 0)  %(line2(i)-line1(i) < 0)
%         %img_wrinkle([line2(i):line1(i)],i,1) = 255;
%     img_wrinkle([line2(i):line1(i)+diff],i,1) = depthColor(i,1);
%     img_wrinkle([line2(i):line1(i)+diff],i,2) = depthColor(i,2);
%     img_wrinkle([line2(i):line1(i)+diff],i,3) = depthColor(i,3);
%
% end
%
% minLine = min(line2);
%
% for i = 1:W
%     diffShift(i) = line2(i)-minLine;
%     newLine1(i) = line1(i)-diffShift(i);
% end
% maxLine = max(newLine1);
%
% for i = 1: W
%     %if(depth(i) > 0)  %(line2(i)-line1(i) < 0)
%         %img_wrinkle([line2(i):line1(i)],i,1) = 255;
%     img_wrinkle_new([minLine:newLine1(i)+diff],i,1) = depthColor(i,1);
%     img_wrinkle_new([minLine:newLine1(i)+diff],i,2) = depthColor(i,2);
%     img_wrinkle_new([minLine:newLine1(i)+diff],i,3) = depthColor(i,3);
% %     img_wrinkle_new([line2(i):line1(i)+diff],i,2) = depthColor(i,2);
% %     img_wrinkle_new([line2(i):line1(i)+diff],i,3) = depthColor(i,3);
% end
% figure()
% imshow(img_wrinkle_new(minLine:maxLine+diff,:,:,:))
% figure()
% imshow(img_wrinkle)
% hold on
% plot(line2)
% plot(line1+diff)











% [row,col,v] = find(img_filt);

%%
% %
% minLine = min(row);
%
% for i = 1:W
%     for j = 1:H
%         diffShift(i,j) = img_filt(:,:)-minLine;
%     end
% end
%
% figure()
% imshow(diffShift)











% for i = 1:W
%     diffShift(i) = line2(i)-minLine;
%     newLine1(i) = line1(i)-diffShift(i);
% end
% maxLine = max(newLine1);
%
% for i = 1: W
%     if(depth(i) > 0)  %(line2(i)-line1(i) < 0)
%         img_wrinkle([line2(i):line1(i)],i,1) = 255;
%     img_wrinkle_new([minLine:newLine1(i)+diff],i,1) = depthColor(i,1);
%     img_wrinkle_new([minLine:newLine1(i)+diff],i,2) = depthColor(i,2);
%     img_wrinkle_new([minLine:newLine1(i)+diff],i,3) = depthColor(i,3);
%     img_wrinkle_new([line2(i):line1(i)+diff],i,2) = depthColor(i,2);
%     img_wrinkle_new([line2(i):line1(i)+diff],i,3) = depthColor(i,3);
% end
% figure()
% imshow(img_wrinkle_new(minLine:maxLine+diff,:,:,:))











% maxLine = max(newLine1);
%
% for i = 1: W
%     %if(depth(i) > 0)  %(line2(i)-line1(i) < 0)
%         %img_wrinkle([line2(i):line1(i)],i,1) = 255;
%     img_wrinkle_new([minLine:newLine1(i)+diff],i,1) = depthColor(i,1);
%     img_wrinkle_new([minLine:newLine1(i)+diff],i,2) = depthColor(i,2);
%     img_wrinkle_new([minLine:newLine1(i)+diff],i,3) = depthColor(i,3);
% %     img_wrinkle_new([line2(i):line1(i)+diff],i,2) = depthColor(i,2);
% %     img_wrinkle_new([line2(i):line1(i)+diff],i,3) = depthColor(i,3);
% end





%img_wrinkle_new = imcrop(img_wrinkle_new,[60 40 100 90]);
%%
% figure()

% hold on
% plot(newLine2)
% plot(mewLine1+diff)
%


%%


%% Roughness

% [peak_height,valley_depth] = peak_valley(ind_line_surf,avg_index,W_crop);
% [avg_rough,rms_rough,max_rough,total_rough] = roughness(peak_height,valley_depth);
%% Visualization
%% Wrinkle
% r1 = img_crop;
% g1 = r1;
% b1 = r1;
% r1(img_line_surf > 0) = 255;
% g1(ref_line > 0) = 255; %%
% b1(img_line_surf > 0) = 0;
% figure('Name','Wrinkle')
% imshow(cat(3,r1,g1,b1),[]);









% figure();
% imshow(img_med,[]);

%diff = 10;
%end

%img_wrinkle(:,:,1) = img_crop;
%img_wrinkle(:,:,2) = img_crop;
%img_wrinkle(:,:,3) = img_crop;


%img_wrinkle([line1(i):H],i,1) = 1;
%img_wrinkle([line1(i):H],i,2) = 1;
%img_wrinkle([line1(i):H],i,3) = 1;

% figure()
% imshow(line2)

%figure(1)
%imshow(img_line2)
%figure(2)
%imshow(img_line_surf)


%figure(1)
%imshow(img_line_surf)
%figure(2)
%imshow(img_line2)

%img_wrinkle(:,:,1) = img_crop;
%img_wrinkle(:,:,2) = img_crop;
%img_wrinkle(:,:,3) = img_crop;
%img_wrinkle_only = zeros(size(img_crop));
%depth = zeros(W);


%figure()


% figure()
% imshow(img_line_surf)
% figure()
% imshow(img_crop)
% hold on;
% plot(line2)













%colormap gray
% for i = 1: Wz
%     %if(line2(i)-line1(i) < 0)
%         %img_wrinkle([line2(i):line1(i)],i,1) = 255;
%     img_wrinkle([line1(i):H],i,1) = 1;
%     img_wrinkle([line1(i):H],i,2) = 1;
%     img_wrinkle([line1(i):H],i,3) = 1;
%     %end
% end

% r1 = img_crop;
% g1 = r1;
% b1 = r1;
% r1(img_line_surf > 0) = 255;
% g1(ref_line > 0) = 255; %%
% b1(img_line_surf > 0) = 0;
% figure('Name','Wrinkle')
% imshow(cat(3,r1,g1,b1),[]);
%f = 0.568; % your float

%cm = colormap % returns the current color map

%colorID = max(1, sum(f > [0:1/length(cm(:,1)):1]));

%myColor = cm(colorID, :) % returns your color



%figure(1)
%imshow(img_wrinkle,[])

%figure(4)
%imshow(img_wrinkle_only,[])

%}