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

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
Previous Post Next Post Home
Google ping Hypersmash.com