root / org.gvsig.toolbox / trunk / org.gvsig.toolbox / org.gvsig.toolbox.algorithm / src / main / java / es / unex / sextante / statisticalMethods / pca / PCAAlgorithm.java @ 59
History | View | Annotate | Download (6.97 KB)
1 | 59 | nbrodin | |
---|---|---|---|
2 | |||
3 | package es.unex.sextante.statisticalMethods.pca; |
||
4 | |||
5 | |||
6 | import java.util.ArrayList; |
||
7 | import java.util.Arrays; |
||
8 | |||
9 | import Jama.EigenvalueDecomposition; |
||
10 | import Jama.Matrix; |
||
11 | import es.unex.sextante.additionalInfo.AdditionalInfoMultipleInput; |
||
12 | import es.unex.sextante.additionalInfo.AdditionalInfoNumericalValue; |
||
13 | import es.unex.sextante.core.GeoAlgorithm; |
||
14 | import es.unex.sextante.core.Sextante; |
||
15 | import es.unex.sextante.dataObjects.IRasterLayer; |
||
16 | import es.unex.sextante.exceptions.GeoAlgorithmExecutionException; |
||
17 | import es.unex.sextante.exceptions.RepeatedParameterNameException; |
||
18 | import es.unex.sextante.outputs.OutputRasterLayer; |
||
19 | import es.unex.sextante.parameters.RasterLayerAndBand; |
||
20 | |||
21 | |||
22 | public class PCAAlgorithm |
||
23 | extends
|
||
24 | GeoAlgorithm { |
||
25 | |||
26 | public static final String RESULT = "RESULT"; |
||
27 | public static final String NBANDS = "NBANDS"; |
||
28 | public static final String INPUT = "INPUT"; |
||
29 | |||
30 | private static double NODATA = -9999999.; |
||
31 | |||
32 | private ArrayList m_RasterLayers; |
||
33 | private IRasterLayer[] m_Windows; |
||
34 | private int m_iBands[]; |
||
35 | private int m_iNX, m_iNY; |
||
36 | private double m_dMean[]; |
||
37 | private int m_iNumBands; |
||
38 | |||
39 | |||
40 | @Override
|
||
41 | public void defineCharacteristics() { |
||
42 | |||
43 | setUserCanDefineAnalysisExtent(true);
|
||
44 | setGroup(Sextante.getText("Statistical_methods"));
|
||
45 | setName(Sextante.getText("Principal_Components_Analysis"));
|
||
46 | |||
47 | try {
|
||
48 | m_Parameters.addMultipleInput(INPUT, Sextante.getText("Input_bands"), AdditionalInfoMultipleInput.DATA_TYPE_BAND, true); |
||
49 | m_Parameters.addNumericalValue(NBANDS, Sextante.getText("Bands_number"),
|
||
50 | AdditionalInfoNumericalValue.NUMERICAL_VALUE_INTEGER, 1, 1, Integer.MAX_VALUE); |
||
51 | addOutputRasterLayer(RESULT, Sextante.getText("Principal_Components_Analysis"),
|
||
52 | OutputRasterLayer.NUMBER_OF_BANDS_UNDEFINED); |
||
53 | } |
||
54 | catch (final RepeatedParameterNameException e) { |
||
55 | Sextante.addErrorToLog(e); |
||
56 | } |
||
57 | |||
58 | } |
||
59 | |||
60 | |||
61 | @Override
|
||
62 | public boolean processAlgorithm() throws GeoAlgorithmExecutionException { |
||
63 | |||
64 | int i;
|
||
65 | int x, y;
|
||
66 | int iBand;
|
||
67 | int iInputBand;
|
||
68 | double dValue;
|
||
69 | double dPCValue;
|
||
70 | |||
71 | m_RasterLayers = m_Parameters.getParameterValueAsArrayList(INPUT); |
||
72 | m_iNumBands = m_Parameters.getParameterValueAsInt(NBANDS); |
||
73 | |||
74 | if ((m_RasterLayers.size() == 0) || (m_iNumBands > m_RasterLayers.size())) { |
||
75 | throw new GeoAlgorithmExecutionException(Sextante.getText("Error_not_enough_bands_selected")); |
||
76 | } |
||
77 | |||
78 | m_iNX = getAnalysisExtent().getNX(); |
||
79 | m_iNY = getAnalysisExtent().getNY(); |
||
80 | |||
81 | final double coVar[][] = getCovarMatrix(); |
||
82 | final Matrix coVarMatrix = new Matrix(coVar); |
||
83 | final EigenvalueDecomposition eigenvalueDecomp = new EigenvalueDecomposition(coVarMatrix); |
||
84 | final Matrix eigenvectors = eigenvalueDecomp.getV();
|
||
85 | final double[] eigenvalues = eigenvalueDecomp.getRealEigenvalues(); |
||
86 | final ValueAndOrder[] vao = new ValueAndOrder[eigenvalues.length]; |
||
87 | for (i = 0; i < eigenvalues.length; i++) { |
||
88 | vao[i] = new ValueAndOrder(eigenvalues[i], i);
|
||
89 | } |
||
90 | |||
91 | Arrays.sort(vao);
|
||
92 | final IRasterLayer pca = getNewRasterLayer(RESULT, Sextante.getText("Principal_Components_Analysis"), |
||
93 | IRasterLayer.RASTER_DATA_TYPE_DOUBLE, m_iNumBands); |
||
94 | |||
95 | for (y = 0; y < m_iNY; y++) { |
||
96 | for (x = 0; x < m_iNX; x++) { |
||
97 | for (iBand = 0; (iBand < m_iNumBands) && setProgress(iBand, m_iNumBands); iBand++) { |
||
98 | dPCValue = 0;
|
||
99 | for (iInputBand = 0; iInputBand < m_RasterLayers.size(); iInputBand++) { |
||
100 | dValue = m_Windows[iInputBand].getCellValueAsDouble(x, y, m_iBands[iInputBand]); |
||
101 | dPCValue += (eigenvectors.get(iInputBand, vao[iBand].getOrder()) * dValue); |
||
102 | } |
||
103 | pca.setCellValue(x, y, iBand, dPCValue); |
||
104 | } |
||
105 | } |
||
106 | } |
||
107 | |||
108 | |||
109 | return !m_Task.isCanceled();
|
||
110 | |||
111 | } |
||
112 | |||
113 | |||
114 | private double[][] getCovarMatrix() { |
||
115 | |||
116 | int i, j;
|
||
117 | |||
118 | final double dCovar[][] = new double[m_RasterLayers.size()][m_RasterLayers.size()]; |
||
119 | m_dMean = new double[m_RasterLayers.size()]; |
||
120 | |||
121 | m_Windows = new IRasterLayer[m_RasterLayers.size()];
|
||
122 | m_iBands = new int[m_RasterLayers.size()]; |
||
123 | |||
124 | for (i = 0; i < m_RasterLayers.size(); i++) { |
||
125 | final RasterLayerAndBand rab = (RasterLayerAndBand) m_RasterLayers.get(i);
|
||
126 | m_Windows[i] = rab.getRasterLayer(); |
||
127 | m_Windows[i].setWindowExtent(this.getAnalysisExtent());
|
||
128 | m_Windows[i].setInterpolationMethod(IRasterLayer.INTERPOLATION_BSpline); |
||
129 | m_dMean[i] = m_Windows[i].getMeanValue(); |
||
130 | m_iBands[i] = rab.getBand(); |
||
131 | } |
||
132 | |||
133 | |||
134 | final int iTotal = (int) (m_RasterLayers.size() * m_RasterLayers.size() / 2.); |
||
135 | int iCount = 0; |
||
136 | for (i = 0; (i < m_RasterLayers.size() - 1) && setProgress(iCount, iTotal); i++) { |
||
137 | dCovar[i][i] = 1.0;
|
||
138 | iCount++; |
||
139 | for (j = i + 1; j < m_RasterLayers.size(); j++) { |
||
140 | dCovar[i][j] = dCovar[j][i] = getCovar(i, j); |
||
141 | iCount++; |
||
142 | } |
||
143 | } |
||
144 | |||
145 | |||
146 | return dCovar;
|
||
147 | |||
148 | } |
||
149 | |||
150 | |||
151 | private double getCovar(final int i, |
||
152 | final int j) { |
||
153 | |||
154 | int x, y;
|
||
155 | int iValues = 0; |
||
156 | double dValuei, dValuej;
|
||
157 | double dSum = 0; |
||
158 | |||
159 | for (y = 0; y < m_iNY; y++) { |
||
160 | for (x = 0; x < m_iNX; x++) { |
||
161 | dValuei = m_Windows[i].getCellValueAsDouble(x, y, m_iBands[i]); |
||
162 | dValuej = m_Windows[j].getCellValueAsDouble(x, y, m_iBands[j]); |
||
163 | if (!m_Windows[i].isNoDataValue(dValuei) && !m_Windows[j].isNoDataValue(dValuej)) {
|
||
164 | dSum += (dValuei - m_dMean[i]) * (dValuej - m_dMean[j]); |
||
165 | iValues++; |
||
166 | } |
||
167 | } |
||
168 | if (iValues > 1) { |
||
169 | return dSum / (iValues - 1); |
||
170 | } |
||
171 | else {
|
||
172 | return NODATA;
|
||
173 | } |
||
174 | } |
||
175 | |||
176 | return NODATA;
|
||
177 | |||
178 | } |
||
179 | |||
180 | private class ValueAndOrder |
||
181 | implements
|
||
182 | Comparable {
|
||
183 | |||
184 | private final int m_iOrder; |
||
185 | private final double m_dValue; |
||
186 | |||
187 | |||
188 | public ValueAndOrder(final double dValue, |
||
189 | final int iOrder) { |
||
190 | |||
191 | m_dValue = dValue; |
||
192 | m_iOrder = iOrder; |
||
193 | |||
194 | } |
||
195 | |||
196 | |||
197 | public int getOrder() { |
||
198 | |||
199 | return m_iOrder;
|
||
200 | |||
201 | } |
||
202 | |||
203 | |||
204 | public double getValue() { |
||
205 | |||
206 | return m_dValue;
|
||
207 | |||
208 | } |
||
209 | |||
210 | |||
211 | public int compareTo(final Object vao) throws ClassCastException { |
||
212 | |||
213 | if (!(vao instanceof ValueAndOrder)) { |
||
214 | throw new ClassCastException(); |
||
215 | } |
||
216 | |||
217 | final double dValue = ((ValueAndOrder) vao).getValue(); |
||
218 | final double dDif = this.m_dValue - dValue; |
||
219 | |||
220 | if (dDif > 0.0) { |
||
221 | return 1; |
||
222 | } |
||
223 | else if (dDif < 0.0) { |
||
224 | return -1; |
||
225 | } |
||
226 | else {
|
||
227 | return 0; |
||
228 | } |
||
229 | |||
230 | } |
||
231 | } |
||
232 | } |