function MineralAnalysis(picture)
% Function "MineralAnalysis" takes the name of a .jpg photo of a rock in thin section
% and allows the user to define Mineral names to differnet color ranges and
% then breaks the image into parts by mineral type and returns a bar graph
% of the percentage composition of the rock.

% Initial settings
tolerance = 0;
numberofcolorstofilter = 15; %defualt
font_size = 15;
inputfile = picture;
photo = 0;
bwmasks = {};
names = {};
colors = {};
mineralmasks = {};
mineralparts = {};
pixelmineralsize = [];
cc = {};
connectivitymins = {};
graindata = {};

% Load original image
original_image = imread(inputfile);
imshow(original_image)
title(inputfile);
set(gcf,'name','Original Image','numbertitle','off');
%saveas(gcf, ['processing_', sprintf('%03d',photo), '.png']);
photo = photo + 1;
imagedimensions = size(original_image);

% Ask how many minerals can be identified
num_minerals = str2num(input('How many different minerals are there in the rock sample? : ','s'));

% Ask for the names of the minerals
for current_mineral=1:num_minerals
    prompt = sprintf('Name of %i Mineral Species : ',current_mineral);
    names{current_mineral} = input(prompt,'s');
end

% Estimate colors seen to create a filterd image that pushes all colors to specific RGB values
% + 2 because people forget black and white!
numberofcolorstofilter = str2num(input('Please estimate how many different distinct mineral colors there\nare in the rock sample in order to filter the image : ','s')) + 2;

% Create the filtered image
fprintf('Creating filtered image for easier processing \n');
[indbase mapbase] = rgb2ind(original_image, numberofcolorstofilter, 'nodither');
filtered_image = ind2rgb(indbase, mapbase);
figure
imshow(filtered_image)
title('Filtered Image');
set(gcf,'name','Filtered Image','numbertitle','off');
%saveas(gcf, ['processing_', sprintf('%03d',photo), '.png']);
photo = photo + 1;

% Select the characteristic mineral colors from the photo for the number of
% minerals that are observed
for current_mineral = 1:num_minerals
    fprintf('Please select the characteristic colors for %s and then press "Enter"\n', names{current_mineral});
    colors{current_mineral} =  impixel;
end

% Initiliaze some properties and images we will use later
for current_mineral = 1:num_minerals
    pixelmineralsize(current_mineral) = 0;
    mineralmasks{current_mineral} = zeros(imagedimensions(1), imagedimensions(2), imagedimensions(3), 'uint8');
    mineralparts{current_mineral} = zeros(imagedimensions(1), imagedimensions(2), imagedimensions(3), 'uint8');
end

fprintf('Creating masks for images \n');
% Iterate through the filtered image to find pixels of a color selected for
% each given mineral. record the locations of these colors by marking them
% white in our mineral mask 
for current_mineral=1:num_minerals
    matrixofcolors = size(colors{1,current_mineral});
    for imagex = 1:imagedimensions(1) 
        for imagey = 1:imagedimensions(2)
            for testpoint = 1:matrixofcolors(1)
                if checkrgb(filtered_image, colors{1,current_mineral}, testpoint, imagex, imagey , tolerance) > 0 
                    mineralmasks{current_mineral}(imagex, imagey, 1) = 255;
                    mineralmasks{current_mineral}(imagex, imagey, 2) = 255;
                    mineralmasks{current_mineral}(imagex, imagey, 3) = 255;
                end
            end
        end
    end
end

% Make bw image masks from the mineralmasks array
figure
for current_mineral = 1:num_minerals
    n = 3;
    h = ones(n,n) / (n*n);
    level = graythresh(mineralmasks{current_mineral});
    level = 4/256;
    bwmasks{current_mineral} = im2bw(mineralmasks{current_mineral}, level);
    bwmasks{current_mineral} = bwareaopen(bwmasks{current_mineral}, 300);
    bwmasks{current_mineral} = imfilter(bwmasks{current_mineral},h);
    subplot(num_minerals, 1, current_mineral);
    imshow(bwmasks{current_mineral});
    title(['Mineral Mask of ',names{current_mineral}],'fontsize',12);
    set(gcf,'name','Mineral Masks','numbertitle','off');
