%This is a script to demo the statistical image feature
%Author: Jin, Tao

clear, close all

% 1. read and display rice image
I=imread('rice.tif');
figure('Name', 'Original image'), imshow(I)
I=im2double(I);

% perform block processing to approximate the background
bkg=blkproc(I, [32,32], 'min(x(:))');
bkg=imresize(bkg, [256,256], 'biliniear');

% subtract the background from the original image
I2 = I - bkg;
I2 = max(min(I2,1),0);

% adjust image constrast
I3 = imadjust(I2, [0 max(I2(:))], [0 1]);

% apply thresholding to the image
bw = I3>0.2;

% label components
[L, numObj] = bwlabel(bw, 4); 

map = hot(numObj+1); %create a colormap
figure('Name', 'labeled binary image'), imshow(L+1, map);

% Compute feature measurements of objects in the image
FeatureData = imfeature(L, 'all');

%perimeter
P=bwperim(L);
figure('Name','perimeter'), imshow(P)

%Area 
Array = [FeatureData.Area];

Min = min(Array);
Max = max(Array);
Mean = mean(Array);
Std = std(Array);

figure('Name', 'histogram for area'), hist(Array, 20);

idx = find(Array < (Mean-Std));
L2 = ismember(L, idx);
L2 = I.*L2;
figure('Name', 'Objects below mean-std :area'), imshow(L2);

idx = find(Array > (Mean+Std));
L2 = ismember(L, idx);
L2 = I.*L2;
figure('Name', 'Objects above mean+std :area'), imshow(L2);

idx = find((Array < (Mean+Std)) & (Array > (Mean-Std)));
L2 = ismember(L, idx);
L2 = I.*L2;
figure('Name', 'Objects between mean-std and mean+std :area'), imshow(L2);

% %Second moments
% %eccentricity, orientation, MajorAxisLength, MinorAxisLength
Array = [FeatureData.Orientation];

idx = find((Array < 0) & (Array > -90));
L2 = ismember(L, idx);
L2 = I.*L2;
figure('Name', 'Objects between -90 to 0 degree'), imshow(L2);

idx = find((Array > 0) & (Array < 90));
L2 = ismember(L, idx);
L2 = I.*L2;
figure('Name', 'Objects between 0 and 90 degree'), imshow(L2);

idx = find((Array < 180) & (Array > 90));
L2 = ismember(L, idx);
L2 = I.*L2;
figure('Name', 'Objects between 90 and 180 degree'), imshow(L2);


