Question about hair removal matlab code.

Hello all! I am doing a program about hair removal. And I found a code online but part of this code I couln't understand. The final part(the loop part) was really confused me.

img_b is a binary image and after dilation we could get the img_b binary image, which the black part is skin and the white part is hair. Thus, the 'for' part in the final loop is to decide if this pixel is skin or hair, if the pixel is in the skin part then it keeps the same pixel value as the original image; if the pixel is in the hair part, then do 'else' and calculate it for removing hair.

So here comes the question, I do not understand the calculation part in 'else' and the meaning of parameter r.When I change the value of r into a bigger one, the running speed gets slower. So I wonder could anyone give me a guide about this part please?when I change the value of r into a bigger one, the running speed gets slower. Thank you!

Images are listed below: the first image is which we need to deal with and the second image is the processed image.

Updated on Mar. 7th 2019: Sorry guys forgot to put the otsu function. Already put it in case anyone need it.

if true
clear;
clc;
r=7;
se1_para = 3;
se2_para = 2;
img_old = imread('D:\MATLAB1\bin\picsam\03-2.jpg');
figure(1),imshow(img_old);
[x,y,z] = size(img_old);
img = rgb2gray(img_old);
se1 = strel('disk',se1_para);
img_c = imclose(img,se1)
figure(2), imshow(img_c,[]);
img_fur = double(img_c) - double(img);
figure(3),imshow(img_fur,[]);
[X Y]=meshgrid(1:x);
tt=(X-280).^2+(Y-280).^2<280^2;
thresh = otsu(img_fur(tt),sum(tt(:)));
img_b = (img_fur>thresh);
figure(4),imshow(img_b);
se2 = strel('disk',se2_para);
img_b =imdilate(img_b,se2);
figure(5),imshow(img_b);
img_new = uint8(zeros(x,y,z));
for i = 1:x 
  for j = 1:y        
      if img_b(i,j) == false  
          img_new(i,j,:) = img_old(i,j,:);
      else           
          ttt = img_old(max(1,i-r):min(i+r,x),max(j-r,1):min(j+r,y),:);
          no_efficient_pix = cat(3,img_b(max(1,i-r):min(i+r,x),max(j-r,1):min(j+r,y)),not(tt(max(1,i-r):min(i+r,x),max(j-r,1):min(j+r,y))));
          no_efficient_pix = any(no_efficient_pix,3);
          ttt = ttt.*repmat(uint8(not(no_efficient_pix)),[1,1,3]);
          efficient_pix_num = (2*r+1)^2-sum(no_efficient_pix(:));
          img_new(i,j,:) = uint8(sum(sum(ttt))./efficient_pix_num);           
      end       
  end   
end
figure(6),imshow(img_new,[]);
end

function thresh = otsu(data,pix_num)
img_var = zeros(256,1);
for i=1:256
     w0 = sum(sum(data<=i-1))./pix_num;
     w1 = 1-w0;
     u0 = sum(sum(data.*double(data<=i-1)))./(w0*pix_num);
     u1 = sum(sum(data.*double(data>i-1)))./(w1*pix_num);
     img_var(i) = w0.*w1.*((u0-u1).^2);
end
[~,I] = max(img_var);
thresh = I-1;
end

[Image Analyst]

I suggest you contact the author. It's the author's fault for not commenting his code. SHAME on him!

[Amirreza Mahbod]

Hi Wei!

I can not run this code. I got error in this line: (otsu function is undefined)

thresh = otsu(img_fur(tt),sum(tt(:)));

I tried to replace it with:

thresh = graythresh(img_fur);

but I got another error in this line:(Index exceeds matrix dimensions error)

no_efficient_pix = cat(3,img_b(max(1,i-r):min(i+r,x),max(j-r,1):min(j+r,y)),not(tt(max(1,i-r):min(i+r,x),max(j-r,1):min(j+r,y))));

Any suggestion?

[Asmaa Mohamed]

Did you solve this error of undefined Otsu function ?? I have the same error!!

[SARAH LONER]

hai

Do you have this error?

no_efficient_pix = cat(3,img_b(max(1,i-r):min(i+r,x),max(j-r,1):min(j+r,y)),not(tt(max(1,i-r):min(i+r,x),max(j-r,1):min(j+r,y))));

Please help to solve this.

[Image Analyst]

J, attach your image and code so we can run it and fix it.

[SARAH LONER]

this are my images sir nd i have used the above code for removing hair from skin. please help me on on this

Answer 1

[Asghar ali]

clc;

close all;

clear all;

[y,Fs] = audioread('Bear.wav');

x = y(:,1);

clear y;

N = length(x);

time = (1:N).*(1/Fs);

han = plot(time,x);

xlab = xlabel('Seconds');

ylab = ylabel('Amplitude');

grid on;

set(han,'LineWidth', 2);

set([xlab, ylab],'FontSize', 24, 'FontName', 'Times');

set(gca,'FontSize',20,'FontName','Times','Fontweight','Bold')

% wavwrite(x,Fs,bits,'test.wav');

[X,f] = centeredFFT(x,Fs);

Y = fftshift(X,Fs);

