Lets Learn together... Happy Reading

" Two roads diverged in a wood, and I,
I took the one less traveled by,
And that has made all the difference "-Robert Frost

K means clustering on RGB image

                        I assume the readers of this post have enough knowledge on K means clustering method and it’s not going to take much of your time to revisit it again.


Let’s start with a simple example, consider a RGB image as shown below.


Let’s choose the number of clusters = 2. Basically we are going to separate the background (first cluster) and the flower (second cluster).

In the second step, let’s choose two random RGB pixel values.

Since the initial pixel values are completely random, we can clearly see that both the initial cluster points are on the flower and not on the background.
The pixel value of the cluster 1 (R1,G1,B1)=[255,191,59]
The pixel value of the cluster 2 (R2,G2,B2) = [255,245,11]

In the third step, find the Euclidean distance between the initial points and all the pixels in the image matrix.
We will see how to calculate the Euclidean distance between the initial points and the pixel value at (1,1) . The pixel value at (1,1) is [24,64,186]





Repeat the above procedure for all the pixel values.

The next step is finding the minimum value among the two clusters and fetch their corresponding cluster index. With reference to the pixel position at (1,1), the minimum value is 292.6072 and it belongs to the cluster 1. So update the matrix ‘ClusterMap’ with 1 at the position (1,1).
Similarly, find the index of the minimum value and update the ‘ClusterMap’ matrix for all the pixel positions.



The visual analysis of the first iteration shows that the clustering is not good and it still needs to be improved.
It’s time to look for new cluster points. Let’s see how to find new cluster points to perform better segmentation.
In order to obtain the new cluster points, compute the mean of the pixel values with respect to each cluster. Once we obtain the mean value for the Red, Green and Blue channels separately, find the difference between the new values and the current values.
In our chosen example,
The new values for the Red, Green and Blue channels are (R1_new,G1_new,B1_new)=[ 93.2942  106.6650  143.5914] for the first cluster
And (R2_new,G2_new,B2_new)=[ 252.1034  224.2628    3.8519] for the second cluster

The difference between the new and the current values with respect to different color channels are
[  161.7058   84.3350   84.5914] for first cluster (R1-R1_new,G1-G1_new,B1-B1_new)
 [  2.8966   20.7372    7.1481] for the second cluster (R2-R2_new,G2-G2_new,B2-B2_new)




If the difference is more than the threshold value then assign the new values to the current values and continue from the third step i.e calculating the Euclidean distance else stop the process and display the clustered image.

In our example, after the fifth iteration the difference between the new and the current values of the
  1st cluster [0.1615    0.1734    0.5675] and 2nd cluster [0.4701    0.5028    0.0364]  indicates that the values are less than 1 for each color channel individually.
The final segmented images are shown below.
%K means clustering
%Author: Aaron Angel

clear all
close all
clc
%READ A RGB IMAGE
A = imread('flower8.jpg');
figure,subplot(121),imagesc(A);title('RGB Image');hold on;

A = double(A);

%CARTESIAN GRID FOR THE 2D IMAGE
[Xp,Yp] = meshgrid(1:size(A,2),1:size(A,1));

%NUMBER OF CLUSTERS
num_clusters = 4;

%THRESHOLD VALUE FOR EACH CHANNEL
Tval = 1;

%GLOBAL THRESHOLD VALUE
Global_Tval = 3;

%RANDOM IMAGE POSITIONS
RX = randi(size(A,1),1,num_clusters);
RY = randi(size(A,2),1,num_clusters);



%PREALLOCATE THE VECTORS
RGB_val = zeros([num_clusters,3]);
RGB_new = zeros([num_clusters,3]);

%FETCH THE RGB PIXEL VALUE OF THE RANDOM IMAGE POSITIONS
for j = 1:num_clusters
   RGB_val(j,:) = A(RX(j),RY(j),:); 
end


myvox = zeros([size(A,1) size(A,2) num_clusters]);
flag = 1;

Rcomp = A(:,:,1); %RED CHANNEL
Gcomp = A(:,:,2); %GREEN CHANNEL
Bcomp = A(:,:,3); %BLUE CHANNEL
it=0;


while flag==1
 
 %EUCLIDEAN DISTANCE BETWEEN THE RANDOM RGB PIXEL VALUE
 %AND THE PIXEL VALUES IN THE IMAGE
 for j = 1: num_clusters
     myvox(:,:,j) = sqrt((A(:,:,1)-RGB_val(j,1)).^2+(A(:,:,2)-RGB_val(j,2)).^2+(A(:,:,3)-RGB_val(j,3)).^2);
 end

%FIND THE MINIMUM VALUE(Y) AND THE CORRESPONDING CLUSTER NUMBERS(ClusterMap)
[Y,ClusterMap] = min(myvox,[],3);


