BrainStrokeSgementation.rar

function [out] = GLCM_Features(glcmin,pairs) warning off; % If 'pairs' not entered: set pairs to 0 if ((nargin > 2) || (nargin == 0)) error('Too many or too few input arguments. Enter GLCM and pairs.'); elseif ( (nargin == 2) ) if ((size(glcmin,1) <= 1) || (size(glcmin,2) <= 1)) error('The GLCM should be a 2-D or 3-D matrix.'); elseif ( size(glcmin,1) ~= size(glcmin,2) ) error('Each GLCM should be square with NumLevels rows and NumLevels cols'); end elseif (nargin == 1) % only GLCM is entered pairs = 0; % default is numbers and input 1 for percentage if ((size(glcmin,1) <= 1) || (size(glcmin,2) <= 1)) error('The GLCM should be a 2-D or 3-D matrix.'); elseif ( size(glcmin,1) ~= size(glcmin,2) ) error('Each GLCM should be square with NumLevels rows and NumLevels cols'); end end format long e if (pairs == 1) newn = 1; for nglcm = 1:2:size(glcmin,3) glcm(:,:,newn) = glcmin(:,:,nglcm) + glcmin(:,:,nglcm+1); newn = newn + 1; end elseif (pairs == 0) glcm = glcmin; end size_glcm_1 = size(glcm,1); size_glcm_2 = size(glcm,2); size_glcm_3 = size(glcm,3); % checked out.autoc = zeros(1,size_glcm_3); % Autocorrelation: [2] out.contr = zeros(1,size_glcm_3); % Contrast: matlab/[1,2] out.corrm = zeros(1,size_glcm_3); % Correlation: matlab out.corrp = zeros(1,size_glcm_3); % Correlation: [1,2] out.cprom = zeros(1,size_glcm_3); % Cluster Prominence: [2] out.cshad = zeros(1,size_glcm_3); % Cluster Shade: [2] out.dissi = zeros(1,size_glcm_3); % Dissimilarity: [2] out.energ = zeros(1,size_glcm_3); % Energy: matlab / [1,2] out.entro = zeros(1,size_glcm_3); % Entropy: [2] out.homom = zeros(1,size_glcm_3); % Homogeneity: matlab out.homop = zeros(1,size_glcm_3); % Homogeneity: [2] out.maxpr = zeros(1,size_glcm_3); % Maximum probability: [2] out.sosvh = zeros(1,size_glcm_3); % Sum of sqaures: Variance [1] out.savgh = zeros(1,size_glcm_3); % Sum average [1] out.svarh = zeros(1,size_glcm_3); % Sum variance [1] out.senth = zeros(1,size_glcm_3); % Sum entropy [1] out.dvarh = zeros(1,size_glcm_3); % Difference variance [4] %out.dvarh2 = zeros(1,size_glcm_3); % Difference variance [1] out.denth = zeros(1,size_glcm_3); % Difference entropy [1] out.inf1h = zeros(1,size_glcm_3); % Information measure of correlation1 [1] out.inf2h = zeros(1,size_glcm_3); % Informaiton measure of correlation2 [1] %out.mxcch = zeros(1,size_glcm_3);% maximal correlation coefficient [1] %out.invdc = zeros(1,size_glcm_3);% Inverse difference (INV) is homom [3] out.indnc = zeros(1,size_glcm_3); % Inverse difference normalized (INN) [3] out.idmnc = zeros(1,size_glcm_3); % Inverse difference moment normalized [3] % correlation with alternate definition of u and s %out.corrm2 = zeros(1,size_glcm_3); % Correlation: matlab %out.corrp2 = zeros(1,size_glcm_3); % Correlation: [1,2] glcm_sum = zeros(size_glcm_3,1); glcm_mean = zeros(size_glcm_3,1); glcm_var = zeros(size_glcm_3,1); % http://www.fp.ucalgary.ca/mhallbey/glcm_mean.htm confuses the range of % i and j used in calculating the means and standard deviations. % As of now I am not sure if the range of i and j should be [1:Ng] or % [0:Ng-1]. I am working on obtaining the values of mean and std that get % the values of correlation that are provided by matlab. u_x = zeros(size_glcm_3,1); u_y = zeros(size_glcm_3,1); s_x = zeros(size_glcm_3,1); s_y = zeros(size_glcm_3,1); % % alternate values of u and s % u_x2 = zeros(size_glcm_3,1); % u_y2 = zeros(size_glcm_3,1); % s_x2 = zeros(size_glcm_3,1); % s_y2 = zeros(size_glcm_3,1); % checked p_x p_y p_xplusy p_xminusy p_x = zeros(size_glcm_1,size_glcm_3); % Ng x #glcms[1] p_y = zeros(size_glcm_2,size_glcm_3); % Ng x #glcms[1] p_xplusy = zeros((size_glcm_1*2 - 1),size_glcm_3); %[1] p_xminusy = zeros((size_glcm_1),size_glcm_3); %[1] % checked hxy hxy1 hxy2 hx hy hxy = zeros(size_glcm_3,1); hxy1 = zeros(size_glcm_3,1); hx = zeros(size_glcm_3,1); hy = zeros(size_glcm_3,1); hxy2 = zeros(size_glcm_3,1); %Q = zeros(size(glcm)); for k = 1:size_glcm_3 % number glcms glcm_sum(k) = sum(sum(glcm(:,:,k))); glcm(:,:,k) = glcm(:,:,k)./glcm_sum(k); % Normalize each glcm glcm_mean(k) = mean2(glcm(:,:,k)); % compute mean after norm glcm_var(k) = (std2(glcm(:,:,k)))^2; for i = 1:size_glcm_1 for j = 1:size_glcm_2 out.contr(k) = out.contr(k) + (abs(i - j))^2.*glcm(i,j,k); out.dissi(k) = out.dissi(k) + (abs(i - j)*glcm(i,j,k)); out.energ(k) = out.energ(k) + (glcm(i,j,k).^2); out.entro(k) = out.entro(k) - (glcm(i,j,k)*log(glcm(i,j,k) + eps)); out.homom(k) = out.homom(k) + (glcm(i,j,k)/( 1 + abs(i-j) )); out.homop(k) = out.homop(k) + (glcm(i,j,k)/( 1 + (i - j)^2)); % [1] explains sum of squares variance with a mean value; % the exact definition for mean has not been provided in % the reference: I use the mean of the entire normalized glcm out.sosvh(k) = out.sosvh(k) + glcm(i,j,k)*((i - glcm_mean(k))^2); %out.invdc(k) = out.homom(k); out.indnc(k) = out.indnc(k) + (glcm(i,j,k)/( 1 + (abs(i-j)/size_glcm_1) )); out.idmnc(k) = out.idmnc(k) + (glcm(i,j,k)/( 1 + ((i - j)/size_glcm_1)^2)); u_x(k) = u_x(k) + (i)*glcm(i,j,k); % changed 10/26/08 u_y(k) = u_y(k) + (j)*glcm(i,j,k); % changed 10/26/08 % code requires that Nx = Ny % the values of the grey levels range from 1 to (Ng) end end out.maxpr(k) = max(max(glcm(:,:,k))); end % glcms have been normalized: % The contrast has been computed for each glcm in the 3D matrix % (tested) gives similar results to the matlab function for k = 1:size_glcm_3 for i = 1:size_glcm_1 for j = 1:size_glcm_2 p_x(i,k) = p_x(i,k) + glcm(i,j,k); p_y(i,k) = p_y(i,k) + glcm(j,i,k); % taking i for j and j for i if (ismember((i + j),[2:2*size_glcm_1])) p_xplusy((i+j)-1,k) = p_xplusy((i+j)-1,k) + glcm(i,j,k); end if (ismember(abs(i-j),[0:(size_glcm_1-1)])) p_xminusy((abs(i-j))+1,k) = p_xminusy((abs(i-j))+1,k) +... glcm(i,j,k); end end end % % consider u_x and u_y and s_x and s_y as means and standard deviations % % of p_x and p_y % u_x2(k) = mean(p_x(:,k)); % u_y2(k) = mean(p_y(:,k)); % s_x2(k) = std(p_x(:,k)); % s_y2(k) = std(p_y(:,k)); end % marginal probabilities are now available [1] % p_xminusy has +1 in index for matlab (no 0 index) % computing sum average, sum variance and sum entropy: for k = 1:(size_glcm_3) for i = 1:(2*(size_glcm_1)-1) out.savgh(k) = out.savgh(k) + (i+1)*p_xplusy(i,k); % the summation for savgh is for i from 2 to 2*Ng hence (i+1) out.senth(k) = out.senth(k) - (p_xplusy(i,k)*log(p_xplusy(i,k) + eps)); end end % compute sum variance with the help of sum entropy for k = 1:(size_glcm_3) for i = 1:(2*(size_glcm_1)-1) out.svarh(k) = out.svarh(k) + (((i+1) - out.senth(k))^2)*p_xplusy(i,k); % the summation for savgh is for i from 2 to 2*Ng hence (i+1) end end % compute difference variance, difference entropy, for k = 1:size_glcm_3 % out.dvarh2(k) = var(p_xminusy(:,k)); % but using the formula in % http://murphylab.web.cmu.edu/publications/boland/boland_node26.html % we have for dvarh for i = 0:(size_glcm_1-1) out.denth(k) = out.denth(k) - (p_xminusy(i+1,k)*log(p_xminusy(i+1,k) + eps)); out.dvarh(k) = out.dvarh(k) + (i^2)*p_xminusy(i+1,k); end end % compute information measure of correlation(1,2) [1] for