figure;

han1 = plot(f,abs(X));

axis([-8000,8000,0,max(abs(X))]);

grid on;

xlab1=xlabel('Frequency(Hz)');

ylab1=ylabel('|X[k]|');

set(han1,'LineWidth', 2);

set([xlab1, ylab1],'FontSize', 24, 'FontName', 'Times');

set(gca,'FontSize',20,'FontName','Times','Fontweight','Bold')

% sound(x,Fs,bits);

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

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

fc = 2000;

fNorm = fc / ( Fs/2);

[b,a] = butter(10, fNorm, 'low');

x_Low = filtfilt(b, a, x);

figure;

freqz(b,a,bits,Fs/2);

grid on;

% sound(x_Low,Fs,bits);

% wavwrite(x_Low,Fs,bits,'x_Low_.wav');

[X_LOW,f_low] = centeredFFT(x_Low,Fs);

figure;

han2=plot(f_low,abs(X_LOW));

axis([-.8e4,.8e4,0,max(abs(X))]);

grid on;

xlab2=xlabel('Frequency(Hz)');

ylab2=ylabel('|X[k]|)');

set(han2,'LineWidth', 2);

set([xlab2, ylab2],'FontSize', 24, 'FontName', 'Times');

set(gca,'FontSize',20,'FontName','Times','Fontweight','Bold')

figure;

han3=plot(time,x_Low);

grid on;

xlab3=xlabel('Frequency(Hz)');

ylab3=ylabel('|X[k]|)');

set(han3,'LineWidth', 2);

set([xlab3, ylab3],'FontSize', 24, 'FontName', 'Times');

set(gca,'FontSize',20,'FontName','Times','Fontweight','Bold')

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

%%%%%%%%%THIS IS A COOL CUSTOM FFT PLOTTING FUNCTION BY MATHWORKS%%%%%%%%%%

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

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

% % function [X,freq]=centeredFFT(x,Fs)

% % %this is a custom function that helps in plotting the two-sided spectrum

% % %x is the signal that is to be transformed

% % %Fs is the sampling rate

% %

% % N=length(x);

% %

% % %this part of the code generates that frequency axis

% % if mod(N,2)==0

% % k=-N/2:N/2-1; % N even

% % else

% % k=-(N-1)/2:(N-1)/2; % N odd

% % end

% % T=N/Fs;

% % freq=k/T; %the frequency axis

% %

% % %takes the fft of the signal, and adjusts the amplitude accordingly

% % X=fft(x)/N; % normalize the data

% % X=fftshift(X); %shifts the fft data so that it is centered

kindly remove error from above code plz do fast

Answer 2

[Luis Enciso]

I fixed the given code in this way, it works fine in squared cropped images:

close all
clear;
clc;
r=7;
se1_para = 3;
se2_para = 2;
img_old = imread('image.jpg');
figure(1),imshow(img_old);

% crop
[x,y,z] = size(img_old);
d = y-x;
im_crop = img_old(:,round(d/2)+1:y-(d-round(d/2)),:);

[x,y,z] = size(im_crop);
img = rgb2gray(im_crop);
se1 = strel('disk',se1_para);
img_c = imclose(img,se1)
figure(2), imshow(img_c,[]);
img_fur = double(img_c) - double(img);
figure(3),imshow(img_fur,[]);
[X Y]=meshgrid(1:x);
tt=(X-x/2).^2+(Y-y/2).^2<(x+y)^2;
thresh = otsu(img_fur(tt),sum(tt(:)));
img_b = (img_fur>thresh);
figure(4),imshow(img_b);
se2 = strel('disk',se2_para);
img_b =imdilate(img_b,se2);
figure(5),imshow(img_b);
img_new = uint8(zeros(x,y,z));
for i = 1:x 
  for j = 1:y        
      if img_b(i,j) == false  
          img_new(i,j,:) = im_crop(i,j,:);
      else           
          ttt = im_crop(max(1,i-r):min(i+r,x),max(j-r,1):min(j+r,y),:); % region 
          no_efficient_pix = cat(3,img_b(max(1,i-r):min(i+r,x),max(j-r,1):min(j+r,y)),not(tt(max(1,i-r):min(i+r,x),max(j-r,1):min(j+r,y))));
          no_efficient_pix = any(no_efficient_pix,3);
          ttt = ttt.*repmat(uint8(not(no_efficient_pix)),[1,1,3]);
          efficient_pix_num = (2*r+1)^2-sum(no_efficient_pix(:));
          img_new(i,j,:) = uint8(sum(sum(ttt))./efficient_pix_num);           
      end       
  end   
end
figure(6),imshow(img_new,[]);

function thresh = otsu(data,pix_num)
img_var = zeros(256,1);
for i=1:256
     w0 = sum(sum(data<=i-1))./pix_num;
     w1 = 1-w0;
     u0 = sum(sum(data.*double(data<=i-1)))./(w0*pix_num);
     u1 = sum(sum(data.*double(data>i-1)))./(w1*pix_num);
     img_var(i) = w0.*w1.*((u0-u1).^2);
end
[~,I] = max(img_var);
thresh = I-1;
end