function M = FindSkin(D, NumberOfPoint, DataLength, TotalMaxPoint)

%Composed by Hao Zhang, March 07, 2005
%input data is the B scan data that has be reshaped from the direct
%FastFrame acquisition from the scope.

%Output is the calcaulted starting point of the skin surface. In order to 
%totally remove the skin, a certain depth should be included, such as 35
%points after the starting point.

%************************************************************************
%If you use our program and publish your work, please cite H.F. Zhang, 
%K. Maslov, and L.V. Wang, ¡°Automatic algorithm for skin profile 
%detection in photoacoustic microscopy,¡± %Journal of Biomedical Optics 
%14, 024050  (2009).
%************************************************************************

Edge = zeros(NumberOfPoint, 1) ;

for t = 1:NumberOfPoint
    
    %Loop controling var for the max values we are looking for
    max_count = 1 ;  
    d = D(:,t) ;
    tmp = sort(d, 'descend') ;
    
    loc = zeros(TotalMaxPoint, 1) ;
    m_value = zeros(TotalMaxPoint, 1) ;
    
    index = 1 ;
    
    m_value(1) = tmp(1) ;
    
    while max_count<=TotalMaxPoint
        tmpv = max(d) ;
        tmploc = min(find(d == tmpv)) ;
        m_value(max_count) = tmpv ;
        loc(max_count) = tmploc ;
        l = tmploc-5 ;
        if l<1
            l=1 ;
        end ;
        r = tmploc+5 ;
        if r>DataLength
            r = DataLength ;
        end ;
        d(l:r) = 0 ;
        max_count = max_count + 1 ;
    end ;
        
    loc1 = sort(loc) ;
    
    space = zeros((TotalMaxPoint-1), 1) ;
    space_l = 1 ;
    max_space = loc1(space_l+1) - loc1(space_l) ;
    for p = 1:(TotalMaxPoint-1)
        if (loc1(p+1) - loc1(p)) > max_space
            max_space = loc1(p+1) - loc1(p) ;
            space_l = p ;
        end ;
    end ;

    %W1(1:(loc1(space_l)+25),t) = 0 ;
    Edge(t) = loc1(space_l) ;
end ;

M = Edge ;