%% Looking at our results with hist and scatter
% We'll start with stats2, where we removed those things that were
% definitely too small.
load rice_stats

areas2 = [stats.Area];
intens2 = [stats.MeanIntensity];

figure
subplot(2,1,1)
histogram(areas2)
title('Rice areas')
subplot(2,1,2)
histogram(intens2)
title('Rice intensities')

% A second argument in the hist command will let us set the number of bins
figure
subplot(2,1,1)
histogram(areas2,20)
title('Histogram with 20 bins')
subplot(2,1,2)
histogram(areas2,length(areas2))
title('Histogram with 100% bins')

% Scatter will let us compare two quantities togther
figure
scatter(areas2,intens2)
xlabel('Rice Areas')
ylabel('Rice Intensitiy')

% There are two populations in the intensity.  Lets see if we can look at
% those below 0.62.

dim_grains = find(intens2 < 0.62)

% Comes back with a list of 20 grains.  We'll use the trick of a for loop
% adding into a new mask to generate the dim grains.
nDim = length(dim_grains);
dimMask = zeros(size(img));
for j = 1:nDim
    tmpMask = (img_label == dim_grains(j));
    dimMask = dimMask + tmpMask;
end

% Display all the dim grains
figure
imshow([dimMask imgDouble])


%% Fitting curves in MATLAB
% First load in our data on the rice grains
clear
close all
load rice_stats

areas = [stats.Area];
tooSmall = find(areas < 50);

stats2 = stats;
stats2(tooSmall) = [];

areas = [stats2.Area];
intens = [stats2.MeanIntensity];

figure
subplot(2,2,1)
hist(areas,20)
title('Rice areas')
subplot(2,2,2)
hist(intens,20)
title('Rice Intensity')
subplot(2,2,[3 4])
scatter(areas,intens)
xlabel('Areas')
ylabel('Intensity')

% If we think that the rice grain sizes are normally distributed, we can
% fit to a gaussian.  First we extract the x and y values of the histogram.

[areaYvals, areaXvals] = hist(areas,20);
[intensYvals, intensXvals] = hist(intens,20);

figure
plot(areaXvals,areaYvals)

% We can use the command fit to fit to a single guassian.
gaussFit = fit(areaXvals',areaYvals','gauss1')
fitCenter = gaussFit.b1
fitSpread = gaussFit.c1

figure
plot(areaXvals,areaYvals,'x')
hold all
plot(gaussFit)
hold off

% We can also try double gaussian
% We can use the command fit to fit to a single guassian.
gaussFit2 = fit(areaXvals',areaYvals','gauss2')
fitCenter2 = gaussFit2.b1
fitSpread2 = gaussFit2.c1

figure
plot(areaXvals,areaYvals,'x')
hold all
plot(gaussFit2)
hold off

% The residuals are how far off the fit is at each point.  It's useful to
% see if we might need another term for the fit.  We can calculate out the
% fit values using the feval command
fitData = feval(gaussFit,areaXvals);
resid = areaYvals - fitData';

figure
plot(areaXvals,resid)
hold all
line([areaXvals(1) areaXvals(end)],[0 0],...
    'Linestyle','--',...
    'Color','k')
hold off


% The curve fitting tool allows you to do a lot of these things in a
% nice graphical interface.

cftool

areasFit(areaXvals,areaYvals)

