Hyperspectra Image Analysis

HSI

 clc
 clear

 
 int lines
 int sample
 int bands
 lines=400;  
 sample=502;
 bands=128;
int a;
int k;
int b;
int band_selected;

a=1;
k=1;
b=1;
band_selected=1; %read pics from a specific band


int min_gray_value_of_a_pixel;
int max_gray_value_of_a_pixel;
min_gray_value_of_a_pixel=1;
max_gray_value_of_a_pixel=255;

% cwd = pwd;
% cd(tempdir);
% pack
% cd(cwd)
 
%%%%%%%Read in data from HIS's file%%%%%%%%%%%%%%% 

x=multibandread('C:\Documents and Settings\Administrator\Desktop\Hyperspectra_Zhang\July_23rd_2013\sample with test soil from home depot\310B-4_310B-7.IMG',[lines,sample,bands], 'uint16',0,'bil','ieee-le');
%%%%%%pre-scaling process, change integer numbers into floating numbers for display%%%%%
%figure(5) ;colormap(gray);
%imagesc(x(:,:,band_selected))
x=double(x);
%%%%%rescaling of the value of pixels in to range of 1 to 255%%%%%%%%%%%%%
min_gray_value_of_a_pixel=min(min(min(x)));
max_gray_value_of_a_pixel=max(max(max(x)));
x=x-min_gray_value_of_a_pixel;
x=x*255/max_gray_value_of_a_pixel;



%%%%Show the image and select ROI%%%%%%%%%%%
% band_selected %%%read from selected band 

image(:,:) = x(:,:,band_selected)/mean(mean(x(:,:,band_selected)));

figure(1)
imshow(image(1:1:lines, 1:1:sample))


p=impoint(gca,[])
api=iptgetapi(p)
p1=api.getPosition()
length=10;
p1=floor(p1);
p2(1,1)=p1(1,1)-length/2;
p2(1,2)=p1(1,2)-length/2;
h=imrect(gca,[p2(1,1) p2(1,2) length length]);
p3(1,1)=p2(1,1)+length;
p3(1,2)=p2(1,2)+length;



% ROI=image(p2(1,1):p3(1,1),p2(1,2):p3(1,2));
% figure(2), imshow(ROI); 
% api=iptgetapi(h);





 y=multibandread('C:\Documents and Settings\Administrator\Desktop\Hyperspectra_Zhang\June_28th_2013 (home deport 100 samples with test soil)\140B-4.img',[lines,sample,bands], 'uint16',0,'bil','ieee-le');
     
%%%%%%%%%%%%%%%%%Average%%%%%%%%%%%%%%%%%%%%%%%%


for a=1:1:bands
     background(a)=mean(mean(y(:,:,a)));
     spectra_raw_data(a)=mean(mean(x(p2(1,2):1:p3(1,2),p2(1,1):1:p3(1,1)),a));
     spectra_normalized_agianst_background(a)=mean(mean(x(p2(1,2):1:p3(1,2),p2(1,1):1:p3(1,1),a)))/mean(mean(y(:,:,a)));
end

figure(2)

subplot(2,3,1),  plot([396.600000000000:4.700000000000:995.500000000000], background)
 

subplot(2,3,2), plot([396.600000000000:4.700000000000:995.500000000000], spectra_raw_data)


subplot(2,3,3), plot([396.600000000000:4.700000000000:995.500000000000], spectra_normalized_agianst_background)




%%%%%%%%SVD%%%%%%%%%%%%%%%spectra_svd%%%%%%


%tol = max (size (x(p2(1,2):1:p3(1,2),:,p2(1,1):1:p3(1,1)))) * sigma_max (x(p2(1,2):1:p3(1,2),:,p2(1,1):1:p3(1,1))) * eps ;
[U,S,V]=svd(pinv(x(p2(1,2):1:p3(1,2),p2(1,1):1:p3(1,1),:))/sqrt(bands-1)); % perform the SVD
lambda=diag(S).^2; %produce diagonal variance 
%tol = max (size (X)) * sigma_max (X) * eps ; 
spectra_svd=pinv(U)*x(p2(1,2):1:p3(1,2),p2(1,1):1:p3(1,1),:); % produce the principal component projection 
subplot(2,3,4), plot([396.600000000000:4.700000000000:995.500000000000],spectra_svd)



%%%%%%%%%PCA%%%%%%%%%%%%%%%% 

k=128; %%%%%this has to be selected for certain contaminations after study 


