1+ # [E, IdxOfE, Xpca] = VCA(X, M, r, v)
2+ #
3+ ###############################################################################
4+ #
5+ # A FUNCTION TO CALCULATE ENDMEMBERS USING Vertex Component Analysis (VCA)
6+ # REFERENCE:
7+ # José M. P. Nascimento and José M. B. Dias
8+ # "Vertex Component Analysis: A Fast Algorithm to Unmix Hyperspectral Data"
9+ # IEEE Trans. Geosci. Remote Sensing,
10+ # April, 2005
11+ ###############################################################################
12+ ###
13+ ### INPUTS:
14+ ### X: DATA MATRIX B WAVELENGTHS x N PIXELS
15+ ### M: NUMBER OF ENDMEMBERS TO ESTIMATE
16+ ### OUTPUTS:
17+ ### E: B x M MATRIX OF ESTIMATED ENDMEMBERS
18+ ### IdxOfEinX: INDICES OF ENDMEMBERS IN DATA X.
19+ ### (THE ENDMEMBERS ARE SELECTED FROM X)
20+ ### XM: X PROJECTED ONTO FIRST M COMPONENTS OF PCA
21+ #
22+ ### OPTIONAL INPUTS:
23+ ### 'SNR' r: ESTIMATED SIGNAL TO NOISE RATIO IN dB
24+ ### 'verbose' v: logical TOGGLE TO TURN DISPLAYS ON AND OFF
25+ ###############################################################################
26+ ###
27+ ### Authors: José Nascimento ([email protected] ) 28+ ### José Bioucas Dias ([email protected] ) 29+ ### Copyright (c)
30+ ### version: 2.1 (7-May-2004)
31+ ###############################################################################
32+ ###
33+ ### ORIGINALLY MODIFED BY Darth Gader OCTOBER 2018
34+ ### TRANSLATED TO PYTHON BY Ronald Fick November 2018
35+ ###############################################################################
36+
37+ import numpy as np
38+ import numpy .matlib
39+ from scipy .sparse .linalg import svds
40+ import math
41+ import matplotlib .pyplot as plt
42+
43+ def VCA (X , M = 2 , r = - 1 , verbose = True ):
44+ ##
45+ #############################################
46+ # Initializations
47+ #############################################
48+
49+ if (X .size == 0 ):
50+ raise ValueError ('There is no data' )
51+ else :
52+ B , N = X .shape
53+
54+ if (M < 0 or M > B or M != int (M )):
55+ raise ValueError ('ENDMEMBER parameter must be integer between 1 and B' )
56+
57+ ##
58+ ####################
59+ ### ESTIMATE SNR ###
60+ ####################
61+ if (r == - 1 ):
62+ ############################################################
63+ ### THE USER DID NOT INPUT AN SNR SO SNR IS CALCULATED HERE
64+ ############################################################
65+
66+ #############################
67+ ### CALCULATE PCA AND PCT ###
68+ MuX = np .mean (X , axis = 1 ) #axis 1 is the second dimension
69+ MuX .shape = (MuX .size ,1 )
70+ Xz = X - np .matlib .repmat (MuX , 1 , N )
71+ SigmaX = np .cov (X )
72+ U ,S ,V = np .linalg .svd (SigmaX )
73+ Xpca = np .matmul (np .transpose (U [:,0 :M ]), Xz )
74+ ProjComputed = True
75+
76+ ### ESTIMATE SIGNAL TO NOISE ###
77+ SNR = EstSNR (X ,MuX ,Xpca )
78+
79+ ### PRINT SNR ###
80+ if verbose :
81+ print ('Estimated SNR = %g[dB]' % SNR )
82+ else :
83+ ############################################################
84+ ### THE USER DID INPUT AN SNR SO NO SNR CALCUATION NEEDED
85+ ############################################################
86+
87+ ProjComputed = False
88+ if verbose :
89+ print ('Input SNR = %g[dB]' % SNR )
90+
91+ ### SET THRESHOLD TO DETERMINE IF NOISE LEVEL IS LOW OR HIGH ###
92+ SNRThresh = 15 + 10 * math .log10 (M )
93+
94+ if verbose :
95+ print ('SNRThresh= %f SNR= %f Difference= %f' % (SNR ,SNRThresh ,SNRThresh - SNR ))
96+
97+ ###########################
98+ ### END OF ESTIMATE SNR ###
99+ ###########################
100+ ##
101+ ##################################################################
102+ ### PROJECTION. ###
103+ ### SELECT AND CALCULATE PROJECTION ONTO M DIMS ###
104+ ### ###
105+ ### IF SNR IS LOW,IT IS ASSUMED THAT THERE IS STILL NOISE IN ###
106+ ### THE SIGNAL SO REDUCE DIM A BIT MORE. ###
107+ ### ADD A CONSTANT VECTOR TO KEEP IN DIM M. ###
108+ ### ###
109+ ### IF SNR IS HIGH, PROJECT TO SIMPLEX IN DIM M-1 TO REDUCE ###
110+ ### EFFECTS OF VARIABLE ILLUMINATION ###
111+ ##################################################################
112+
113+ if (SNR < SNRThresh ):
114+ ##########################
115+ ### BEGIN LOW SNR CASE ###
116+ ##########################
117+
118+ ### PRINT MESSAGE ###
119+ if verbose :
120+ print ('Low SNR so Project onto Dimension M-1.' )
121+
122+ ### REDUCE SIZE OF PCT MATRIX TO M-1 ###
123+ Dim = M - 1
124+ MuX = np .mean (X , axis = 1 )
125+ MuX .shape = (MuX .size , 1 )
126+ BigMuX = np .matlib .repmat (MuX , 1 , N )
127+ if ProjComputed :
128+ U = U [:,0 :Dim ]
129+ else :
130+ Xz = X - BigMuX
131+ SigmaX = np .cov (np .transpose (X ))
132+ U ,S ,V = np .linalg .svd (SigmaX )
133+ #U,S,V = np.linalg.svd(SigmaX, full_matrices=True)
134+ Xpca = np .matmul (np .transpose (U ), Xz )
135+
136+ ### REDUCE DIMENSIONALITY IN PCA DOMAIN ###
137+ XpcaReduced = Xpca [0 :Dim ,:]
138+
139+ ### RECONSTRUCT X "WITHOUT NOISE" BY PROJECTING BACK TO ORIGINAL SPACE ###
140+ XNoiseFree = np .matmul (U , XpcaReduced ) + BigMuX
141+
142+ ### CONCATENATE CONSTANT VEC = MAX NORM OF ALL DATA POINTS ###
143+ BiggestNorm = np .sqrt (max (np .sum (np .square (XpcaReduced ),1 )))
144+
145+ YpcaReduced = np .concatenate ([XpcaReduced , BiggestNorm * np .ones ((1 ,N ))])
146+ ########################
147+ ### END LOW SNR CASE ###
148+ ########################
149+ else :
150+ ###########################
151+ ### BEGIN HIGH SNR CASE ###
152+ ###########################
153+ if verbose :
154+ print ('High SNR so project onto dimension M' )
155+
156+ ### CONTRUCT "PCA-LIKE" MATRIX BY DIAGONALIZING CORRELATION MATRIX ###
157+ ### IF SUBTRACT THE MEAN, THEN MUPCA WILL BE 0, SO NO MEAN SUBTRACTION ###
158+ # xb = a: solve b.T x.T = a.T instead
159+ U ,S ,V = np .linalg .svd (np .matmul (X ,np .transpose (X ))/ N )
160+
161+ ### CALC PCT WITHOUT SUBTRACTING MEAN AND THEN RECONSTRUCT ###
162+ XpcaReduced = np .matmul (U .T ,X )
163+
164+ ### RECONSTRUCT X VIA INVERSE XFORM ON REDUCED DIM PCA DATA ###
165+ ### XXX PDG NOT SURE IF THIS IS CORRECT IN ORIGINAL CODE OR HERE. ###
166+ ### I THINK THE MEAN SHOULD BE ADDED LIKE IN LOW SNR CASE ###
167+ XNoiseFree = np .matmul (U , XpcaReduced [0 :B ,:]) # again in dimension L (note that x_p has no null mean)
168+
169+ #################################################################
170+ ### CALCULATE NORMALIZED PROJECTION ###
171+ ### SEE LAST PARAGRAPH OF SECTION A AND FIGURE 4 IN VCA PAPER.###
172+ ### MEAN OF THE PROJECT DATA IS USED AS u FROM THAT SECTION ###
173+ Mupca = np .mean (XpcaReduced ,1 )
174+ Mupca .shape = (Mupca .size , 1 )
175+ Denom = np .sum (np .multiply (XpcaReduced , np .tile (Mupca ,[1 , N ])))
176+ YpcaReduced = np .divide (XpcaReduced , np .tile (Denom , [B , 1 ]))
177+
178+ ##
179+ ###########################################################
180+ # VCA ALGORITHM
181+ ###########################################################
182+ ###
183+ ### INITIALIZE ARRAY OF INDICES OF ENDMEMBERS, IdxOfE
184+ ### INITIALIZE MATRIX OF ENDMEMBERS IN PCA SPACE, Epca
185+ ### DO M TIMES
186+ ### PICK A RANDOM VECTOR w
187+ ### USE w TO FIND f IN NULL SPACE OF Epca
188+ ### PROJECT f ONTO ALL PCA DATA PTS
189+ ### FIND INDEX OF PCA DATA PT WITH MAX PROJECTION
190+ ### USE INDEX TO ADD PCA DATA PT TO Epca
191+ ### END DO
192+ ###
193+ ### USE INDICES TO SELECT ENDMEMBERS IN SPECTRAL SPACE
194+ ###########################################################
195+
196+ VCAsize = YpcaReduced .shape [0 ]
197+
198+ IdxOfE = np .zeros (VCAsize )
199+ IdxOfE = IdxOfE .astype (int )
200+ Epca = np .zeros ((VCAsize ,VCAsize ))
201+ Epca [VCAsize - 1 ,0 ] = 1
202+ for m in range (0 ,VCAsize ):
203+ w = np .random .rand (VCAsize ,1 )
204+ f = w - np .matmul (np .matmul (Epca ,np .linalg .pinv (Epca )), w )
205+ f = f / np .sqrt (np .sum (np .square (f )))
206+ v = np .matmul (f .T ,YpcaReduced )
207+ absV = np .abs (v )
208+ absV = np .squeeze (absV )
209+ IdxOfE [m ] = np .argmax (absV )
210+ v_max = absV [IdxOfE [m ]]
211+ Epca [:,m ] = YpcaReduced [:,IdxOfE [m ]]
212+
213+ E = XNoiseFree [:,IdxOfE ]
214+
215+ return E [:,0 :M ], IdxOfE [0 :M ], Xpca
216+ #############################################
217+ ### END OF VCA FUNCTION
218+ #############################################
219+ #############################################
220+ ##
221+ #############################################
222+ #############################################
223+ ### FUNCTION TO ESTIMATE SNR
224+ #############################################
225+
226+ def EstSNR (X ,MuX ,Xpca ):
227+ #####################################################
228+ ### THE ESTIMATED SNR IS EQUIVALENT TO THE RATIO OF
229+ ### POWER OF SIGNAL TO POWER OF NOISE
230+ ### BUT IS COMPUTED AS
231+ ### POWER OF (SIGNAL + NOISE) TO POWER OF NOISE
232+ #####################################################
233+ B , N = X .shape
234+ M , N = Xpca .shape
235+
236+ ### POWER OF SIGNAL + NOISE, SIGNAL, AND NOISE, RESP. ###
237+ Psn = np .sum (np .square (X .flatten ()))/ N
238+ Ps = np .sum (np .square (Xpca .flatten ()))/ N + np .matmul (np .transpose (MuX ),MuX )
239+ Pn = Psn - Ps
240+ SNR = 10 * np .log10 (Ps / Pn )
241+
242+ ### OLD SNR CODE AS IN VCA PAPER. NOT MUCH DIFFERENT ###
243+ ### BECAUSE (M/B) IS SMALL ###
244+ ### SNR = 10*log10( (Ps - M/B*Psn)/(Psn - Ps) );
245+ ### fprintf('SNR= %8.4f RatSNR= %8.4f\n', SNR, RatSNR);
246+
247+ return SNR
248+ ###############
249+ ### THE END ###
250+ ###############
0 commit comments