wsh_skeleton(bw, thresh, radius)

function [skel, ridges_grad] = wsh_skeleton(bw, thresh, radius)
%effect:   compute the Euclidean skeleton binary image "bw" using watershed algorithm
%Input:
%       bw: binary image
%       thresh: threshold to extract skeleton area
%       radius: the radius of neighborhood of a pixel (1 by default)
%Output:
%       skel: the binary image of the final skeleton
%       rigid_grad: the gradient of the rigides
%Function is written by Su Dongcai on 2009/11/12
%Notes:
%the function work in 4 steps:
%1. compute the inverse eucilidean distance trasform ("inv_distance") of "bw"
%2. compute the ridges(watersheds) of "inv_distance" by watershed algorithm
%3. compute the gradient of ridges "rigid_grad"
%4. compute the final skeleton by Thresholding, dilation and thining "rigid_grad".
%The meaning of "gradient of ridges" is defined due to
%Nicholas R. Howe and Alex Telea: "A point on the skeleton sits at the center of a circle that touches the edge of the figure at multiple points. 
%The intensity in the gray-level skeleton image ("rigides_grad" in this file) is based on the shortest distance you would have to travel around the perimeter of 
%the figure to connect the most distant two points."

%Example:
%read in a binary image:
%bw = imread('h.gif');bw=~bw;
%[skel, ridges_grad] = wsh_skeleton(bw, 5, 1);

%Acknowledgement:
%[1] This work was inspired by Nicholas R. Howe. Please go to his website
%http://maven.smith.edu/~nhowe/research/code/
%for more details.
%[2] Xiaochen He, for the images provided in his file: 
%"A Corner Detector based on Global and Local Curvature Properties"
%http://www.mathworks.com/matlabcentral/fileexchange/7652-a-corner-detector
%-based-on-global-and-local-curvature-properties

%Any suggestions, questions, bug reports is welcome to my mailbox:
[e-mail]


%extract the contour "curve" of "bw"
curve = extract_curve(bw,1);
inv_distance = bwdist(~bw);
inv_distance_org = inv_distance;
lenOcurve = size(curve, 1);
inv_distance(~bw) = Inf;
%extract the whole perimeter, and label it as "-1"
bw_perim = bwperim(bw);
inv_distance(bw_perim) = -1;
%label the odd part of perimeter as different number "i"
for i=1:2:lenOcurve
    coord_r = curve(i, 1);
    coord_c = curve(i, 2);
    inv_distance(coord_r, coord_c) = -2-i/lenOcurve;
end
%start watershed and get the ridges(watershed)
L = watershed(inv_distance);
Null_field = max(L(:))+1;
L(~bw) = Null_field;
ridges = (L==0);
%comput the gradient of the ridges:
ridges_grad = cmp_skelGrad(ridges, L, Null_field, radius);
%thresholding:
skel = ridges_grad>thresh;
%post morphorlogical processing:
se = strel('square', 2);
skel = imdilate(skel, se);
skel = bwmorph(skel, 'thin', Inf);
%display:
figure, 
subplot(2, 2, 1), imshow(bw), title('the original binary image');
subplot(2, 2, 2), imshow(inv_distance_org, []), title('the eucilidean distance map');
subplot(2, 2, 3), imshow(ridges, []), title('the rigides generated by watershed');
subplot(2, 2, 4), imshow(skel), title('the final skeleton');


function ridge_grad = cmp_skelGrad(ridges, L, Null_field, radius)
%effet:     compute the gradient of ridges
%Inputs:
%       ridges: binary image
%       L: the cooresponding labeled image generated by watershed algorithm
%       Null_field: number labeled the outside area
%       radius: the radius of neighbor used to comput "ridge_grad"
%Output:
%       ridge_grad: the gradient information
%function is written by Su Dongcai on 2009/11/12

ridge_grad = zeros(size(ridges));
[r, c] = find(ridges);
for i=1:length(r)
    ridge_grad(r(i), c(i)) = s_cmp_skelGrad(ridges, L, Null_field, [r(i) ; c(i)], radius);
end

function ridge_grad_val = s_cmp_skelGrad(ridges, L, Null_field, point, radius)
%effect:    the sub_function of "cmp_skelGrad"
%function is written by Su Dongcai on 2009/11/12
lenOedge = 2*radius + 1;
vertical_right_up_idx = [ [point(1)+radius :-1: point(1)-radius]', (point(2)+radius)*ones(lenOedge, 1)];
vertical_left_down_idx = [ [point(1)-radius : point(1)+radius]', (point(2)-radius)*ones(lenOedge, 1)];
horizonal_up_left_idx = [ (point(1)-radius)*ones(lenOedge, 1), [point(2)+radius : -1 : point(2)-radius]' ];
horizonal_down_right_idx = [ (point(1)+radius)*ones(lenOedge, 1), [point(2)-radius : point(2)+radius]' ];

perim_idx = [vertical_right_up_idx; horizonal_up_left_idx; vertical_left_down_idx; horizonal_down_right_idx];
adjacent_area = [];

for i=1:size(perim_idx, 1)
    if inRange(perim_idx(i, :), size(L))
        area_idx = L(perim_idx(i, 1), perim_idx(i, 2));
        if (area_idx ~= Null_field) && (area_idx > eps)
            adjacent_area = [adjacent_area, area_idx];
        end
    end
end

ridge_grad_val = cmp_skel_grad(adjacent_area, Null_field);

function ridge_grad_val = cmp_skel_grad(adjacent_area, Null_field)
%function is written by Su Dongcai on 2009/11/12
numOarea = numel(adjacent_area);
if numOarea < 2
    ridge_grad_val = 0;
    return;
end

area_diff = zeros(1, numOarea);
adjacent_area = [adjacent_area, adjacent_area(end)];

for i=1:numOarea
    area_diff(i) = roundDistance(adjacent_area(i+1)-adjacent_area(i), Null_field);
end

[ridge_grad_val, break_idx] = max(area_diff);
break_area_idx = adjacent_area(break_idx);

for i=1:numOarea
    dist = roundDistance(adjacent_area(i)-break_area_idx, Null_field);
    if dist > ridge_grad_val
        ridge_grad_val = dist;
    end
end

function dist = roundDistance(d, len)
%function is written by Su Dongcai on 2009/11/12
d = abs(d);
dist = min(d, len-d);