A Hybrid Skin Lesion Segmentation Tool: Using PCA and Iterative Canny Edge

This is an improved version of of a previous skin lesion segmentation algorithm that I developed.

Below is the script for the Skin Lesion Segmentation Algorithm.
If you use this script in anyway please cite the author.

This script was developed and copyrighted by Tyler L. Coye (2015)

This is a third release of this tool. In this release, I improved the segmentation by applying an iterative canny edge
to the image mask. This has the added benefit of improving contour matching between the mask and original lesion.

The method outlined below is a novel approach to lesion segmentation. No other script, to my knowledge, utilizes principle
component analysis for color to gray conversion or an iterative canny edge. The threshold level calculation is also unique to this script. It is typical to divide the sum of threshold levels by '4.' For this purpose, I have found it better to divide by '2.'

Feel free to leave a comment and let me know how it has worked for you or how you have applied it to your work/research.

-Tyler

%% Hybrid Lesion Detection 2.0 is Copyrighted by Tyler Coye, 2015.
% If you use this script please notify and cite the author.
% If you have any questions regarding this script you can contact me at
% [e-mail].
% This Script Uses:
% -Iterative Median Filtering
% -2-D Wavelet Transformation
% -2-D inverse Wavelet transformation
% -Otsu Thresholding on individual dwt2 levels
% -Canny edge detection
% This is an improved version of the Hybrid Skin Lesion Detection
% Algorithm.
% The following changes were made:
% Added Color to gray via PCA (novel method applied to this type of problem)
% Added morphological closing
% Removed ROI cropping
% Iterative Canny Edge (novel method applied to this type of problem)

% Read image
im = im2double(imread('th3.jpg'));

% Convert RGB to Gray via PCA
lab = rgb2lab(im);
f = 0;
wlab = reshape(bsxfun(@times,cat(3,1-f,f/2,f/2),lab),[],3);
[C,S] = pca(wlab);
S = reshape(S,size(lab));
S = S(:,:,1);
gray = (S-min(S(:)))./(max(S(:))-min(S(:)));

% Morphological Closing
se = strel('disk',1);
close = imclose(gray,se);

% Complement Image
K= imcomplement(close)

%% 2-D wavelet Decomposition using B-Spline
[cA,cH,cV,cD] = dwt2(K,'bior1.1');

%% Otsu thresholding on each of the 4 wavelet outputs
thresh1 = multithresh(cA);
thresh2 = multithresh(cH);
thresh3 = multithresh(cV);
thresh4 = multithresh(cD);

%% Calculating new threshold from sum of the 4 otsu thresholds and dividing by 2
level = (thresh1 + thresh2 + thresh3 + thresh4)/2;

% single level inverse discrete 2-D wavelet transform
X = idwt2(cA,cH,cV,cD,'bior1.1')

% Black and White segmentation
BW=imquantize(X,level);

%% Iterative Canny Edge (Novel Method)
BW1 = edge(edge(BW,'canny'), 'canny');

%% Post-Processing
BW3 = imclearborder(BW1);
CC = bwconncomp(BW3);
S = regionprops(CC, 'Area');
L = labelmatrix(CC);
BW4 = ismember(L, find([S.Area] >= 100));
BW5 = imfill(BW4,'holes');

%% Present Final Image
[B,L,N] = bwboundaries(BW5);
figure; imshow(im); hold on;
for k=1:length(B),
boundary = B{k};
plot(boundary(:,2),...
boundary(:,1),'g','LineWidth',2);
end