clear all;clc;


% enter colorspace option
output='cielab';

% treshold of im2bw
thresh=.65;

% Define well numbers
well=14;


% Define respective depths of samples for each well
well_14_depth=[1 2 3 4 5 6 7 8 9 10 16 17 21 22 26 27 31 40 42 45 50 60 61 70 71 85]';


depths={well_14_depth};

for i=1:length(well);
% necessary to clear variables as they may change in size with each well
clearvars -except well well_12_depth well_13_depth well_14_depth well_15_depth well_16_depth depths thresh output i n

% Name the directory(ies) where your original images are stored
directory=strcat('Images/well_',num2str(well(i)));

% get directory listing
listing = dir(strcat(directory,'/*.JPG'));

% get number of files matching the query
nFiles = numel(listing);

for n=1:nFiles
    j= fprintf([listing(n).name '\n']);
   
    % get image name
    imName = strcat(directory,'/',listing(n).name);
 
    % read image
    rgbImg = imread(imName);
    
    % Seperate out region of interest
    disp('Seperating out region of interest...');
    bw=1-im2bw(rgbImg,thresh);
    cc=bwconncomp(bw, 8);
    graindata=regionprops(cc,'Area');
    grain_areas=[graindata.Area];
    
    [max_area, idx]=max(grain_areas);

    region=false(size(bw));
    region(cc.PixelIdxList{idx})=true;
    fill=imfill(region,'holes');
    
    % use the region mask to extract the color information from the region
    disp('Extracting color information');
    red= immultiply(fill, rgbImg(:,:,1)); 
        avg_red(n)=mean(nonzeros(red));
    green = immultiply(fill, rgbImg(:,:,2));
        avg_green(n)=mean(nonzeros(green));
    blue = immultiply(fill, rgbImg(:,:,3)); 
        avg_blue(n)=mean(nonzeros(blue));

    % create the final image of the region of interest
    roiImg = cat(3, red, green, blue); 
    
    % convert region of interest image into colorspace specified above
    c='Converting image to '; b=[c  output]; disp(b);
    
    D=colorspace(strcat(output,'<-rgb'),roiImg);
    avg_L(n)=mean(nonzeros(D(:,:,1)));
    avg_a(n)=mean(nonzeros(D(:,:,2)));
    avg_b(n)=mean(nonzeros(D(:,:,3)));

end

% Concatonate well data so that it is in one matrix
well_dat=cat(2,depths{i}, avg_red', avg_green', avg_blue', avg_L', avg_a', avg_b');

saving=strcat('well_',num2str(well(i)),'_data');
save(saving,'well_dat')

% Create a scatter image of the avg color as a function of depth
disp('Creating Scatter Image');
scatter_matrix=cat(2,well_dat(:,1),well_dat(:,2),well_dat(:,3),well_dat(:,4));
a=scatter_matrix;
cforscat=a(:,2:4);

disp('Creating Color Intensity Plots');

H=figure;
subplot(1,4,1);hold on;

for k=1:.1:9
scatter(ones(size(a(:,1))).*k,a(:,1),100,cforscat./255,'Filled','Marker','s');
end

set(gca,'XColor','white');
ylabel('Depth (ft)');title(['Well ',num2str(well(i)),' Color Sequence']);
axis([0 10 0 roundn(max(well_dat(:,1)),1)]);
set(gca,'ydir','reverse','XTick',[])
hold off;

% plot the water content
subplot(1,4,2);hold on;
plot(well_dat(:,5),well_dat(:,1),'s','MarkerFaceColor','w','MarkerEdgeColor','k');
set(gca, 'ydir','reverse','XAxisLocation','top');
xlabel('Luminosity');
set(gcf,'Color','White');
% lsline
hold off;

% plot the green/red colorspace with depth
subplot(1,4,3);hold on;
plot(well_dat(:,6),well_dat(:,1),'rs','MarkerFaceColor','r','MarkerEdgeColor','k');
set(gca, 'ydir','reverse','XAxisLocation','top');
xlabel('Green<-------->Red');
% lsline
hold off;


subplot(1,4,4);hold on;%figure;
plot(well_dat(:,7),well_dat(:,1),'ys','MarkerFaceColor','y','MarkerEdgeColor','k');
set(gcf,'Color','White');set(gca, 'ydir','reverse','XAxisLocation','top');
xlabel('Blue<------->Yellow');
hold off;

% Save matrix and image
filename=strcat('well_',num2str(well(i)),'_graph');
saveas(H,filename,'fig');

end