% Min = min(Array);
% Max = max(Array);
% Mean = mean(Array);
% Std = std(Array);
% 
% figure('Name', 'histogram for Orientation'), hist(Array, 20);
% 
% idx = find(Array < (Mean-Std));
% L2 = ismember(L, idx);
% L2 = I.*L2;
% figure('Name', 'Objects below mean-std :Orientation'), imshow(L2);
% 
% idx = find(Array > (Mean+Std));
% L2 = ismember(L, idx);
% L2 = I.*L2;
% figure('Name', 'Objects above mean+std :Orientation'), imshow(L2);
% 
% idx = find((Array < (Mean+Std)) & (Array > (Mean-Std)));
% L2 = ismember(L, idx);
% L2 = I.*L2;
% figure('Name', 'Objects between mean-std and mean+std :Orientation'), imshow(L2);
% 
% %
% Array = [FeatureData.Eccentricity];
% 
% Min = min(Array);
% Max = max(Array);
% Mean = mean(Array);
% Std = std(Array);
% 
% figure('Name', 'histogram for Eccentricity'), hist(Array, 20);
% 
% idx = find(Array < (Mean-Std));
% L2 = ismember(L, idx);
% L2 = I.*L2;
% figure('Name', 'Objects below mean-std :Eccentricity'), imshow(L2);
% 
% idx = find(Array > (Mean+Std));
% L2 = ismember(L, idx);
% L2 = I.*L2;
% figure('Name', 'Objects above mean+std :Eccentricity'), imshow(L2);
% 
% idx = find((Array < (Mean+Std)) & (Array > (Mean-Std)));
% L2 = ismember(L, idx);
% L2 = I.*L2;
% figure('Name', 'Objects between mean-std and mean+std :Eccentricity'), imshow(L2);
% 
% %
% Array = [FeatureData.MajorAxisLength];
% 
% Min = min(Array);
% Max = max(Array);
% Mean = mean(Array);
% Std = std(Array);
% 
% figure('Name', 'histogram for MajorAxisLength'), hist(Array, 20);
% 
% idx = find(Array < (Mean-Std));
% L2 = ismember(L, idx);
% L2 = I.*L2;
% figure('Name', 'Objects below mean-std :MajorAxisLength'), imshow(L2);
% 
% idx = find(Array > (Mean+Std));
% L2 = ismember(L, idx);
% L2 = I.*L2;
% figure('Name', 'Objects above mean+std :MajorAxisLength'), imshow(L2);
% 
% idx = find((Array < (Mean+Std)) & (Array > (Mean-Std)));
% L2 = ismember(L, idx);
% L2 = I.*L2;
% figure('Name', 'Objects between mean-std and mean+std :MajorAxisLength'), imshow(L2);
% 
% %
% Array = [FeatureData.MinorAxisLength];
% 
% Min = min(Array);
% Max = max(Array);
% Mean = mean(Array);
% Std = std(Array);
% 
% figure('Name', 'histogram for MinorAxisLength'), hist(Array, 20);
% 
% idx = find(Array < (Mean-Std));
% L2 = ismember(L, idx);
% L2 = I.*L2;
% figure('Name', 'Objects below mean-std :MinorAxisLength'), imshow(L2);
% 
% idx = find(Array > (Mean+Std));
% L2 = ismember(L, idx);
% L2 = I.*L2;
% figure('Name', 'Objects above mean+std :MinorAxisLength'), imshow(L2);
% 
% idx = find((Array < (Mean+Std)) & (Array > (Mean-Std)));
% L2 = ismember(L, idx);
% L2 = I.*L2;
% figure('Name', 'Objects between mean-std and mean+std :MinorAxisLength'), imshow(L2);
% 
% Centroid
Array = [FeatureData.Centroid];
Array = reshape(Array, 2, numObj);

Min =  [min(Array(1,:));  min(Array(2,:))];
Max =  [max(Array(1,:));  max(Array(2,:))];
Mean = [mean(Array(1,:)); mean(Array(2,:))];
Std =  [std(Array(1,:));  std(Array(2,:))];

figure('Name', 'histogram for Centroid.row'), hist(Array(1,:), 20);
figure('Name', 'histogram for Centroid.col'), hist(Array(2,:), 20);

[row, col] = find( (Array(1,:) < (Mean(1,1)-Std(1,1))) & (Array(2,:) < (Mean(2,1)-Std(2,1))));
L2 = ismember(L, col);
L2 = I.*L2;
figure('Name', 'Objects below mean-std :Centroid'), imshow(L2);

[row, col] = find( (Array(1,:) > (Mean(1,1)+Std(1,1))) & (Array(2,:) > (Mean(2,1)+Std(2,1))));
L2 = ismember(L, col);
L2 = I.*L2;
figure('Name', 'Objects above mean+std :Centroid'), imshow(L2);

[row, col] = find( (Array(1,:) > (Mean(1,1)-Std(1,1))) & (Array(2,:) > (Mean(2,1)-Std(2,1))) & ...
    (Array(1,:) < (Mean(1,1)+Std(1,1))) & (Array(2,:) < (Mean(2,1)+Std(2,1))));
L2 = ismember(L, col);
L2 = I.*L2;
figure('Name', 'Objects between mean-std and mean+std :Centroid'), imshow(L2);