Color Clustering Matlab

by

 

28 Jan 2012 (Updated )

Computes the clusters of Pixels based upon their color.

binary_quant.m
%%Computes the clusters of Pixels based upon their color.
%% The algorithm is based uopon binary tree quantization technique
%% described by Orchard and Bouman.

%% input arguments
%%Image = input image should be RGB.
%% K= Maximum possible clusters required.

%% output arguments
%% C= A cell array containing cluster structures. Each cluster has
%% following fields:
%% point = An integer array containing linear indices of points in
%% clusters.
%% M = Sum of RGB values of all the pixels in cluster
%% N=  Number of Pixels in the cluster
%% R=  The sum of squares of pixels in the cluster
%% vector= Principal eigen vector of cluster

%% References:
%% ORCHARD, M. T., AND BOUMAN, C. A. 1991. Color Quantiza-tion of Images. 
%%IEEE Transactions on Signal Processing 39, 12,26772690.



function C=binary_quant(Image,K)
%% Read the image and convert it into double.
I=imread(Image);
I=im2double(I);
si=1;
R_Comp=I(:,:,1);
G_Comp=I(:,:,2);
B_Comp=I(:,:,3);
%% Treat complete image as single cluster and compute different fields of
%% cluster.
R1=zeros(3,3);
R1(1,1)=sum(sum(R_Comp.*R_Comp));
R1(1,2)=sum(sum(R_Comp.*G_Comp));
R1(1,3)=sum(sum(R_Comp.*B_Comp));
R1(2,1)=sum(sum(G_Comp.*R_Comp));
R1(2,2)=sum(sum(G_Comp.*G_Comp));
R1(2,3)=sum(sum(G_Comp.*B_Comp));
R1(3,1)=sum(sum(B_Comp.*R_Comp));
R1(3,2)=sum(sum(B_Comp.*G_Comp));
R1(3,3)=sum(sum(B_Comp.*B_Comp));
R_Comp=R_Comp(:);
G_Comp=G_Comp(:);
B_Comp=B_Comp(:);
N1=size(I(:),1)/3;
M=sum(sum(I));
M1=zeros(1,3);
M1(1,1)=M(:,:,1);
M1(1,2)=M(:,:,2);
M1(1,3)=M(:,:,3);
%% Compute covariance of the Image.
R1_bar=R1-((M1*M1')/N1);
%% Compute eigen values and eigen vectors of the image.
[V,D]=eig(R1_bar);
point=1:1:size(I,1)*size(I,2);
%% Form the initial cluster.
cluster.R=R1;
cluster.M=M1;
cluster.N=N1;
cluster.point=point;
cluster.vector=V(:,1);
%% Place the cluster in Cell array.
C{si}=cluster;
%% Store eigen value in a seperate array.
val(si)=abs(D(1,1));
for i=2:1:K
    if(max(val(:))==-1)
        break;
    end
    %% Pick the cluster having highest eigen value and split it into two parts.
  index=find(max(val(:))==val);
  vector=C{index}.vector;
  point=C{index}.point;
  point=point(:);
  M=C{index}.M;
  N=C{index}.N;
  Q=M/N;
  R=C{index}.R;
  R1=zeros(3,3);
  M1=zeros(1,3);
         a= vector(1,1)*R_Comp(point(:)) + vector(2,1)*G_Comp(point(:)) + vector(3,1)*B_Comp(point(:));
         %% Pixels of the cluster are divided into two parts based on their
         %% closeness to a plane perpendicular to principal eigen vector
         %% and passing through mean.
         point1=point(find(a<=vector'*Q'));
         point2=point(find(a>vector'*Q'));
         %% Compute various statistical value for the newly formed
         %% clusters.
              R1(1,1)=sum(R_Comp(point1(:)).*R_Comp(point1(:)));
              R1(1,2)=sum(R_Comp(point1(:)).*G_Comp(point1(:)));
              R1(1,3)=sum(R_Comp(point1(:)).*B_Comp(point1(:)));
              R1(2,1)=sum(G_Comp(point1(:)).*R_Comp(point1(:)));
              R1(2,2)=sum(G_Comp(point1(:)).*G_Comp(point1(:)));
              R1(2,3)=sum(G_Comp(point1(:)).*B_Comp(point1(:)));
              R1(3,1)=sum(B_Comp(point1(:)).*R_Comp(point1(:)));
              R1(3,2)=sum(B_Comp(point1(:)).*G_Comp(point1(:)));
              R1(3,3)=sum(B_Comp(point1(:)).*B_Comp(point1(:)));
              M1(1,1)=sum(R_Comp(point1(:)));
              M1(1,2)=sum(G_Comp(point1(:)));
              M1(1,3)=sum(B_Comp(point1(:)));
  R2=R-R1;
  M2=M-M1;
  N1=size(point1,1);
  N2=size(point2,1);
  %% Place the clusters into cell array and delete the parent cluster.
  if(N1>0 && N2>0)
  cluster1.point=point1;
  cluster1.R=R1;
  cluster1.M=M1;
  cluster1.N=N1;
  R1_bar=R1-((M1*M1')/N1);
  [V1 D1]=eig(R1_bar);
  cluster1.vector=V1(:,1);
  cluster2.R=R2;
  cluster2.N=N2;
  cluster2.M=M2;
  cluster2.point=point2;
  R2_bar=R2-((M2*M2')/N2);
  [V2 D2]=eig(R2_bar);
  cluster2.vector=V2(:,1);
  val(si+1)=abs(D2(1,1));
  val(index)=abs(D1(1,1));
  C{index}=cluster1;
  C{si+1}=cluster2;
  si=si+1;
  else
     val(index)=-1;
     i=i-1;
  end
end
%% display original Image.
figure,imshow(I);
%% display newly formed clusters.
for i=1:1:si;
    temp=C{i}.point;
    I1=zeros(size(I,1),size(I,2),3);
       [y_p x_p]=ind2sub([size(I,1) size(I,2)],temp);
  for t=1:1:size(y_p,1)
           y=y_p(t);
           x=x_p(t);
        I1(y,x,1)=I(y,x,1);
        I1(y,x,2)=I(y,x,2);
        I1(y,x,3)=I(y,x,3);
    end
    figure,imshow(I1);
end
end 

Contact us