for a=1:1:bands
     
     spectra_pca(a)=mean(mean(x(p2(1,2):1:p3(1,2),p2(1,1):1:p3(1,1),a)))/mean(mean(x(:,:,a)));
     
end

Sigma=spectra_pca'*spectra_pca;
[U,S,V]=svd(Sigma);
U_reduce(:,1:k)=U(:,1:k);
spectra_pca=pinv(U_reduce)*pinv(spectra_pca);

 
subplot(2,3,5), plot([396.600000000000:4.700000000000:995.500000000000],spectra_pca)

%%%%%%%%%%ICA, define the weight for each component first [1 0; 0,1] for test soil and clean background%%%%%%%%%%%%%%%%%

%original image show up, one is  test soil, the other is colorful background
%with test soil
IMG_1=multibandread('C:\Documents and Settings\Administrator\Desktop\Hyperspectra_Zhang\June_28th_2013 (home deport 100 samples with test soil)\140B-4.img',[lines,sample,bands], 'uint16',0,'bil','ieee-le'); 
%colorful background
IMG_2=multibandread('C:\Documents and Settings\Administrator\Desktop\Hyperspectra_Zhang\June_28th_2013 (home deport 100 samples with test soil)\140B-4.img',[lines,sample,bands], 'uint16',0,'bil','ieee-le'); 

IMG_1=double(IMG_1);
IMG_2=double(IMG_2);

for a=1:1:bands
     S1(a)=mean(mean(IMG_1(:,a,:))); %spectra of background with test soil
     S2(a)=mean(mean(IMG_2(:,a,:))); %spectra of clean colorful background
end


figure(3)
subplot(3,2,1), plot([396.600000000000:4.700000000000:995.500000000000],S1)
subplot(3,2,2), plot([396.600000000000:4.700000000000:995.500000000000],S2)

%Mixing Matrics: how sensitive specrtra of images to the data
beta1=1;
beta2=0;
beta3=0;
beta4=1;

A=[beta1 beta2;beta3 beta4];
S1=double(A(1,1)*S1+A(1,2)*S2);
S2=double(A(2,1)*S1+A(2,2)*S2);

subplot(3,2,3), plot(396.600000000000:4.700000000000:995.500000000000,S1)
subplot(3,2,4), plot(396.600000000000:4.700000000000:995.500000000000,S2)

%[m,n]=size(X1);
% X1=reshape(X1,m*n,1);
% X2=reshape(X2,m*n,1);

theta0=0.5*atan(-2*sum(S1.*S2)/sum(S1.^2-S2.^2))
%rotation matrix
Us=[cos(theta0) sin(theta0); -sin(theta0) cos(theta0)]  
Sig1=sum((S1*cos(theta0)+S2*sin(theta0)).^2)
Sig2=sum((S1*cos(theta0-pi/2)+S2*sin(theta0-pi/2)).^2)
Sigma_ICA=[1/sqrt(Sig1) 0; 0 1/sqrt(Sig2)];

%X1bar=Sigma_ICA(1,1)*(Us(1,1)*X1+Us(1,2)*X2);
%X2bar=Sigma_ICA(2,2)*(Us(2,1)*X1+Us(2,2)*X2);
%X1bar=reshape(X1bar,m*n,1);
%X2bar=reshape(X2bar,m*n,1);

S1bar=Sigma_ICA(1,1)*(Us(1,1)*S1+Us(1,2)*S2);
S2bar=Sigma_ICA(2,2)*(Us(2,1)*S1+Us(2,2)*S2);

Phi0=0.25*atan(-sum(S1bar.^3).*S2bar-2*S1bar.*(S2bar.^3)/sum(3*(S1bar.^2).*(S2bar.^2)-0.5*(S1bar.^4)-0.5*(S2bar.^4)))

V=[cos(Phi0) sin(Phi0); -sin(Phi0) cos(Phi0)];
S1bar=V(1,1)*S1bar+V(1,2)*S2bar;
S2bar=V(2,1)*S1bar+V(2,2)*S2bar;

%rescaling of vectors
min1=min(min(min(S1bar)));
S1bar=S1bar-min1;
max1=max(max(max(S1bar)));
S1bar=S1bar*(255/max1);

min2=min(min(min(S2bar)));
S2bar=S2bar-min2;
max2=max(max(max(S2bar)));
S2bar=S2bar*(255/max2);

subplot(3,2,5),plot(396.600000000000:4.700000000000:995.500000000000,S1bar);
subplot(3,2,6),plot(396.600000000000:4.700000000000:995.500000000000,S2bar);