%MEAN VALUE PIXEL VALUES IN EACH CHANNEL WITH RESPECT TO THE CLUSTERS
for j = 1:num_clusters
  
   RGB_new(j,1) = mean(Rcomp(ClusterMap(:)==j));
   RGB_new(j,2) = mean(Gcomp(ClusterMap(:)==j));
   RGB_new(j,3) = mean(Bcomp(ClusterMap(:)==j));
end

%DIFFERENCE BETWEEN THE NEW VALUE AND THE OLD VALUE
DiffVal = abs(RGB_val-RGB_new);


%IF THE DIFFERENCE IS LESS,STOP THE ITERATION ELSE
%ASSIGN THE NEW VALUE AND CONTINUE
if(sum(DiffVal(:)) < Global_Tval)
    flag = 0;
else
  if(sum(DiffVal(:,1))>Tval)
    RGB_val(:,1) = RGB_new(:,1);
  end
  if(sum(DiffVal(:,2))>Tval)
    RGB_val(:,2) = RGB_new(:,2);
  end
  if(sum(DiffVal(:,3))>Tval)
     RGB_val(:,3) = RGB_new(:,3);
  end
end




%NUMBER OF ITERATIONS
it=it+1

end;
subplot(122),imagesc(ClusterMap);title('Clusters');colormap(jet);


%DISPLAY THE SEGMENTED IMAGE BASED ON COLORS
m=2;
n = round(num_clusters/m);
for k=1:num_clusters
F = repmat(logical(ClusterMap==k),1,1,3).*double(A);
figure(2),subplot(n,m,k),imagesc(uint8(F));hold on;
end



EXPLANATION:


In this example, the number of clusters are 4, so the yellow, pink and blue flowers are defined as clusters separately while the remaining details are combined as one cluster.



Hope you all enjoyed the post. Feel free to comment and join us in Facebook and twitter for more updates.
like button Like "IMAGE PROCESSING" page

Principal Component Analysis


Principal component analysis is a statistical technique that is used in finding patterns and reducing the dimensions of multi-dimensional data. There is an excellent tutorial by Lindsay I Smith on this topic so I will be focusing more on the application part in this post.




Let’s start with the image compression application in this post.

EXAMPLE:



MATLAB CODE:
clear all
clc

%READ A GRAYSCALE IMAGE
I = imread('cameraman.tif');

figure(1),subplot(121),imshow(I);title('INPUT IMAGE');

A = double(I);

%MEAN OF EACH COLUMN IN THE IMAGE
meanCols = mean(A,1);

%SUBTRACT THE MEAN VALUE OF EACH COLUMN
for k = 1:size(A,1)
    A(k,:) = A(k,:)-meanCols;
end


%COMPUTE COVARIANCE MATRIX
covmat = cov(A);
%OBTAIN EIGEN VALUES
[coeff,D] = eig(covmat);
coeff = fliplr(coeff);
figure(2),imagesc(coeff);colormap(jet);
FV = coeff(:,1:50)'; %PRINCIPAL COMPONENT
Res = FV*A';