end
%saveas(gcf, ['processing_', sprintf('%03d',photo), '.png']);
photo = photo + 1;


% Make regionprops to retrieve data (not used in this function could
%be utilized if I continued to improve this script)
for current_mineral = 1:num_minerals
    cc{current_mineral} = bwconncomp(bwmasks{current_mineral});
    graindata{current_mineral} = regionprops(cc{current_mineral}, 'all');
    connectivitymins{current_mineral} = cc{current_mineral}.Connectivity;
end

% Make mineral parts images
for current_mineral = 1:num_minerals
    matrixofcolors = size(colors{1,current_mineral});
    fprintf('Creating isolated images of crystals of %s using the image mask... \n', names{current_mineral});
    for imagex = 1:imagedimensions(1)
        for imagey = 1:imagedimensions(2)
            if bwmasks{current_mineral}(imagex,imagey) == 1
                mineralparts{current_mineral}(imagex,imagey,1) = original_image(imagex, imagey, 1);
                mineralparts{current_mineral}(imagex,imagey,2) = original_image(imagex, imagey, 2);
                mineralparts{current_mineral}(imagex,imagey,3) = original_image(imagex, imagey, 3);
                pixelmineralsize(current_mineral) = pixelmineralsize(current_mineral) + 1;
            end
        end
    end
end

% Plot mineral part images
figure
for current_mineral = 1:num_minerals
    subplot(num_minerals,1,current_mineral);
    imshow(mineralparts{current_mineral});
    title(['Minerals of ',names{current_mineral}, ' within the rock'],'fontsize',12);
    set(gcf,'name','Mineral Sub-Parts','numbertitle','off');
end
%saveas(gcf, ['processing_', sprintf('%03d',photo), '.png']);
photo = photo + 1;

% Make a colored bar graph showing % area covered of total image by each
% selected mineral
figure
for current_mineral = 1:num_minerals
    percentageofrock(current_mineral) = (pixelmineralsize(current_mineral) / (imagex * imagey)) * 100;
    colorselectionsize = size(colors{1,current_mineral});
    if colorselectionsize(1) > 1
        barColorMap(current_mineral,:) = mean(colors{1,current_mineral});
    else
        barColorMap(current_mineral,:) = colors{1,current_mineral};
    end
    
    barnumber(current_mineral) = bar(current_mineral, percentageofrock(current_mineral), 'BarWidth', 0.5);
    set(barnumber(current_mineral), 'FaceColor', barColorMap(current_mineral,:));
    percent = num2str(ceil(percentageofrock(current_mineral)));
    barTopper = sprintf('%s',percent,' %');
    text(current_mineral, percentageofrock(current_mineral) + 3, barTopper, 'FontSize', 15);
    hold on
end

ylim([0 100]);
title('Percentage of minerals, determined by color, contained in the thin section', 'FontSize', font_size)
xlabel('Mineral Type', 'FontSize', font_size);
ylabel('Portion of Rock (%)', 'FontSize', font_size)
set(gca, 'XTick', 1:num_minerals, 'XtickLabel', names);
set(gcf,'name','Histogram of Mineral Percentages','numbertitle','off');
%saveas(gcf, ['processing_', sprintf('%03d',photo), '.png']);
photo = photo + 1;

    function truefalse = checkrgb(original_image, pointselections, testpoint, imagex, imagey,tolerance)
        % Checks to see if the rgb of a given pixel in the original_image matches
        % one of the selections the user made. If the colors match returns a +1,
        % otherwise returns a -1
        if (abs(original_image(imagex,imagey,1) - pointselections(testpoint,1)) <= tolerance) && (abs(original_image(imagex,imagey,2) - pointselections(testpoint,2)) <= tolerance) && (abs(original_image(imagex,imagey,3) - pointselections(testpoint,3)) <= tolerance)
            truefalse = +1;
        else
            truefalse = -1;
        end
    end
end

