Function for Counting Colors in a Skin Lesion

This function calculates the Color (C) score for ABCD calculation of skin lesions.

Below is the script for the function Cscore which I developed. Cscore segments and extracts a skin lesion, then outputs the final color score based on the ABCD criteria.
Copyrighted by Tyler L. Coye (2015)
If you use this script in your research in anyway, please cite me as the author.

YOU WILL NEED THE imoverlay.m that can be found HERE: https://www.mathworks.com/matlabcentral/fileexchange/50839-a-novel-retinal-blood-vessel-segmentation-algorithm-for-fundus-images/content/sample/imoverlay.m

function [Finalscore, A] = Cscore(image)

***PART I-SEGMENTATION***

% Read image
I = imread(image);
I = imresize(I, [200 200])
im = im2double(I);

% 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));
BW51 = imfill(BW4,'holes');
BW5 = imcomplement(BW51)

% overlay with a green outer mask
out = imoverlay(im, BW5, [0 1 0])

%% *** PART 2 - Lesion Color Scoring Method***

im=im2double(out)% Convert to floating points
r=im(:,:,1);
g=im(:,:,2);
b=im(:,:,3);

% Calculate number of pixels with a given color
Rblack = sum(sum(im(:,:,1)<.2 & im(:,:,2)<.2 & im(:,:,3)<.2));
Rwhite = sum(sum(im(:,:,1)>.8 & im(:,:,2)>.8 & im(:,:,3)>.8));
Rred = sum(sum(im(:,:,1) >.8 & im(:,:,2)<.2 & im(:,:,3)<.2));
Rlightbrown = sum(sum(im(:,:,1)>.6 & im(:,:,1)<1 & im(:,:,2)>0.32 & im(:,:,2)< 0.72 & im(:,:,3)>.05 & im(:,:,3)<.45));
Rdarkbrown = sum(sum(im(:,:,1)>.2 & im(:,:,1)<.6 & im(:,:,2)>0.06 & im(:,:,2)< 0.46 & im(:,:,3)>0 & im(:,:,3)<.33));
Rbluegray =sum(sum(im(:,:,1)<.2 & im(:,:,2)>0.32 & im(:,:,2)< 0.72 & im(:,:,3)>.34 & im(:,:,3)<.74));

% Calculate total pixels
[rows columns numberOfColorChannels] = size(im);
numberOfPixels = rows*columns;

% Calculate number of green pixels (this is everything outside of the
% lesion)
Rgreen = sum(sum(im(:,:,1)==0 & im(:,:,2)==1 & im(:,:,3)==0));

% Take the difference between total and green to get total lesion pixels
tlp = numberOfPixels - Rgreen

% Score individual colors
if (Rblack/tlp)*100>5
a=1
else
a=0
end

if (Rwhite/tlp)*100>5
b=1
else
b=0
end

if (Rred/tlp)*100>5
c=1
else
c=0
end

if (Rlightbrown/tlp)*100>5
d=1
else
d=0
end

if (Rdarkbrown/tlp)*100>5
e=1
else
e=0
end

if (Rbluegray/tlp)*100>5
f=1
else
f=0
end

% Build Data table
Color = {'Black';'White';'Red';'Light Brown';'Dark Brown';'Blue Gray'};
Color_Count = [Rblack;Rwhite;Rred;Rlightbrown;Rdarkbrown;Rbluegray];
Percentage = [(Rblack/tlp)*100;(Rwhite/tlp)*100;(Rred/tlp)*100;(Rlightbrown/tlp)*100;(Rdarkbrown/tlp)*100;(Rbluegray/tlp)*100];
Score = [a;b;c;d;e;f]

A = table(Color_Count,Percentage,Score,'RowNames',Color)

% Determine final score
G = [a, b, c, d, e, f];
Finalscore = sum(G>0)

%% Present Image
[B,L,N] = bwboundaries(BW51);
figure; imshow(im); hold on;
for k=1:length(B),
boundary = B{k};
plot(boundary(:,2),...
boundary(:,1),'g','LineWidth',2);
hold on
himage = imshow(im);
set(himage, 'AlphaData', 0.5);
text(0,-10,strcat('\color{red}\fontsize{12}Estimated Color Score:',num2str(Finalscore)))
end
end