%imageProcess4
% This is the code I do most of my project with.  If you are going to run
% it, you need my hard drive, or a file with all 129 LED images....
% Then you would need to change the directory line.
% I would suggest commenting out all the figure, imshows that are now
% saved, or commenting out the saves stuff as well, to go faster.

clf
clear
close all

cd /Users/cmhamiltonian/Documents/Matlab

%control is an image I made which I use in my program to select the
%area....in which I will analyze the image data.
% the selection of this area is a key factor to the processing of ....these
% images because for reasons we don't understand ,they don't perfectly line
% up.  We also don't want to analyze data outside of the lense.
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-4

%imagedimensions(standard for all these CFE images)

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

for i=105:4:N
    LEDcleanI = bmpImages(i).name;
    cleanImage = imread(LEDcleanI);  %origonal LED image
        figure, imshow(image);
    
    
    cleanname = strcat('Inew',int2str(i));%to be used when saving
    cleanbw = strcat('Ibw',int2str(i));%
    IcleanGray = rgb2gray(cleanImage);
        figure, imshow(IcleanGray);
    IcSelect = imsubtract(controlG,IcleanGray);
    figure, imshow(IcSelect);
    saveas(gcf,[cleanname,'.png'])
    
     bwClean = im2bw(IcSelect,level);
     figure, imshow(bwClean);
     saveas(gcf,[cleanbw,'.png'])
     
     POCclean = bwlabeln(bwClean,4);
     i
    
    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
    
    figure, imshow(bwClean);
    
    %# find connected components (4-connected neighborhood)
    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
    X = k
        
        Iname = strcat('Inew',int2str(k));
        InameBW = strcat('Ibw',int2str(k));
        Ifinal = strcat('Ifinal',int2str(k));
        
        LEDsampleI = bmpImages(k).name;
        sample = imread(LEDsampleI);
        
        IsampleGray = rgb2gray(sample);
        
               figure, imshow(IsampleGray);
              saveas(gcf,[Iname,'First.png'])
        
        %takes previously present particles out and sets analyzing radius
        IsSelect = imsubtract(IcSelect, IsampleGray);
        figure, imshow(IsSelect);
        saveas(gcf,[Iname,'.png']);
        q = imtool(IsSelect)
        

        bwSample = im2bw(IsSelect,level);
        figure, imshow(bwSample);
        saveas(gcf,[InameBW,'.png']);
        figure, imshow(bwSelect)
% wipes out everything in the imaged outside of the lense radius so that attributes 
%and imperfections of the lense are not counted as particles
%this radius had to be figured out by hand over many many painstaking hours
%using imageJ  and imtool
%had to origionally experiment by writing a different program and
%experimenting with it... cally 'makeControl.m' I will include this so you
%can see it too
        
        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
        
        figure, imshow(bwSample);
        saveas(gcf,[Ifinal,'.png']);
        
                POCsample = bwlabeln(bwSample,4);
        
        
%         # find connected components (4-connected neighborhood)
        CC = bwconncomp(bwSample, 4);
        
%         # find blocks with two or more connected 1's
        idx = ( cellfun(@numel,CC.PixelIdxList) > 1 );
        POCsample = sum(idx) - POCclean;
        f(k)=POCsample;
        
    end

    g = plot(f,'pb');
    hold on
    
end