%RECONSTRUCTION
Org = (FV'*Res)';

%ADD THE MEAN VALUES OF THE COLUMN
Output = zeros(size(A));
for k = 1:size(A,1)
    Output(k,:) = Org(k,:)+meanCols;
end

Output = uint8(Output);
figure(1),subplot(122),imshow(Output);title('RESTORED IMAGE');







EXAMPLE 2:

Let’s consider the LANDSAT-8 dataset. I have chosen the following bands.
Band 2– BLUE,
Band 3 – Green,
Band 4 – Red,
Band 5 – NIR
Band 6 - SWIR-1
Band 7 – SWIR-2









MATLAB CODE:
clear all

% READ ALL THE TIF IMAGE NAMES
fnames = dir('*.TIF');
dim1=400;
dim2=950;
C=[];
inc=1;

%CREATE A VOXEL OF IMAGES
for k = [4,5,6,7,8,9];

     A = imread(fnames(k).name);
     A = A(1901:2300,2251:3300);
     inc=inc+1;
     C =[C,A(:)];
 end

 C1 = double(C);
 %SUBTRACT THE MEAN
 meanval = zeros([1,size(C,2)]);
 for k = 1:size(C,2)
     meanval(1,k) = mean(squeeze(C(:,k)));
    C1(:,k) = C1(:,k)- meanval(1,k);
 end

%PERFORM PCA ANALYSIS
coeff = pca(C1,'Algorithm','eig');
figure,imagesc(coeff);colormap(jet);
FV=coeff(:,[1:3])'; %PRINCIPAL COMPONENTS


Res = FV*C1';


%RECONSTRUCTION
Org = (FV'*Res)';

%2-Blue, 3-Green 4-Red, 5-NIR, 6 and 7 SWIR-1 and 2
%NATURAL COLOR
%ADD THE MEAN VALUE
R1 = Org(:,3)+mean(C(:,3));
G1 = Org(:,2)+mean(C(:,2));
B1 = Org(:,1)+mean(C(:,1));

R1= reshape(R1,[size(A,1) size(A,2)]);
G1= reshape(G1,[size(A,1) size(A,2)]);
B1= reshape(B1,[size(A,1) size(A,2)]);

%ENHANCE THE CONTRAST
R1 = BCET16(0,65535,20000,R1);
G1 = BCET16(0,65535,20000,G1);
B1 = BCET16(0,65535,20000,B1);

Output = zeros([size(A,1) size(A,2) 3]);
Output(:,:,1)=R1;
Output(:,:,2)=G1;
Output(:,:,3)=B1;

figure(8),imagesc(uint16(Output));title('Natural Color');



%6-SWIR-2,4-NIR,2-GREEN
%ADD THE MEAN VALUE
R1 = Org(:,6)+mean(C(:,6));
G1 = Org(:,4)+mean(C(:,4));
B1 = Org(:,2)+mean(C(:,2));

R1= reshape(R1,[size(A,1) size(A,2)]);
G1= reshape(G1,[size(A,1) size(A,2)]);
B1= reshape(B1,[size(A,1) size(A,2)]);

%ENHANCE THE CONTRAST
R1 = BCET16(0,65535,25000,R1);
G1 = BCET16(0,65535,25000,G1);
B1 = BCET16(0,65535,25000,B1);

Output1 = zeros([size(A,1) size(A,2) 3]);
Output1(:,:,1)=R1;
Output1(:,:,2)=G1;
Output1(:,:,3)=B1;
figure(4),imagesc(uint16(Output1));title('Vegetation and Water');



%SWIR-1, SWIR-2,Red
%Urban
%ADD THE MEAN VALUE
R1 = Org(:,6)+mean(C(:,6));
G1 = Org(:,5)+mean(C(:,5));
B1 = Org(:,3)+mean(C(:,3));

R1= reshape(R1,[size(A,1) size(A,2)]);
G1= reshape(G1,[size(A,1) size(A,2)]);
B1= reshape(B1,[size(A,1) size(A,2)]);

%ENHANCE THE CONTRAST
R1 = BCET16(0,65535,25000,R1);
G1 = BCET16(0,65535,25000,G1);
B1 = BCET16(0,65535,25000,B1);

Output2 = zeros([size(A,1) size(A,2) 3]);
Output2(:,:,1)=R1;
Output2(:,:,2)=G1;
Output2(:,:,3)=B1;
figure(5),imagesc(uint16(Output2));title('Urban Environment');



%SWIR-1,NIR,RED(4,3,2)
%Color IR
%ADD THE MEAN VALUE
R1 = Org(:,4)+mean(C(:,4));
G1 = Org(:,3)+mean(C(:,3));
B1 = Org(:,2)+mean(C(:,2));

R1= reshape(R1,[size(A,1) size(A,2)]);
G1= reshape(G1,[size(A,1) size(A,2)]);
B1= reshape(B1,[size(A,1) size(A,2)]);

%ENHANCE THE CONTRAST
R1 = BCET16(0,65535,25000,R1);
G1 = BCET16(0,65535,25000,G1);
B1 = BCET16(0,65535,25000,B1);

Output3 = zeros([size(A,1) size(A,2) 3]);
Output3(:,:,1)=R1;
Output3(:,:,2)=G1;
Output3(:,:,3)=B1;
figure(6),imagesc(uint16(Output3));title('Color IR');


EXPLANATION:

The total number of bands chosen is 6. PC1, PC2 and PC3 are used to reconstruct the 6 bands. PC1 is the largest contributor. By combining different bands together we obtain natural color, urban environment, vegetation and color IR.

NOTE: The dataset is obtained from https://landsatlook.usgs.gov/

The function BCET16 is used to enhance the contrast. Save the below function as BCET16.m and make sure the m file is available in the current working directory.


function y=BCET16(Gmin,Gmax,Gmean,x)
    %BALANCE CONTRAST ENHANCEMENT TECHNIQUE
    x    = double(x); % INPUT IMAGE
    Lmin = min(x(:)); % MINIMUM OF INPUT IMAGE
    Lmax = max(x(:)); % MAXIMUM OF INPUT IMAGE
    Lmean = mean(x(:)); %MEAN OF INPUT IMAGE
    LMssum = mean(x(:).^2); %MEAN SQUARE SUM OF INPUT IMAGE



    bnum = Lmax.^2*(Gmean-Gmin) - LMssum*(Gmax-Gmin) + Lmin.^2*(Gmax-Gmean);
    bden = 2*(Lmax*(Gmean-Gmin)-Lmean*(Gmax-Gmin)+Lmin*(Gmax-Gmean));

    b = bnum/bden;

    a = (Gmax-Gmin)/((Lmax-Lmin)*(Lmax+Lmin-2*b));

    c = Gmin - a*(Lmin-b).^2;

    y = a*(x-b).^2+c; %PARABOLIC FUNCTION
    y = uint16(y);
end


To learn more about BCET check, Balance Contrast Enhancement Technique
like button Like "IMAGE PROCESSING" page
Next Post Home
Google ping Hypersmash.com