How to retro-engineer hyperPCA and obtain the components from the coefficients and original input?

2 次查看（过去 30 天）
Federico
Federico 2022-2-10
编辑： Mahesh Belagal Sarvisetty 2022-2-17
Dear all, I am using the hyperpca matlab function to obtain the principal components for some hyperspectral images. The function is called like this:
[components,coefficients,variance]=hyperpca(hypercube,numcoefficients);
Where hypercube is the hyperspectral image passed as an m*n*p int16 matrix.
If I understand correctly PCA is a function reprojecting the input data (hypercube) into the components space (components) and therefore a linear function.
My expectations would be that multiplying the input data by the coefficients would return the principal components. I proceeded as follows:
reshaped_inputs=reshape(hypercube,size(hypercube,1)*size(hypercube,2),size(hypercube,3));
reshaped_new_components=single(reshaped_inputs)*coefficients; % need to cast the inputs from int16 to single otherwise mtimes won't work
new_components=reshape(reshaped_new_components,size(hypercube,1),size(hypercube,2),size(coefficients,2));
I would then expect that new_components is == components, but instead the difference between those two matrices is far from zero.
Why is that? Is there some sort of scaling/centering that I'm missing?
Thanks for any inputs!

请先登录，再回答此问题。

采纳的回答

Mahesh Belagal Sarvisetty
By default, hyperpca mean-centering the input data in order to compute principal components. So, you can set MeanCentered value as false to match new_components with the components from the hyperpca function or you can mean center the input data before multiplying with the coefficients.
Here is quick snippet with MeanCentered value false
% Read Hypercube
hCube = hypercube("paviaU.dat");
% Perform PCA
[princComp, coeff] = hyperpca(hCube,20,"MeanCentered",false);
% Extract and reshape the DataCube
dataCube = (hCube.DataCube);
sz = size(dataCube);
dataCube = reshape(dataCube,[sz(1)*sz(2) sz(3)]);
% Generate the new principal component
newPrincComp = dataCube * coeff;
sz = size(princComp);
newPrincComp = reshape(newPrincComp,[sz(1) sz(2) sz(3)]);
% Calculate max absolute difference
diff = imabsdiff(princComp,newPrincComp);
max(diff(:))
  2 个评论
Federico
Federico 2022-2-16
Dear @Mahesh Belagal Sarvisetty, thank you for your answer!
I followed your example, applied it to my datacube, but, unfortunately, there's still something off when I attempt it.
The maximum difference I see between the original principal components and the back-calculated ones is quite huge (2.0621e+04).
I don't know what could be the cause, my hypothesis are:
  • Could too many zeroes on the image border affect the PCA computation? The Pavia image had some zero pixels, but not that many...
  • I have to convert my datacube from int16 to single. I've checked some numbers manually (since many functions, such as the simple A-B do not work on int16 data) and I can't find any difference between the int16 and the single version of the datacube.
  • I'm asking something that's not really possible. For example, the same toolbox has the inverseProjection function which should reconstruct the full hyperspectral datacube from the PCA components and coefficients. Even on the Indian Pines example data the minimum difference between the reconstructed cube and the old datacube is >960
% MathWorks Example
hcube=hypercube('indian_pines.dat');
[pcDataCube,coeff]=hyperpca(hcube,10);
reconstructedData=inverseProjection(pcDataCube,coeff);
% Difference calculation
diff=imabsdiff(hcube.DataCube,reconstructedData);
disp(['Minimum difference: ',sprintf('%.2f',min(diff(:)))]);
disp(['Average difference: ',sprintf('%.2f',mean(diff(:)))]);
disp(['Maximum difference: ',sprintf('%.2f',max(diff(:)))]);
% Displayed output:
'Minimum difference: 961.11'
'Average difference: 2503.44'
'Maximum difference: 9315.66'
Any idea or comment?
Mahesh Belagal Sarvisetty
Mahesh Belagal Sarvisetty 2022-2-17
编辑：Mahesh Belagal Sarvisetty 2022-2-17
inverseProjection input 'coeff' contains coefficients for only 10 principal components. In order to reconstruct original data we have use coefficients of all principal components and also set the MeanCentered value as false.
Take a look at this example, reconstructedData is same as original datacube
hcube = hypercube('indian_pines.dat');
[pcDataCube,coeff] = hyperpca(hcube,220,"MeanCentered",false);
reconstructedData = inverseProjection(pcDataCube,coeff);
% Difference calculation
diff = imabsdiff(hcube.DataCube,reconstructedData);
disp(['Minimum difference: ',sprintf('%.2f',min(diff(:)))]);
disp(['Average difference: ',sprintf('%.2f',mean(diff(:)))]);
disp(['Maximum difference: ',sprintf('%.2f',max(diff(:)))]);

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Hyperspectral Image Processing 的更多信息

标签

产品


版本

R2021a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by