svn-gvsig-desktop / trunk / extensions / extGeoreferencing / src / com / iver / cit / gvsig / GeoOperations.java @ 4840
History | View | Annotate | Download (9.67 KB)
1 |
package com.iver.cit.gvsig; |
---|---|
2 |
|
3 |
import java.awt.geom.Point2D; |
4 |
import java.io.BufferedOutputStream; |
5 |
import java.io.DataOutputStream; |
6 |
import java.io.File; |
7 |
import java.io.FileOutputStream; |
8 |
import java.io.IOException; |
9 |
import java.util.zip.DataFormatException; |
10 |
|
11 |
import com.iver.cit.gvsig.fmap.layers.FLyrPoints; |
12 |
import com.iver.cit.gvsig.utils.AffineT; |
13 |
import com.iver.cit.gvsig.utils.MathUtils; |
14 |
import com.iver.cit.gvsig.utils.PointT; |
15 |
|
16 |
/**
|
17 |
* @author Nacho Brodin (brodin_ign@gva.es)
|
18 |
*/
|
19 |
public class GeoOperations{ |
20 |
|
21 |
private FLyrPoints lyrPoints = null; |
22 |
private int order = 1; |
23 |
private PointT[] src = null; |
24 |
private PointT[] dst = null; |
25 |
private AffineT affine = null; |
26 |
|
27 |
/**
|
28 |
* Constructor
|
29 |
*/
|
30 |
public GeoOperations(){}
|
31 |
|
32 |
/**
|
33 |
* Constructor. Crea las estructuras de puntos y calcula la transformaci?n af?n.
|
34 |
* @param lyr
|
35 |
*/
|
36 |
public GeoOperations(FLyrPoints lyr){
|
37 |
this.lyrPoints = lyr;
|
38 |
src = new PointT[lyr.getCountPoints()];
|
39 |
dst = new PointT[lyr.getCountPoints()];
|
40 |
int nPoint = 0; |
41 |
for(int i = 0; i<lyr.getCountPoints(); i++){ |
42 |
if(lyr.getPoint(i).active == true){ |
43 |
src[nPoint] = new PointT();
|
44 |
dst[nPoint] = new PointT();
|
45 |
src[nPoint].setX(lyr.getPoint(i).pixelPoint.getX()); |
46 |
src[nPoint].setY(lyr.getPoint(i).pixelPoint.getY()); |
47 |
src[nPoint].setI(lyr.getCountPoints()); |
48 |
dst[nPoint].setX(lyr.getPoint(i).mapPoint.getX()); |
49 |
dst[nPoint].setY(lyr.getPoint(i).mapPoint.getY()); |
50 |
dst[nPoint].setI(lyr.getCountPoints()); |
51 |
nPoint++; |
52 |
} |
53 |
} |
54 |
|
55 |
if(lyr.getCountPoints() >= 3 * 10) |
56 |
order = 3;
|
57 |
else if(lyr.getCountPoints() >= 3 * 6) |
58 |
order = 2;
|
59 |
else
|
60 |
order = 1;
|
61 |
|
62 |
|
63 |
affine = new AffineT();
|
64 |
|
65 |
//Calcular la transformaci?n afin por m?nimos cuadrados
|
66 |
affine = computeLeastSquaresAffine(affine, src, dst, order); |
67 |
} |
68 |
|
69 |
/**
|
70 |
* A partir de la transformaci?n af?n creada en el construtor
|
71 |
* genera un fichero de georreferenciaci?n para la imagen.
|
72 |
* @param lyr Capa de puntos
|
73 |
* @param order Orden del sistema de ecuaciones
|
74 |
* @param widthPx Ancho en pixeles de la imagen a georreferenciar
|
75 |
* @param heightPx Alto en pixeles de la imagen a georreferenciar
|
76 |
* @param file Nombre del fichero raster a georreferenciar
|
77 |
*/
|
78 |
public void createWorldFile(int widthPx, int heightPx, String file){ |
79 |
try{
|
80 |
File f = new File(file); |
81 |
String nameWorldFile = f.getPath().substring(0, f.getPath().lastIndexOf(".")) + getExtensionWorldFile(file); |
82 |
createWorldFile(affine, widthPx, heightPx, nameWorldFile); |
83 |
}catch(IOException ex){ |
84 |
System.err.println("Can't create WorldFile"); |
85 |
} |
86 |
} |
87 |
|
88 |
/**
|
89 |
* A partir de la transformaci?n af?n calculada en el contructor
|
90 |
* transforma los puntos pasados como par?metros.
|
91 |
* @param lyr Capa de puntos
|
92 |
* @param list Lista de puntos a transformar
|
93 |
* @return Lista de puntos transformados
|
94 |
*/
|
95 |
public Point2D[] transformPoints(Point2D[] list){ |
96 |
Point2D[] result = new Point2D[list.length]; |
97 |
for(int i = 0; i < list.length;i++){ |
98 |
double[] pto = transformPoint((int)list[i].getX(), (int)list[i].getY(), affine); |
99 |
result[i] = new Point2D.Double(pto[0], pto[1]); |
100 |
} |
101 |
return result;
|
102 |
} |
103 |
|
104 |
/**
|
105 |
* Crea un fichero de georreferenciaci?n (worldfile) a partir de la transformaci?n af?n. Para
|
106 |
* esto necesita obtener las coordenadas reales de la coordenada en pixels (0,0) y calcular
|
107 |
* el tama?o de pixel.
|
108 |
* @param affine Transformaci?n
|
109 |
* @param widthPx Ancho en pixeles de la imagen a georreferenciar
|
110 |
* @param heightPx Alto en pixeles de la imagen a georreferenciar
|
111 |
* @param nameWordlFile Nombre del fichero de georreferenciaci?n
|
112 |
* @throws IOException
|
113 |
*/
|
114 |
private void createWorldFile(AffineT affine, int widthPx, int heightPx, String nameWordlFile) throws IOException { |
115 |
StringBuffer data = new StringBuffer(); |
116 |
double[] min = transformPoint(0, 0, affine); |
117 |
double[] max = transformPoint(widthPx, heightPx, affine); |
118 |
File f = new File(nameWordlFile); |
119 |
DataOutputStream dos = new DataOutputStream( new BufferedOutputStream(new FileOutputStream(f)) ); |
120 |
|
121 |
data.append((max[0] - min[0])/(widthPx - 1)+"\n"); |
122 |
data.append("0.0\n");
|
123 |
data.append("0.0\n");
|
124 |
data.append((max[1] - min[1])/(heightPx - 1)+"\n"); |
125 |
data.append(""+min[0]+"\n"); |
126 |
data.append(""+min[1]+"\n"); |
127 |
|
128 |
dos.writeBytes(data.toString()); |
129 |
dos.close(); |
130 |
} |
131 |
|
132 |
/**
|
133 |
* Obtiene la extensi?n del fichero de georreferenciaci?n
|
134 |
* @return String con la extensi?n del fichero de georreferenciaci?n dependiendo
|
135 |
* del valor del formato obtenido del servidor. Por defecto asignaremos un .wld
|
136 |
*/
|
137 |
private String getExtensionWorldFile(String file){ |
138 |
String ext = file.substring(file.lastIndexOf(".") + 1); |
139 |
String extWorldFile = ".wld"; |
140 |
if(ext.equals("tif") || ext.equals("tiff")) |
141 |
extWorldFile = ".tfw";
|
142 |
if(ext.equals("jpeg") || ext.equals("jpg")) |
143 |
extWorldFile = ".jpgw";
|
144 |
return extWorldFile;
|
145 |
} |
146 |
|
147 |
/**
|
148 |
* Calcula la transformaci?n af?n a partir de los puntos introducidos por m?nimos cuadrados.
|
149 |
* @param af Transformaci?n af?n
|
150 |
* @param src Puntos de entrada en coordenadas pixel
|
151 |
* @param dst Puntos de destino en coordenadas del mundo
|
152 |
* @param order Orden del sistema de ecuaciones
|
153 |
*/
|
154 |
public AffineT computeLeastSquaresAffine(AffineT af, PointT[] src, PointT[] dst, int order){ |
155 |
double[] b = new double[dst.length]; |
156 |
int i,cofs;
|
157 |
|
158 |
af.setOrder(order); |
159 |
cofs = order <= 2 ? order*3 : 10; |
160 |
af.mallocCofs(cofs); |
161 |
|
162 |
//First compute the X cofs
|
163 |
for(i = 0; i < dst.length; i++) { |
164 |
b[i] = dst[i].getX(); |
165 |
} |
166 |
|
167 |
af.setCofX(singleLeastSquaresAffine(af.getXcofs(), src, b, order)); |
168 |
|
169 |
//Now compute the Y cofs
|
170 |
for(i = 0; i < dst.length; i++) { |
171 |
b[i] = dst[i].getY(); |
172 |
} |
173 |
|
174 |
af.setCofY(singleLeastSquaresAffine(af.getYcofs(), src, b, order)); |
175 |
|
176 |
return af;
|
177 |
} |
178 |
|
179 |
private double[] singleLeastSquaresAffine(double[] a, PointT[] src, double[] dst, int order){ |
180 |
int i,j;
|
181 |
int n = dst.length;
|
182 |
|
183 |
int points = order <= 2 ? order * 3 : 10; |
184 |
double[][] combined = new double[points][points + 1]; |
185 |
double[][] A = new double[n + 1][points]; |
186 |
double[][] b = new double[n + 1][1]; |
187 |
|
188 |
for(i = 0; i < n; i++) { |
189 |
A[i][0] = 1.0; |
190 |
A[i][1] = src[i].getX();
|
191 |
A[i][2] = src[i].getY();
|
192 |
if(order > 1) { |
193 |
A[i][3] = src[i].getX() * src[i].getX();
|
194 |
A[i][4] = src[i].getX() * src[i].getY();
|
195 |
A[i][5] = src[i].getY() * src[i].getY();
|
196 |
} |
197 |
if(order > 2) { |
198 |
A[i][6] = src[i].getX() * src[i].getX() * src[i].getX();
|
199 |
A[i][7] = src[i].getX() * src[i].getX() * src[i].getY();
|
200 |
A[i][8] = src[i].getX() * src[i].getY() * src[i].getY();
|
201 |
A[i][9] = src[i].getY() * src[i].getY() * src[i].getY();
|
202 |
} |
203 |
b[i][0] = dst[i];
|
204 |
} |
205 |
|
206 |
double[][] Atrans = MathUtils.transpose(A, n, points); |
207 |
double[][] left = MathUtils.multmatrix(Atrans,A,points,n,points); |
208 |
double[][] right = MathUtils.multmatrix(Atrans,b,points,n,1); |
209 |
|
210 |
for(i = 0; i < points; i++) { |
211 |
combined[i][0] = right[i][0]; |
212 |
for(j = 0; j < points; j++) { |
213 |
combined[i][j+1] = left[i][j];
|
214 |
} |
215 |
} |
216 |
|
217 |
try{
|
218 |
combined = solve(combined, points, points+1);
|
219 |
}catch(DataFormatException ex){ |
220 |
System.err.println("Can't solve matrix"); |
221 |
} |
222 |
|
223 |
for(i = 0; i < points; i++) { |
224 |
a[i] = combined[i][0];
|
225 |
} |
226 |
return a;
|
227 |
} |
228 |
|
229 |
|
230 |
private double[][] solve(double[][] mat,int rows,int cols)throws DataFormatException{ |
231 |
int i,j,k;
|
232 |
double d,tmp;
|
233 |
int big;
|
234 |
|
235 |
for(i = 0; i < rows; i++) { |
236 |
// Find largest row
|
237 |
big = i; |
238 |
for(j = i; j < rows; j++) {
|
239 |
if(Math.abs(mat[j][i+1]) > Math.abs(mat[big][i+1])) { |
240 |
big = j; |
241 |
} |
242 |
} |
243 |
// swap row i and row big
|
244 |
for(k = 0; k < cols ; k++) { |
245 |
tmp = mat[i][k]; |
246 |
mat[i][k] = mat[big][k]; |
247 |
mat[big][k] = tmp; |
248 |
} |
249 |
if(mat[i][i+1] == 0) |
250 |
throw new DataFormatException(); |
251 |
|
252 |
d = 1.0 / mat[i][i+1]; |
253 |
for(j = 0; j < cols ; j++) { |
254 |
mat[i][j] *= d; |
255 |
//assert(!isnan(mat[i][j]));
|
256 |
} |
257 |
for(k = 0; k < rows; k++) { |
258 |
if(k == i)
|
259 |
continue;
|
260 |
if(mat[k][i+1] != 0.0) { |
261 |
d = mat[k][i+1] / mat[i][i+1]; |
262 |
for(j = 0; j < cols ; j++) { |
263 |
mat[k][j] -= d*mat[i][j]; |
264 |
//assert(!isnan(mat[k][j]));
|
265 |
} |
266 |
} |
267 |
} |
268 |
} |
269 |
return mat;
|
270 |
} |
271 |
|
272 |
/**
|
273 |
* Obtiene las coordenadas de un punto en coordenadas del mundo a partir de una transformaci?n
|
274 |
* y el punto de entrada en coordenadas de la imagen.
|
275 |
* @param sx Coordenada X del punto de entrada
|
276 |
* @param syz Coordenada Y del punto de entrada
|
277 |
* @param af Transformaci?n
|
278 |
* @return Punto transformado
|
279 |
*/
|
280 |
public double[] transformPoint(int sx, int sy, AffineT af){ |
281 |
double[] p = new double[2]; |
282 |
if(af.getOrder() == 1) { |
283 |
p[0] = af.getCofX(0) + af.getCofX(1) * sx + af.getCofX(2) * sy; |
284 |
p[1] = af.getCofY(0) + af.getCofY(1) * sx + af.getCofY(2) * sy; |
285 |
} |
286 |
else if(af.getOrder() == 2) { |
287 |
p = quadTransformPoint(sx, sy, af); |
288 |
} |
289 |
else {
|
290 |
p = cubicTransformPoint(sx, sy, af); |
291 |
} |
292 |
return p;
|
293 |
} |
294 |
|
295 |
private static double[] quadTransformPoint(int sx, int sy, AffineT af){ |
296 |
double[] p = new double[2]; |
297 |
p[0] = af.getCofX(0) + af.getCofX(1) * sx + af.getCofX(2) * sy + |
298 |
af.getCofX(3) * sx * sx + af.getCofX(4) * sx * sy + af.getCofX(5) * sy * sy; |
299 |
p[1] = af.getCofY(0) + af.getCofY(1) * sx + af.getCofY(2) * sy + |
300 |
af.getCofY(3) * sx * sx + af.getCofY(4) * sx * sy + af.getCofY(5) * sy * sy; |
301 |
return p;
|
302 |
} |
303 |
|
304 |
private static double[] cubicTransformPoint(int sx, int sy, AffineT af){ |
305 |
double[] p = new double[2]; |
306 |
p[0] = af.getCofX(0) + af.getCofX(1) * sx + af.getCofX(2) * sy + |
307 |
af.getCofX(3) * sx * sx + af.getCofX(4) * sx * sy + af.getCofX(5) * sy * sy + |
308 |
af.getCofX(6) * sx * sx * sx + af.getCofX(7) * sx * sx * sy + |
309 |
af.getCofX(8)*sx*sy*sy + af.getCofX(9)*sy*sy*sy; |
310 |
p[1] = af.getCofY(0) + af.getCofY(1) * sx + af.getCofY(2) * sy + |
311 |
af.getCofY(3) * sx * sx + af.getCofY(4) * sx * sy + af.getCofY(5) * sy * sy + |
312 |
af.getCofY(6) * sx * sx * sx + af.getCofY(7) * sx * sx * sy + |
313 |
af.getCofY(8) * sx * sy * sy + af.getCofY(9) * sy * sy * sy; |
314 |
return p;
|
315 |
} |
316 |
} |