Great news. Thanks. Michael
On Jun 25, 2012, at 12:51 PM, Markus Metz wrote: > On Mon, Jun 25, 2012 at 8:21 PM, Michael Barton <[email protected]> > wrote: >> On Jun 25, 2012, at 12:03 PM, Markus Metz wrote: >>> >>> Here is an example based on on the landsat imagery in mapset landsat, >>> North Carolina sample dataset: >>> >>> i.pca without rescaling of the output: >>> i.pca input=lsat7_2000_10,lsat7_2000_20,lsat7_2000_30 >>> output_prefix=lsat7_2000_pc rescale=0,0 >>> >>> The eigenvectors are >>> 0.4494, 0.5128, 0.7315 >>> 0.6726, 0.3447,-0.6548 >>> 0.5879,-0.7863, 0.1901 >>> >>> this is a square 3x3 matrix. R gives the inverse matrix as >>> >>> 0.4493372, 0.6726549, 0.5879235 >>> 0.5128129, 0.3446144,-0.7862660 >>> 0.7315070,-0.6548316, 0.1899992 >>> > [...] >>> >>> Without using R, the inverse of a 3x3 matrix can be manually >>> calculated as follows: >>> >>> original matrix: >>> a b c >>> d e f >>> g h k >>> >>> inverse matrix: >>> (ek - fh) (ch - bk) (bf - ce) >>> (fg - dk) (ak - cg) (cd - af) >>> (dh - eg) (gb - ah) (ae - bd) >> >> Your example from R above *looks* like the inverse of matrix... >> >> a b c >> d e f >> g h k >> >> is (given rounding errors) ... >> >> a d g >> b e h >> c f k >> >> That is, transpose rows and columns. >> Is this just a coincidence or a reasonable approximation? > > Oops, right. A PCA uses an orthogonal transformation, in which case > the inverse of the transformation matrix is identical to the transpose > of the transformation matrix. Easy solution. >> >>> >>> each entry in this inverse matrix needs to be multiplied with >>> >>> 1 / (a(ek - fh) - b(kd - fg) + c(dh - eg)) >> >> by a(ek - fh)... , you mean a = factor score 1 of a cell, b = factor score 2 >> of a cell, c = factor score 3 of a cell (not eigenvectors a,b,c or the >> matrix), right? > > No, a,b,c are here Eigenvector values. But this formula is the general > case for the inverse of a 3x3 matrix and not needed here, because the > inverse and the transpose are the same. > > Markus M > >> >>> >>> in order to get the transformation matrix to be applied to the pc scores. >>> >>> The formula is more complex for matrices with more than 3 dimensions, >>> e.g. for i.pca with 6 input bands. >>> >>>> >>>> So my questions are: >>>> >>>> 1) Do I have the inversion correct in terms of how it needs to be >>>> calculated in GRASS? >>> >>> I don't think so. Using your formula and testing for the difference gives me >>> min=-546.255208609669 >>> max=174.771853180844 >>> range=721.027061790513 >>> mean=79.9249815357317 >>> >>> a bit too much IMHO. >> >> Right. Looks like I'll need to use r.univar to calculate the mean of each >> band to add back in. Thanks again. I will try this out manually and test >> before updating the pan sharpening routine. >> >>> >>>> 2) Even though the mean of all bands seems to be subtracted from each >>>> factor score, does inverting the matrix mean that the mean of *each* band >>>> is added back to its transformed value? Adding either the mean of all >>>> original bands or 1 original band seems to produce values that are much >>>> higher than the original, and so need to be rescaled. Maybe this is OK. >>> >>> The mean of each band is added back to each recalculated band. See >>> above r.mapcalc formula for band 10. >>> >>>> 3) I do not normalize or rescale in i.pca. This seems to make it easier to >>>> do the inverse PCA with fewer calculations. Is there any reason I *should* >>>> rescale and/or normalize? >>> >>> Normalization applies to the input of i.pca and is by default not >>> done. It needs to be done for heterogenous input data such as e.g. >>> rainfall, temperature, NDVI, etc. Rescaling is automatically applied >>> to the output of i.pca unless explicitly disabled with rescale=0,0. >> >> Thanks. That is what I thought. >> >> Michael >> >>> >>> Markus M >>> >>>> >>>> Thanks much >>>> Michael >>>> >>>> >>>> On Jun 25, 2012, at 10:11 AM, Markus Metz wrote: >>>> >>>>> An inverse PCA can be regarded as the inverse of a transformation >>>>> using matrix notation. PC scores are calculated with >>>>> >>>>> b = A a >>>>> >>>>> with A being the transformation matrix composed of the Eigenvectors, a >>>>> being the vector of the original values and b the PC scores. What you >>>>> now need is inverse of A, A^-1. The original values can then be >>>>> retrieved with >>>>> >>>>> A^-1 b = a >>>>> >>>>> A^-1 is the inverse of the transformation matrix A which you can get >>>>> in R with solve(A). >>>>> >>>>> For a PCA, the original values are usually shifted to the mean and >>>>> optionally scaled to stddev before computing the Eigenvectors. The >>>>> mean shift is always performed by i.pca, scaling is optional. That >>>>> means that A^-1 b gives the original values shifted to the mean and >>>>> maybe scaled, and the mean of each original band needs to be added to >>>>> get the original values used as input to i.pca. With scaling applied, >>>>> the shifted values need to be multiplied by the stddev for each >>>>> original band. >>>>> >>>>> HTH, >>>>> >>>>> Markus M >>>>> >>>>> >>>>> On Tue, Jun 19, 2012 at 12:46 AM, Michael Barton <[email protected]> >>>>> wrote: >>>>>> The constant (i.e., the band mean) was in the pca primer that someone >>>>>> sent me a link to in this discussion. Using the Eigenvectors resulting >>>>>> from i.pca, I can achieve the results of i.pca using my formula below. >>>>>> This is essentially the same as your formula minus the constant--which >>>>>> doesn't really make much (of any) difference in the final result. >>>>>> >>>>>> However, my question is about performing an *inverse pca*--getting back >>>>>> to the original values from the principal components maps. The idea of >>>>>> pca sharpening is that you perform a pca, then do an inverse pca >>>>>> substituting the pan band for pc1 to enhance the resolution. The >>>>>> equations I show below seem to work, but what I've read indicates that >>>>>> it is not possible to *exactly* get the original values back; you can >>>>>> only approximate them. >>>>>> >>>>>> Michael >>>>>> >>>>>> >>>>>> On Jun 17, 2012, at 10:48 AM, Duccio Rocchini wrote: >>>>>> >>>>>>> Dear all, >>>>>>> first, sorry for the delay... >>>>>>> Here are my 2 cents to be added to the discussion. I re-took in my >>>>>>> hands the John Jensen book. >>>>>>> Accordingly >>>>>>> >>>>>>> new brightness values1,1,1 = a1,1*BV1,1,1 +a2,1*BV1,1,2..... + >>>>>>> an,1*BV1,1,m >>>>>>> >>>>>>> where a=eigenvector and BV=original brightness value. >>>>>>> >>>>>>> I found no evidence for the mean term: "- ((b1+b2+b3)/3)" >>>>>>> >>>>>>> Michael: do you have a proof/reference for that? >>>>>>> >>>>>>> P.S. thanks for involving me in this discussion which is really >>>>>>> stimulating! >>>>>>> >>>>>>> Duccio >>>>>>> >>>>>>> 2012/6/7 Michael Barton <[email protected]>: >>>>>>>> >>>>>>>> I think I've figured it out. >>>>>>>> >>>>>>>> If (ev1-1, ev1-2, ev1-3) are the eigenvectors of the first principal >>>>>>>> component for 3 imagery bands (b1, b2, b3), the corresponding factor >>>>>>>> scores of the PC1, PC2, and PC3 maps (fs1, fs2, fs3) are calculated as: >>>>>>>> >>>>>>>> fs1 = (ev1-1*b1) + (ev1-2*b2) + (ev1-3*b3) - ((b1+b2+b3)/3) >>>>>>>> fs2 = (ev2-1*b1) + (ev2-2*b2) + (ev2-3*b3) - ((b1+b2+b3)/3) >>>>>>>> fs3 = (ev3-1*b1) + (ev3-2*b2) + (ev3-3*b3) - ((b1+b2+b3)/3) >>>>>>>> >>>>>>>> So to do an inverse PCA on the same data you need to do the following: >>>>>>>> >>>>>>>> b1' = (fs1/ev1-1) + (fs2/ev2-1) + (fs3/ev3-1) >>>>>>>> b2' = (fs1/ev1-2) + (fs2/ev2-2) + (fs3/ev3-2) >>>>>>>> b3' = (fs1/ev1-3) + (fs2/ev2-3) + (fs3/ev3-3) >>>>>>>> >>>>>>>> Adding the constant back on doesn't really seem to matter because you >>>>>>>> need to rescale b1' to b1, b2' to b2, and b3' to b3 anyway. >>>>>>>> >>>>>>>> Michael >>>>>>>> >>>>>>>> On Jun 7, 2012, at 1:55 AM, Markus Neteler wrote: >>>>>>>> >>>>>>>>> Hi Duccio, >>>>>>>>> >>>>>>>>> On Wed, Jun 6, 2012 at 11:39 PM, Michael Barton >>>>>>>>> <[email protected]> wrote: >>>>>>>>>> On Jun 6, 2012, at 2:20 PM, Markus Neteler wrote: >>>>>>>>>>> On Wed, Jun 6, 2012 at 5:09 PM, Michael Barton >>>>>>>>>>> <[email protected]> wrote: >>>>>>>>> ... >>>>>>>>>>>> I'm working on a pan sharpening script that will combine your >>>>>>>>>>>> i.fusion.brovey with options to do other pan sharpening methods. >>>>>>>>>>>> An IHS transformation is easy. I think that a PCA sharpening is >>>>>>>>>>>> doable too if I can find an equation to rotate the PCA channels >>>>>>>>>>>> back into unrotated space--since i.pca does provide the >>>>>>>>>>>> eigenvectors. >>>>>>>>>>> >>>>>>>>>>> Maybe there is material in (see m.eigenvector) >>>>>>>>>>> http://grass.osgeo.org/wiki/Principal_Components_Analysis >>>>>>>>>> >>>>>>>>>> This has a lot of good information and ALMOST has what I need. >>>>>>>>>> Everything I read suggests that there is a straightforward way to >>>>>>>>>> get the original values from the factor scores if you have the >>>>>>>>>> eigenvectors. But no one I've yet found provides the equation or >>>>>>>>>> algorithm to do it. >>>>>>>>> >>>>>>>>> @Duccio: any idea about this by chance? >>>>>>>>> >>>>>>>>> thanks >>>>>>>>> Markus >>>>>>>> >>>>>>>> _____________________ >>>>>>>> C. Michael Barton >>>>>>>> Visiting Scientist, Integrated Science Program >>>>>>>> National Center for Atmospheric Research & >>>>>>>> University Corporation for Atmospheric Research >>>>>>>> 303-497-2889 (voice) >>>>>>>> >>>>>>>> Director, Center for Social Dynamics & Complexity >>>>>>>> Professor of Anthropology, School of Human Evolution & Social Change >>>>>>>> Arizona State University >>>>>>>> www: http://www.public.asu.edu/~cmbarton, http://csdc.asu.edu >>>>>>>> >>>>>>> >>>>>>> >>>>>>> >>>>>>> -- >>>>>>> Duccio Rocchini, PhD >>>>>>> >>>>>>> http://gis.cri.fmach.it/rocchini/ >>>>>>> >>>>>>> Fondazione Edmund Mach >>>>>>> Research and Innovation Centre >>>>>>> Department of Biodiversity and Molecular Ecology >>>>>>> GIS and Remote Sensing Unit >>>>>>> Via Mach 1, 38010 San Michele all'Adige (TN) - Italy >>>>>>> Phone +39 0461 615 570 >>>>>>> [email protected] >>>>>>> [email protected] >>>>>>> skype: duccio.rocchini >>>>>> >>>>>> _____________________ >>>>>> C. Michael Barton >>>>>> Visiting Scientist, Integrated Science Program >>>>>> National Center for Atmospheric Research & >>>>>> University Consortium for Atmospheric Research >>>>>> 303-497-2889 (voice) >>>>>> >>>>>> Director, Center for Social Dynamics & Complexity >>>>>> Professor of Anthropology, School of Human Evolution & Social Change >>>>>> Arizona State University >>>>>> www: http://www.public.asu.edu/~cmbarton, http://csdc.asu.edu >>>>>> >>>>>> >>>>>> >>>>>> >>>>>> >>>>>> _______________________________________________ >>>>>> grass-dev mailing list >>>>>> [email protected] >>>>>> http://lists.osgeo.org/mailman/listinfo/grass-dev >>>> >>>> _____________________ >>>> C. Michael Barton >>>> Visiting Scientist, Integrated Science Program >>>> National Center for Atmospheric Research & >>>> University Corporation for Atmospheric Research >>>> 303-497-2889 (voice) >>>> >>>> Director, Center for Social Dynamics & Complexity >>>> Professor of Anthropology, School of Human Evolution & Social Change >>>> Arizona State University >>>> www: http://www.public.asu.edu/~cmbarton, http://csdc.asu.edu >>>> >> >> _____________________ >> C. Michael Barton >> Visiting Scientist, Integrated Science Program >> National Center for Atmospheric Research & >> University Corporation for Atmospheric Research >> 303-497-2889 (voice) >> >> Director, Center for Social Dynamics & Complexity >> Professor of Anthropology, School of Human Evolution & Social Change >> Arizona State University >> www: http://www.public.asu.edu/~cmbarton, http://csdc.asu.edu >> _____________________ C. Michael Barton Visiting Scientist, Integrated Science Program National Center for Atmospheric Research & University Corporation for Atmospheric Research 303-497-2889 (voice) Director, Center for Social Dynamics & Complexity Professor of Anthropology, School of Human Evolution & Social Change Arizona State University www: http://www.public.asu.edu/~cmbarton, http://csdc.asu.edu _______________________________________________ grass-dev mailing list [email protected] http://lists.osgeo.org/mailman/listinfo/grass-dev
