%vs_pattern1
%% identify objects by pattern recognition
% using euclidian and Mahalanobis distance

A=im2double(imread('orgimage.png'));
 
[L N]=bwlabel(A);
s=regionprops(L,'area');

%% discard small objects. Create a simplified image.
areas=[s.Area];
idx=find(areas > 200);  %200 is indicative; adjust with others
F=ismember(L,idx); %Give a binary matrix based on c only where vaues from idx exist
B=A;
B(~F)=0;

%% Recalculate interesting properties for selected objects
[L N]=bwlabel(B);
s=regionprops(L,'area','centroid','eccentricity', 'solidity', 'perimeter');
are=[s.Area];
ecc=[s.Eccentricity];
sol=[s.Solidity];
per=[s.Perimeter];
ffa=4*pi*are./per./per;
% create the patterns
Patts=[ffa' sol' ecc']';

%% Plot the patterns just to see
imshow(B); hold on;
for n=1:N
    mstr= num2str( Patts(:,n), 3);
    text( s(n).Centroid(1)+10, s(n).Centroid(2), mstr, 'Color', [1 0 1]);
    text( s(n).Centroid(1)-30, s(n).Centroid(2), num2str(n), 'Color', [0 1 0], 'BackgroundColor',[0 0 0]);
end

%% now check all againts two specif objects know to be a category each
pA=Patts(:,1);  %%object A (pá)
pB=Patts(:,4);  %%object B (serrote)
dA=[]; dB=[];
for n=1:N
    dA = [dA norm(Patts(:,n)-pA)]; %#ok<*AGROW>
    dB = [dB norm(Patts(:,n)-pB)];
end

% just for display
Res=[Patts' dA' dB'];

%% Now use mahalanobis distance
% must perform some stats analysis before

samplesA=[1 12 14 16 17 18];  %all pás
samplesB=[4 6 19];   % all serrotes


PattsA=Patts(:, samplesA);
PattsB=Patts(:, samplesB);

ssA=cov(PattsA');
ssB=cov(PattsB');

% ----
% dAmm=[]; dBmm=[];
% %Use next line to use as reference the mean of all A and B patterns
% % or comment to use only one of them as reference
% pA=mean(PattsA,2); pB=mean(PattsA,2);
% 
% for n=1:N
%      dAmm = [dAmm norm((Patts(:,n)-pA)'*inv(ssA)*(Patts(:,n)-pA) )]; %#ok<*AGROW>
%      dBmm = [dBmm norm((Patts(:,n)-pB)'*inv(ssB)*(Patts(:,n)-pB) )];
% end

Y=Patts';
X=PattsA';
dAm = MAHAL(Y,X) ;%returns the Mahalanobis distance (in squared units) of each observation (point) in Y from the sample data in X, i.e.,

X=PattsB';
dBm = MAHAL(Y,X) ;
%% -------------


% just for display
Res=[Patts' dA' dB' dAm/max(max(dAm)) dBm/max(max(dBm)) ]

