function [CAN] = detect(name_im,color)
%
%
%      This program is free software; you can redistribute it and/or modify
%      it under the terms of the GNU General Public License as published by
%      the Free Software Foundation. This program is distributed as 
%      test tool and is not optimized. 
%      For any problem or question, please contact:
%      bruni@iac.rm.cnr.it
%      
%
%
% function [CAN] = detect(name_im,color)
%
% Input variables
%     name_im = scratched image
%     color is 1 for detecting black scratches whilt it is 0 for
%           white scratches.
%
%
% Output variables
%     CAN = vector containing detected scratches locations (column number)   
%
% 
%
A = imread(name_im);

tic

[m,n] = size(A);
heig = round(m/3);
step = 10;
CAN = [];
row = 1;
for i = 1:1
    fasc = row + heig - 1; 
    if fasc > m
       fasc = m; 
    end 
    Im = A(row:fasc,:);
    [can] = select3(Im,step,color);
    toc   
    CAN(i,1:length(can)) = can;
    if (size(CAN,2))  ~= length(can)
       CAN(i,length(can) + 1:size(CAN,2)) = 0;
    end
    row = fasc + 1; 
end

% plot
mat = 255 * ones(20,n);
mat(:,can)=0;
mat = [mat ; double(A)/256];

figure
imshow(mat)
title('detected scratches')





%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%% S U B R O U T I N E S %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


function [can] = select3(I1,step,color)
%
% Input variables
%       I1 = scratched image
%       step = parameter for cross-section computation
%       colore = 1 for black scratches 
%                0 for white scratches
%
% Output variables
%       can = vector containing detected scratches locations (column number) 
% 


% image luminance cross-section computation
cross = cross_sec(I1,step);

% width and heigth selections
[enmed,enmedlat,en,mxm,primicand,bmed] = select2(cross,color);

