--- a
+++ b/Image features calculation code/Working/haralick/cooccur.m
@@ -0,0 +1,60 @@
+function h = cooccur( m,mask,dir,dist )
+%COOCCUR Calculate a co-occurrence matrix
+%
+%   Mike Lam
+%   27 June 2006
+%
+%   Calculates the co-occurrence matrix of m taken dist units along 
+%   the direction given by dir.
+    
+    % convert data array to 16-bit signed integers
+    m = int16(m);
+    
+    % convert direction to offset vector [dr,dc]
+    % take the negative of the row element because
+    % rows are the inverse of the y-axis
+    radians = dir/180*pi;
+    dr = int16(-round(sin(radians)));
+    dc = int16(round(cos(radians)));
+    
+    % find size of matrix and range of values
+    mcols = int16(size(m,2));
+    mrows = int16(size(m,1));
+        
+    % find size of matrix and range of values (again)
+    mmin = min(m(:));
+    mmax = max(m(:));
+    
+    % initialize co-occurrence matrix
+    crows = int16(mmax-mmin)+1;
+    ccols = int16(mmax-mmin)+1;
+    com = int16(zeros(crows,ccols));
+    
+    % build co-occurrence matrix
+    for r = 1:int16(mrows)
+        for c = 1:int16(mcols)  
+            %PFS 7/15/2013
+            %guard on nodule extraction, only allow pixels in the
+            %segmentation to be used
+            if( mask(r,c) )
+                for d = 1:dist
+                    % given direction
+                    nr = r+(dr*int16(d));
+                    nc = c+(dc*int16(d));
+                    if (nr > 0 && nr <= mrows && nc > 0 && nc <= mcols)
+                        i = m(r,c)-int16(mmin)+1;
+                        j = m(nr,nc)-int16(mmin)+1;
+                        com(i,j) = com(i,j) + 1;
+                    end
+                end
+            end
+        end
+    end    
+    
+    %convert to floating-point data for analysis
+    com = double(com);
+    
+    % return feature vector
+    h = calcHaralick(com); 
+    
+end