% ==================================================================================================
% Module:   SK.m
%
% Usage:    Compute the morphological skeleton of an input image
%
% Purpose: 
%     Find a morphological skeleton using Maragos' serial decomposition documented in
%     Pitas and Venetsanopoulos (P&V)
%
%     
% Input Variables:
%     X      input image (binary data)
%     B      structuring element (binary data)
%     
% Returned Results:
%     SK(X)     Skeleton of input image (see P&V pp. 190-191)
%     Si(X)     Matrix of each step in skeletonization (can be used for
%               reconstruction per P&V 6.8.10)
%
% Restrictions/Notes:
%     This function is geared for binary images. 
%
% See Also:
% 
% References:
% 
% Author:       Justin Ford
% Date:         2/14/04
% Revisions:    none
% ==================================================================================================

function [unions, Si] = SK(X,B);

disp(sprintf('\n\tSK.m: initializing...\n'))

% Initialize matrix to hold Si(X)'s
Si = [];

% Compute S0(X)     P&V 6.8.7
disp(sprintf('\tSK.m: computing S0(X)...'))
X1_erode = erosion(X,B);
X1 = dilation(X1_erode,B);
X1 = X1((1+floor((size(B,1)/2))):(size(X1,1)-floor((size(B,1)/2))),(1+floor((size(B,1)/2))):(size(X1,2)-floor((size(B,1)/2))));     % Resize dilation result to 256x256
Si = [Si setdifference(X,X1)];

% Compute S1(X) - SN(X)     P&V 6.8.8-6.8.9
step = 0;

while ( max(max(Si(1:256,(step*256+1):((step+1)*256)))) ~= 0 )
    step = step + 1;
    disp(sprintf('\tSK.m: computing S%d(X)...',step))
    X2_erode = erosion(X1_erode,B);
    X2 = dilation(X2_erode,B);
    X2 = X2((1+floor((size(B,1)/2))):(size(X2,1)-floor((size(B,1)/2))),(1+floor((size(B,1)/2))):(size(X2,2)-floor((size(B,1)/2))));     % Resize dilation result to 256x256
    Si = [Si setdifference(X2,X1_erode)];
    X1_erode = X2_erode;
end

disp(sprintf('\tSK.m: ...S%d(X) empty, removing from matrix Si...',step))
Si = Si(1:256,1:((step)*256));

% Compute SK(X)=U Si(X), i=0,1,...,N    P&V 6.8.5
unions = Si(1:256,1:256);
disp(sprintf('\n\tSK.m: computing SK(X)...'))
for i = 1:(step-1)
    unions = union(unions,Si(1:256,(i*256+1):((i+1)*256)));
end

disp(sprintf('\n\tSK.m: done.'))