%%% Hunter's "Fracture Orientation and Density" code%%%

%---- See "fracinfoREADME.doc" for instructions ----%

%close all and clear all
close all;
clear all;

%read in image
imName = 'FracSmall.png';
I = imread(imName);
figure; imshow(I);

%convert to bw
bw = im2bw(I);
bw = 1-bw;              %matlab likes white grains
figure; imshow(bw);

%% 1) use regionprops to get orientation data

%identify all the grains
cc = bwconncomp(bw, 8);

%get grain data using regionprops
graindata = regionprops(cc, 'all');
grain_orient = graindata.Orientation;
avg_orient = mean(grain_orient);
avg_strike = 90 - avg_orient  %converts to degrees strike measured cw from N


%% 2) use a moving average to compute fracture density; choose appropriate NPix

NPix=50;                             %sets size of filter
h = ones(NPix,NPix) / (NPix*NPix);
filtbw = imfilter(bw,h);
figure, imshow(filtbw); title('filtered image (averaged)');


%% 3) use alternate method of computing fracture density (lab 5 RM method)

%use original bw image. Use "cc display image" routine to isolate each
%grain and index them with numbers (1,2,3,4,5,etc).  With the indexed
%image, use the "HW5 region merging" routine to check how many fractures
%are in the viscinity of each fracture.

%create FracNear (duplicate of bw) and get size number of fractures
FracNear = bw;
x=size(cc.PixelIdxList);
NumFrac = x(2);             %retrieves number of fractures

%label each fracture with increasing index (1, 2, 3, 4, 5, etc)
for i = 1:NumFrac
    bw(cc.PixelIdxList{i}) = i;
end

figure, imshow(bw);  title('indexed image');

%thicken by a  "Nthick" parameters to see how many neighbors each has
W=bw;
Nthick = 30;  %sets the maximum distance to "nearby fractures"

%use "bwmorph 'thicken'" and a double for-loop to count nearby fractures
for i = 1:NumFrac
    
     fprintf(strcat('Checking Region ', num2str(i), '\n'));
     
         WW = W;
         WW(W~=i) = 0;                  %creates a matrix of i and zeroes, isolates single grain
         BWWW = im2bw(WW, 0);           %converts array to bw with x>0 set to 1
         thickBWWW = bwmorph(BWWW, 'thicken', Nthick);      %thicken by Nthick pixels
         
    % find overlapping regions (loop over all regions again)
    for j = 1:NumFrac       
        %% skip the rest if  i and j are equal
        if j == i
            continue
        end
      
        
        %% if the grains are adjacent, +1 to the "ith" grain's bin
        MM = W;     
        MM(W~=j) = 0;                 %create a matrix of j and zeroes

        BWMM = im2bw(MM, 0);          %converts array to bw, x>0 set to 1
        BWcom = thickBWWW + BWMM;     %combines thickened bws, if they overlap there
                                      %should be some 2's in addition to 1's and 0's
        if BWcom ~= 2                               %no 2's = not adjacent
            continue
        else
            FracNear(W==i) = FracNear(W==i) + 1;    %keeps track of how many fractures are nearby the "i"th fracture
        end
        
        
  end % end of inner loop
end % end of outer loop

figure, imshow(FracNear); title('Nearby Fracures');
impixelinfo;           %use "impixelinfo" to see how many fractures are 
                       %adjacent to each individual fracture
                       