X
3316

How to use the correlation matrix and specify a subset in PC_ROTATE

This Help Articles provides a code example of how to perform forward and inverse Principal Components Analysis (PCA) rotations using the correlation matrix in ENVI's PC_ROTATE procedure. The code also includes how you can remove unwanted PC bands for the inverse rotation.
 
ENVI's forward and inverse PCA rotations can be performed using the covariance or correlation matricies. This is very straight-forward in the ENVI user interface, but not apparent when using ENVI's PC_ROTATE procedure, which requires the user to calculate statistics for input to the routine. The below code example shows how the correlation matrix can be calculated and used for the forward and inverse PCA rotations and how to subset out unwanted PC bands when performing the inverse PCA.
---------------------------------------------
PRO example_subset_pc_rotate

;select and open a file
file = 'bhtmref.img'
envi_open_file, file, r_fid=fid

;query the file for relevant information
ENVI_FILE_QUERY, fid, dims=dims, NS = ns, NL = nl, NB = nb
pos  = LINDGEN(nb)
out_name = 'test_forward_pca.img'

;calculate statistics for the forward PCA rotation
ENVI_DOIT, 'ENVI_STATS_DOIT', $
   FID = fid, POS = pos, DIMS = dims, cov=cov, $
   MEAN = avg, COMP_FLAG = 5

;convert covariance to correlation
cor = cov
    for i=0L, n_elements(pos)-1 do $
      for j=0L, n_elements(pos)-1 do begin
        div = sqrt(cov[i,i]) * sqrt(cov[j,j])
        if (div gt 0) then cor[j,i] = (cor[i,j] = cov[i,j] / div) $
        else cor[j,i] = (cor[i,j] = 0.)
      endfor
    ; compute evec and eval based on this correlation matrix
    evec = cor
    trired, evec, eval, temp, /double
    triql, eval, temp, evec, /double
    ord = reverse(sort(eval))
    eval = eval[ord]
    evec = evec[*,ord]

;perform the forward PCA rotation
ENVI_DOIT, 'PC_ROTATE', $
   FID = fid, POS = pos, DIMS = dims, $
   MEAN = avg, EVAL = eval, EVEC = evec, $
   OUT_DT = 4, OUT_NAME = out_name, $
   OUT_NB = 6, R_FID = pca_fid, $
   /FORWARD

;again calculate statistics for the inverse PCA
;to get covariance, correlation, evec, and eval
ENVI_DOIT, 'ENVI_STATS_DOIT', $
   FID = fid,  pos=pos, DIMS = dims, cov=cov, $
   MEAN = avg, COMP_FLAG = 5

cor = cov
    for i=0L, n_elements(pos)-1 do $
      for j=0L, n_elements(pos)-1 do begin
        div = sqrt(cov[i,i]) * sqrt(cov[j,j])
        if (div gt 0) then cor[j,i] = (cor[i,j] = cov[i,j] / div) $
        else cor[j,i] = (cor[i,j] = 0.)
      endfor    
    evec = cor
    trired, evec, eval, temp, /double
    triql, eval, temp, evec, /double
    ord = reverse(sort(eval))
    eval = eval[ord]
    evec = evec[*,ord]

;perform the inverse PCA rotation
;For the inverse, set pos = 0 for only the first PCA band
;set a variable called 'limit_nb' to 6 
;This sets the output limit to 6 bands
out_name2 = 'test_inverse_pca.img'
ENVI_DOIT, 'PC_ROTATE', $
   FID = pca_fid, POS = 0, DIMS = dims, $
   MEAN = avg, EVAL = eval, EVEC = evec, $
   OUT_DT = 4, OUT_NAME = out_name2, $
   R_FID = r_fid, limit_nb=6

 END

Review on 12/31/2013 MM