% energy selection (Weber's law)
lc = length(primicand);
n = size(I1,2);
candid = [];
cont = 0;
for i = 1:lc
    cp = primicand(i);
    [ok] = sel_web(cross,mxm,cp,color,enmed,enmedlat,en,bmed);
    if ok == 1
       candid = [candid cp];
       cont = cont + 1;
    end
end
can = mxm(candid);


%
% plot
%
figure
plot(1:n,cross,'b-',mxm(candid),cross(mxm(candid)),'r.')
axis([1 n -40 40])
set(gca,'XTick',20:50:length(cross)-5)
set(gca,'YTick',-30:30:30)
xlabel('column number')
ylabel('gray level')
title('detected scratches locations')






function [a] = cross_sec(I,s)
%
% image luminance cross-section computation
%
% Input variables
%       I   = scratched image
%       s   = step parameter for local mean computation
%       
%
% Output variables
%       a   = cross-section of image I       


sig = mean(I);
lsig = length(sig);
h = 0;
for i = 1:lsig
    lims = i + s; 
    if (i + s) > lsig
        lims = lsig;
    end  
    a(i) = sig(i) - mean(sig(i - h:lims));
    if i <= s
       h = h + 1;
    end
end
  




function [enmed,enmedlat,en,j,jind,bmed] = select2(cross,color)
%
% width and height selections
%
% 
% 
%

x = diff(cross);
l = length(x);
y = sign(x);
if y(1) == 0
   k = 1;
else
  j(1) = 1;
  k = 2;
end
for h = 1:l - 1
    if (y(h + 1) == 0)
        y(h + 1) = y(h);
    end            
    if (y(h) ~= y(h + 1))  
       j(k) = h + 1;
       k = k + 1;
    end
end  
lj = length(j);

%
% if color = 1 minima of cross-section are computed
% if color = 0 maxima of cross-section are computed
%

if color == 1
   if cross(j(2)) > cross(j(1))
      m = 3:2:lj;
   else
      m = 4:2:lj;
   end  
else 
   if cross(j(2)) < cross(j(1))
      m = 3:2:lj;
   else
      m = 4:2:lj; 
   end    
end
lm = length(m);
k = 0; 
k1 = 0;
if m(lm) == lj
   lm = lm - 1;
end
bmed = mean(cross(j(m)));
for h = 1:lm
    if (((j(m(h) + 1) - j(m(h) - 1)) >= 3) & ((j(m(h) + 1) - j(m(h) - 1)) <= 10))
       k = k + 1;  
       indc(k) = m(h);
       if k1 == 0
          indc1(1) = m(h) - 1; 
          k1 = 2;
       elseif indc1(k1 - 1)  ~= (m(h) - 1)
              indc1(k1) = m(h) - 1;
              k1 = k1 + 1;
        end
       indc1(k1) = m(h) + 1; 
       k1 = k1 + 1;     
    end
end
lindc = length(indc);
lindc1 = length(indc1);
en = zeros(1,length(cross));
f = lindc;
if indc(lindc) == lj 
   f = lindc - 1;
end
for h = 1:f
    en(j(indc(h))) = ((abs(cross(j(indc(h))) - cross(j(indc(h) - 1))) + abs(cross(j(indc(h) + 1)) - cross(j(indc(h)))))) / 2;  
end
enmed = mean(en(j(indc)));
ind = find(en(j(indc)) >= enmed);
jind = indc(ind);
f = lindc1;
if indc1(lindc1) == lj 
   f = lindc1 - 1;
end
for h = 1:f
    en(j(indc1(h))) = ((abs(cross(j(indc1(h))) - cross(j(indc1(h) - 1))) + abs(cross(j(indc1(h) + 1)) - cross(j(indc1(h)))))) / 2;  
end
enmedlat = mean(en(j(indc1)));




function [ok] = sel_web(x,mxm,cp,color,enmed,enmedlat,en,bmed)
%
%  energy selection based on Weber's law
%  
%  Output variables
%        ok = vector containing detected scratches locations
%

bp = x(mxm(cp));
lmxm = length(mxm);
if mxm(cp) <= 24
   mcps = mxm(1); 
else 
   mcps = mxm(cp) - 24;
end
if mxm(cp) > (mxm(lmxm) - 24)
   mcpd = mxm(lmxm);
else
   mcpd = mxm(cp) + 24;
end
trova = find(mxm(1:cp) > mcps);
cps = trova(1);
trova1 = find(mxm(cp:lmxm) >= mcpd);
cpd = trova1(1) + cp - 1;


%
%
% bpw = mimimum value of bp satisfying Weber's law 
% t3 = energy threshold using Weber's law
%
%
bpw = bmed / 0.98;
M = (x(mxm(cp - 1)) + x(mxm(cp + 1))) / 2;
dislob = floor((mxm(cp + 1) - mxm(cp - 1)) / 2);
Erette = abs(bp - M) * dislob;
alfa = 10;
gamma = enmed / (alfa * abs(bp));
coeffsm = (abs(M) / abs(bp)) ^ (1 / dislob);
if coeffsm > 1
   coeffsm = 1;
end
t3 = abs(bpw - M) * dislob;



%
% errors computation
%
man1 = find(en(mxm(cp - 2:-2:cps)) <= enmed);
if length(man1) == 0
   if color == 1
      n1 = -enmed / 2;
   else
      n1 = enmed / 2;
   end 
else
   maxn1 = x(mxm(cp - 2 * man1));
   if color == 1
      n1 = min(maxn1);
   else
      n1 = max(maxn1);
   end
end
man2 = find(en(mxm(cp - 3:-2:cps)) <= enmedlat);
if length(man2) == 0
   if color == 1
      n2 = enmedlat / 2;
   else
      n2 = -enmedlat / 2;
   end 
else
   maxn2 = x(mxm(cp - 2 * man2 - 1));
   if color == 1
      n2 = max(maxn2); 
   else
      n2 = min(maxn2);
   end
end
vl = abs(n1 - n2) / 6;
mand1 = find(en(mxm(cp + 2:2:cpd)) <= enmed);
if length(mand1) == 0
   if color == 1
      nd1 = -enmed / 2;
   else
      nd1 = enmed / 2;
   end
else
   maxnd1 = x(mxm(cp + 2 * mand1));
   if color == 1
      nd1 = min(maxnd1); 
   else
      nd1 = max(maxnd1);
   end  
end
mand2 = find(en(mxm(cp + 3:2:cpd)) <= enmedlat);
if length(mand2) == 0
   if color == 1 
      nd2 = enmedlat / 2;
   else
      nd2 = -enmedlat / 2;
   end
else
   maxnd2 = x(mxm(cp + 2 * mand2 + 1));
   if color == 1
      nd2 = max(maxnd2); 
   else
      nd2 = min(maxnd2);
   end
end
vr = abs(nd1 - nd2) / 6;

errx1 = 6 * (1 - (1 - gamma) * (coeffsm ^(abs(mxm(cp - 1) - mxm(cp))))) * vl;
errx2 = 6 * (1 - (1 - gamma) * (coeffsm ^(abs(mxm(cp + 1) - mxm(cp))))) * vr; 
errb = (1 - gamma) * (vl + vr) / 2;
Epr = sqrt(((dislob ^2) / 4) * (errx1 ^2 + errx2 ^2) + ((dislob ^2) * (errb ^2)));
Eprw = sqrt(((dislob ^2) / 4) * (errx1 ^2 + errx2 ^2));
Ew = t3 + (Eprw / 2);
E = Erette - (Epr / 2);
if E < 0
   E = Erette;
end




%
% energy selection
%
bo = (abs(bp - x(mxm(cp-1)))<abs(bpw)) | (abs(bp-x(mxm(cp+1)))<abs(bpw));
if ((E >= Ew)&(~bo))
   ok = 1;
else 
      ok = 0;
end   






