clc; clear all; close all;

a = imread('osteo2.jpg');
b = imcrop(a);
P = impixel(b);
close;

[m,n,p] = size(b);

Pav = mean(P,1);
imf = (im2double(reshape(b, m*n, 3))).*255;

cutoff = .3; % standard deviations

q = repmat(Pav,[size(imf,1),1]);

% Find the difference from the reference for each [R,G,B]
distance = sum((imf - q).^2, 2);

% Figure out cases where the distance is within the cutoff
tf = distance <= cutoff.*std(distance);

tra = reshape(tf, m, n);

figure(2);
imshow(b);

figure(3);
imshow(tra);
%
% minRegionSize = 10
%
% [lbl, num] = bwlabel(tra, 8);
%
% for i=1:max(lbl)
%     if sum (lbl==i) < minRegionSize
%         lbl(lbl==i) = 0;
%     end
% end
%
% num

% minRegionSize = 4;
% connStruct = bwconncomp (tra);
% regions = connStruct.PixelIdxList;
%
% ct = 0;
%
% for i=1:connStruct.NumObjects
%     if numel (regions{i-ct}) < minRegionSize
%         regions{i-ct} = []
%         ct = ct+1;
%     end
% end