% imageprocess5
%This is the code I make my graphs with

clf
clear
close all

cd /Users/cmhamiltonian/Documents/Matlab

Control = imread('Control.bmp');
%  figure, imshow(control);
controlG = rgb2gray(Control);


cd /Volumes/CFE/001-CFE001-82/20120531/001-CFE001-82-LED
bmpImages = dir('*.bmp')


n = numel(bmpImages)
cleanerCounter = 0;
sampleCounter = 0;
N=n

%imagedimensions(standard for all these CFE images)

height =1900;
width  = 1800;
H = 980;
K = 1000;
Radius = 942;
level = .8
X = 0;


% CFE
maxX = 20;
maxY = 100;
width = 10;
depth = 2;
A = zeros(maxX, maxY);
A(width,depth) = 1;

for i=1:4:N-1
    graphname = strcat('POCgraph',int2str(i));%
    LEDcleanI = bmpImages(i).name;
    cleanImage = imread(LEDcleanI);  %origonal LED image    
    IcleanGray = rgb2gray(cleanImage);
    IcSelect = imsubtract(controlG,IcleanGray);   
    bwClean = im2bw(IcSelect,level);
    POCclean = bwlabeln(bwClean,4);
    for y = 1:height
        for x = 1:width
            r = sqrt(((y-K)^2)+((x-H)^2));
            if (r > Radius)
                bwClean(y,x) = 1;
            end
        end
    end
    cc = bwconncomp(bwClean, 4);
    %# find blocks with two or more connected 1's
    idx = ( cellfun(@numel,cc.PixelIdxList) > 1 );
    POCclean = sum(idx);
    
    cleanerCounter = cleanerCounter +1;%sanity check
    for k=i+1:i+3
       
        LEDsampleI = bmpImages(k).name;
        sample = imread(LEDsampleI);
        
        IsampleGray = rgb2gray(sample);
       %takes previously present particles out and sets analyzing radius
        IsSelect = imsubtract(IcSelect, IsampleGray);      
        bwSample = im2bw(IsSelect,level);
        for y = 1:height
            for x = 1:width
                r = sqrt(((y-K)^2)+((x-H)^2));
                if (r > Radius)
                    bwSample(y,x) = 1;
                end
                %
            end
        end
        POCsample = bwlabeln(bwSample,4);
        CC = bwconncomp(bwSample, 4);
 
        idx = ( cellfun(@numel,CC.PixelIdxList) > 1 );
        POCsample = sum(idx) - POCclean;  %num of grains
        %         f(k)=POCsample;
        bw = 1 - bwSample;
        %
        labeled = labelmatrix(CC);  %give all grains rand RGB color
        RGB_label = label2rgb(labeled, @spring, 'c', 'shuffle');
        %         figure, imshow(RGB_label)
        
        numPixels = cellfun(@numel,CC.PixelIdxList);
        [biggest,idx] = max(numPixels);
        bwSample(cc.PixelIdxList{idx}) = 0;
        %         figure, imshow(bw);
        
        CC = bwconncomp(bw,4);
        graindata = regionprops(CC, 'basic');
        POC = size(graindata);
        graindata(100).Area;
        grain_areas = [graindata.Area];
        nbins = 100;
        %         figure, hist(grain_areas, nbins)
        %         title('Histogram of Grain Area');
        
        totArea =sum(grain_areas);
        eul = bweuler(bw,8);
        
        f(k) = totArea;
        
        %
        %        %---------------code for CFE
        
        
    end
    
    if((i+20<129)&&(rem(i,5)==0))
        
        depth = depth+5;
        A(width,depth) = 1;
        subplot(1,2,1)
        imagesc(A');
        colormap(jet);
        axis equal;
        axis tight;
        title('CFE Position')
        ylabel('Depth as a factor of 15 ft')
        %         xlabel('CFE maintains same azimuthal coord')
        drawnow
        
        subplot(1,2,2)
        plot(f,'.-k');
        hold on
        title('POC Area')   %'Numerical Particle Count'
        xlabel('Image #')
        ylabel('Area of POC')
        
         saveas(gcf,[graphname,'.png'])
    end
end


