root / trunk / extensions / extGeoreferencing / src / org / gvsig / georeferencing / process / geotransform / GeoTransformProcess.java @ 20914
History | View | Annotate | Download (11 KB)
| 1 | 18530 | nbrodin | /* gvSIG. Sistema de Informaci?n Geogr?fica de la Generalitat Valenciana
|
|---|---|---|---|
| 2 | *
|
||
| 3 | * Copyright (C) 2006 Instituto de Desarrollo Regional and Generalitat Valenciana.
|
||
| 4 | *
|
||
| 5 | * This program is free software; you can redistribute it and/or
|
||
| 6 | * modify it under the terms of the GNU General Public License
|
||
| 7 | * as published by the Free Software Foundation; either version 2
|
||
| 8 | * of the License, or (at your option) any later version.
|
||
| 9 | *
|
||
| 10 | * This program is distributed in the hope that it will be useful,
|
||
| 11 | * but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||
| 12 | * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||
| 13 | * GNU General Public License for more details.
|
||
| 14 | *
|
||
| 15 | * You should have received a copy of the GNU General Public License
|
||
| 16 | * along with this program; if not, write to the Free Software
|
||
| 17 | * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307,USA.
|
||
| 18 | *
|
||
| 19 | * For more information, contact:
|
||
| 20 | *
|
||
| 21 | * Generalitat Valenciana
|
||
| 22 | * Conselleria d'Infraestructures i Transport
|
||
| 23 | * Av. Blasco Iba?ez, 50
|
||
| 24 | * 46010 VALENCIA
|
||
| 25 | * SPAIN
|
||
| 26 | *
|
||
| 27 | * +34 963862235
|
||
| 28 | * gvsig@gva.es
|
||
| 29 | * www.gvsig.gva.es
|
||
| 30 | *
|
||
| 31 | * or
|
||
| 32 | *
|
||
| 33 | * Instituto de Desarrollo Regional (Universidad de Castilla La-Mancha)
|
||
| 34 | * Campus Universitario s/n
|
||
| 35 | * 02071 Alabacete
|
||
| 36 | * Spain
|
||
| 37 | *
|
||
| 38 | * +34 967 599 200
|
||
| 39 | */
|
||
| 40 | |||
| 41 | package org.gvsig.georeferencing.process.geotransform; |
||
| 42 | 18673 | amunoz | |
| 43 | 18530 | nbrodin | import org.gvsig.raster.RasterProcess; |
| 44 | 18611 | nbrodin | import org.gvsig.raster.datastruct.GeoPoint; |
| 45 | import org.gvsig.raster.util.RasterToolsUtil; |
||
| 46 | 18530 | nbrodin | |
| 47 | import Jama.Matrix; |
||
| 48 | |||
| 49 | /**
|
||
| 50 | * Clase que representa una transformacion polinomial para la convertir las
|
||
| 51 | * coordenadas de una imagen en la imagen rectificada.
|
||
| 52 | *
|
||
| 53 | *
|
||
| 54 | * @author Alejandro Mu?oz Sanchez (alejandro.munoz@uclm.es)
|
||
| 55 | * @version 20/1/2008
|
||
| 56 | **/
|
||
| 57 | public class GeoTransformProcess extends RasterProcess { |
||
| 58 | |||
| 59 | 19175 | dguerrero | // Lista de puntos de control.
|
| 60 | 18530 | nbrodin | private GeoPoint gpcs[]= null; |
| 61 | |||
| 62 | // porcentage del proceso global de la tarea
|
||
| 63 | private int percent=0; |
||
| 64 | |||
| 65 | // orden del polinomio aproximador
|
||
| 66 | protected int orden = 0; |
||
| 67 | |||
| 68 | // numero minimo de puntos
|
||
| 69 | protected int minGPC=0; |
||
| 70 | |||
| 71 | // Lista con los valores de las potencias de x e y
|
||
| 72 | private int exp[][] = null; |
||
| 73 | |||
| 74 | // rms total en las x
|
||
| 75 | private double rmsXTotal=0; |
||
| 76 | |||
| 77 | // rms total en las y
|
||
| 78 | private double rmsYTotal=0; |
||
| 79 | |||
| 80 | // rms total
|
||
| 81 | private double rmsTotal=0; |
||
| 82 | |||
| 83 | GeoTransformDataResult resultData= null;
|
||
| 84 | |||
| 85 | 18673 | amunoz | //coeficientes polinomios conversion coordenadas reales a coordenadas pixel
|
| 86 | private double mapToPixelCoefX[]= null; |
||
| 87 | private double mapToPixelCoefY[]=null; |
||
| 88 | |||
| 89 | |||
| 90 | //coeficientes polinomio conversion coordenadas pixel a coordenadas reales
|
||
| 91 | private double pixelToMapCoefX[]= null; |
||
| 92 | private double pixelToMapCoefY[]=null; |
||
| 93 | |||
| 94 | 18530 | nbrodin | /**
|
| 95 | * Metodo que recoge los parametros del proceso de transformaci?n
|
||
| 96 | * <LI>gpcs: lista de puntos de control</LI>>
|
||
| 97 | * <LI> orden: orden del polinomio de transformacion </LI>
|
||
| 98 | */
|
||
| 99 | public void init() { |
||
| 100 | gpcs = (GeoPoint[])getParam("gpcs"); |
||
| 101 | orden= (int)getIntParam("orden"); |
||
| 102 | minGPC = (orden+1)*(orden+2)/2; |
||
| 103 | 18611 | nbrodin | exp = new int[minGPC][2]; |
| 104 | resetErrors(); |
||
| 105 | 18530 | nbrodin | resultData= new GeoTransformDataResult();
|
| 106 | // Chequear si el numero de puntos es suficiente para determinar la transformacion de orden n.
|
||
| 107 | if(!enoughPoints()){
|
||
| 108 | minGPC=0;
|
||
| 109 | } |
||
| 110 | } |
||
| 111 | |||
| 112 | /**
|
||
| 113 | 18611 | nbrodin | * Inicializa los valores de los errores a 0
|
| 114 | */
|
||
| 115 | private void resetErrors() { |
||
| 116 | for (int i = 0; i < gpcs.length; i++) |
||
| 117 | gpcs[i].resetErrors(); |
||
| 118 | } |
||
| 119 | |||
| 120 | /**
|
||
| 121 | 18530 | nbrodin | * @return true si se proporciona el numero minimo de puntos de entrada
|
| 122 | * para la transformaci?n de orden seleccionado. false en caso contrario.
|
||
| 123 | * */
|
||
| 124 | public boolean enoughPoints() { |
||
| 125 | return (gpcs.length>=minGPC);
|
||
| 126 | } |
||
| 127 | |||
| 128 | 18590 | nbrodin | /**
|
| 129 | * Obtiene el resultado del proceso.
|
||
| 130 | * @return
|
||
| 131 | */
|
||
| 132 | public Object getResult() { |
||
| 133 | return resultData;
|
||
| 134 | } |
||
| 135 | 18530 | nbrodin | |
| 136 | /**
|
||
| 137 | * Proceso
|
||
| 138 | **/
|
||
| 139 | public void process() { |
||
| 140 | 18611 | nbrodin | if(minGPC > 0) { |
| 141 | 19022 | amunoz | // Obtencion los polinomios de transformaci?n
|
| 142 | calculatePolinomialCoefs(); |
||
| 143 | 18530 | nbrodin | // calculo de los resultados
|
| 144 | calculateRMSerror(); |
||
| 145 | // Se almacenan los resultados en dataResult
|
||
| 146 | resultData.setGpcs(gpcs); |
||
| 147 | 18673 | amunoz | resultData.setPixelToMapCoefX(pixelToMapCoefX); |
| 148 | resultData.setPixelToMapCoefY(pixelToMapCoefY); |
||
| 149 | resultData.setMapToPixelCoefX(mapToPixelCoefX); |
||
| 150 | resultData.setMapToPixelCoefY(mapToPixelCoefY); |
||
| 151 | 18530 | nbrodin | resultData.setRmsXTotal(rmsXTotal); |
| 152 | resultData.setRmsYTotal(rmsYTotal); |
||
| 153 | resultData.setRmsTotal(rmsTotal); |
||
| 154 | resultData.setPolynomialOrden(orden); |
||
| 155 | |||
| 156 | if(externalActions!=null) |
||
| 157 | externalActions.end(resultData); |
||
| 158 | 18590 | nbrodin | return;
|
| 159 | 18530 | nbrodin | } |
| 160 | 18590 | nbrodin | resultData = null;
|
| 161 | 18530 | nbrodin | } |
| 162 | |||
| 163 | |||
| 164 | |||
| 165 | /**
|
||
| 166 | 19022 | amunoz | * Calculo de los coeficientes de los polinomios de aproximacion.
|
| 167 | 18530 | nbrodin | * */
|
| 168 | 19022 | amunoz | public void calculatePolinomialCoefs(){ |
| 169 | double matrixD[][]= new double [minGPC][minGPC]; |
||
| 170 | double matrixD2[][]= new double [minGPC][minGPC]; |
||
| 171 | 18530 | nbrodin | double result[]= new double[minGPC]; |
| 172 | 18673 | amunoz | double result2[]= new double[minGPC]; |
| 173 | 19022 | amunoz | double result3[]= new double[minGPC]; |
| 174 | double result4[]= new double[minGPC]; |
||
| 175 | 18530 | nbrodin | int k=-1; |
| 176 | // Obtencion de la primera fila de la matriz
|
||
| 177 | for (int filas=0; filas<minGPC;filas++) |
||
| 178 | { k=-1;
|
||
| 179 | for (int i=0; i<=orden; i++) |
||
| 180 | for(int j=0; j<=i; j++){ |
||
| 181 | k++; |
||
| 182 | 18673 | amunoz | for(int v=0; v<gpcs.length;v++){ |
| 183 | 19022 | amunoz | if(gpcs[v].active){
|
| 184 | matrixD[filas][k]+=(Math.pow(gpcs[v].pixelPoint.getX(),i-j)* Math.pow(gpcs[v].pixelPoint.getX(),exp[filas][0])) |
||
| 185 | * (Math.pow(gpcs[v].pixelPoint.getY(),j)* Math.pow(gpcs[v].pixelPoint.getY(),exp[filas][1])); |
||
| 186 | matrixD2[filas][k]+=(Math.pow(gpcs[v].mapPoint.getX(),i-j)* Math.pow(gpcs[v].mapPoint.getX(),exp[filas][0])) |
||
| 187 | * (Math.pow(gpcs[v].mapPoint.getY(),j)* Math.pow(gpcs[v].mapPoint.getY(),exp[filas][1])); |
||
| 188 | } |
||
| 189 | 18673 | amunoz | |
| 190 | } |
||
| 191 | 18530 | nbrodin | // Para la fila 0 se guardan los exponentes
|
| 192 | if(filas==0){ |
||
| 193 | exp[k][0]=i-j;
|
||
| 194 | exp[k][1]=j;
|
||
| 195 | |||
| 196 | // Se calcula el resultado de !!!!!
|
||
| 197 | 18673 | amunoz | for(int v=0; v<gpcs.length;v++){ |
| 198 | 19022 | amunoz | if(gpcs[v].active){
|
| 199 | result[k]+=(Math.pow(gpcs[v].pixelPoint.getX(),i-j)* Math.pow(gpcs[v].pixelPoint.getX(),exp[filas][0])) |
||
| 200 | 18673 | amunoz | * (Math.pow(gpcs[v].pixelPoint.getY(),j)* Math.pow(gpcs[v].pixelPoint.getY(),exp[filas][1]))*gpcs[v].mapPoint.getY(); |
| 201 | |||
| 202 | 19022 | amunoz | result2[k]+=(Math.pow(gpcs[v].mapPoint.getX(),i-j)* Math.pow(gpcs[v].mapPoint.getX(),exp[filas][0])) |
| 203 | 18673 | amunoz | * (Math.pow(gpcs[v].mapPoint.getY(),j)* Math.pow(gpcs[v].mapPoint.getY(),exp[filas][1]))*gpcs[v].pixelPoint.getY(); |
| 204 | 19022 | amunoz | |
| 205 | result3[k]+=(Math.pow(gpcs[v].pixelPoint.getX(),i-j)* Math.pow(gpcs[v].pixelPoint.getX(),exp[filas][0])) |
||
| 206 | * (Math.pow(gpcs[v].pixelPoint.getY(),j)* Math.pow(gpcs[v].pixelPoint.getY(),exp[filas][1]))*gpcs[v].mapPoint.getX(); |
||
| 207 | |||
| 208 | result4[k]+=(Math.pow(gpcs[v].mapPoint.getX(),i-j)* Math.pow(gpcs[v].mapPoint.getX(),exp[filas][0])) |
||
| 209 | * (Math.pow(gpcs[v].mapPoint.getY(),j)* Math.pow(gpcs[v].mapPoint.getY(),exp[filas][1]))*gpcs[v].pixelPoint.getX(); |
||
| 210 | |||
| 211 | |||
| 212 | } |
||
| 213 | 18673 | amunoz | } |
| 214 | 18530 | nbrodin | } |
| 215 | } |
||
| 216 | } |
||
| 217 | 20876 | amunoz | |
| 218 | 19022 | amunoz | Matrix matrixResult= new Matrix(matrixD);
|
| 219 | Matrix matrixResult2= new Matrix(matrixD2);
|
||
| 220 | 18788 | amunoz | pixelToMapCoefY=solveSystem(matrixResult,result); |
| 221 | mapToPixelCoefY= solveSystem(matrixResult2,result2); |
||
| 222 | 19022 | amunoz | pixelToMapCoefX=solveSystem(matrixResult,result3); |
| 223 | mapToPixelCoefX= solveSystem(matrixResult2,result4); |
||
| 224 | 18530 | nbrodin | } |
| 225 | 18673 | amunoz | |
| 226 | |||
| 227 | 18760 | amunoz | |
| 228 | 18530 | nbrodin | /**
|
| 229 | * @return array con la solucion al sistema de ecuadiones.
|
||
| 230 | * */
|
||
| 231 | 18788 | amunoz | public double[] solveSystem(Matrix matrix, double columResult[]){ |
| 232 | 18530 | nbrodin | double xCoef []= new double[minGPC]; |
| 233 | double [][]a= new double[columResult.length][1]; |
||
| 234 | for (int i=0; i<columResult.length;i++) |
||
| 235 | a[i][0]=columResult[i];
|
||
| 236 | 20876 | amunoz | Matrix c=null;
|
| 237 | if(matrix.det()==0.0){ |
||
| 238 | 20914 | amunoz | // Resolucion del sistema usando la libreria flanagan
|
| 239 | flanagan.math.Matrix matrixFL= new flanagan.math.Matrix(matrix.getArray());
|
||
| 240 | xCoef= matrixFL.solveLinearSet(columResult); |
||
| 241 | 20876 | amunoz | } |
| 242 | 20914 | amunoz | else{
|
| 243 | 20876 | amunoz | c= matrix.solve(new Matrix(a));
|
| 244 | 20914 | amunoz | for (int i=0; i<columResult.length;i++) |
| 245 | xCoef[i]=c.get(i,0);
|
||
| 246 | } |
||
| 247 | 18530 | nbrodin | return xCoef;
|
| 248 | } |
||
| 249 | |||
| 250 | |||
| 251 | /**
|
||
| 252 | * @return vector con los RMS
|
||
| 253 | * */
|
||
| 254 | 18611 | nbrodin | public void calculateRMSerror() { |
| 255 | 18774 | amunoz | |
| 256 | 18530 | nbrodin | int numgpcs= gpcs.length;
|
| 257 | double rmsxTotal=0; |
||
| 258 | double rmsyTotal=0; |
||
| 259 | 18774 | amunoz | rmsTotal=0;
|
| 260 | 19029 | amunoz | int num=0; |
| 261 | 18530 | nbrodin | |
| 262 | 18611 | nbrodin | for(int im_point = 0; im_point < numgpcs; im_point++) { |
| 263 | 19022 | amunoz | if (gpcs[im_point].active)
|
| 264 | {
|
||
| 265 | double tmp[]= getCoordMap(gpcs[im_point].pixelPoint.getX(),gpcs[im_point].pixelPoint.getY()); |
||
| 266 | double tmp2[]= getCoordPixel(tmp[0],tmp[1]); |
||
| 267 | gpcs[im_point].setEvaluateX(tmp2[0]);
|
||
| 268 | gpcs[im_point].setEvaluateY(tmp2[1]);
|
||
| 269 | gpcs[im_point].setErrorX(Math.pow(gpcs[im_point].getEvaluateX() - gpcs[im_point].pixelPoint.getX(), 2)); |
||
| 270 | rmsxTotal += gpcs[im_point].getErrorX(); |
||
| 271 | gpcs[im_point].setErrorY(Math.pow(gpcs[im_point].getEvaluateY() - gpcs[im_point].pixelPoint.getY(), 2)); |
||
| 272 | rmsyTotal +=gpcs[im_point].getErrorY(); |
||
| 273 | gpcs[im_point].setRms(Math.sqrt(gpcs[im_point].getErrorX() +gpcs[im_point].getErrorY()));
|
||
| 274 | rmsTotal += gpcs[im_point].getRms(); |
||
| 275 | 19029 | amunoz | num++; |
| 276 | 19022 | amunoz | } |
| 277 | 18530 | nbrodin | } |
| 278 | 18774 | amunoz | |
| 279 | 19029 | amunoz | this.rmsTotal =Math.sqrt( rmsTotal / num); |
| 280 | this.rmsXTotal = Math.sqrt(rmsxTotal / num); |
||
| 281 | this.rmsYTotal = Math.sqrt(rmsyTotal / num); |
||
| 282 | 18530 | nbrodin | |
| 283 | } |
||
| 284 | |||
| 285 | /**
|
||
| 286 | * @return error total para la coordenada X
|
||
| 287 | * */
|
||
| 288 | public double getRMSx(){ |
||
| 289 | return rmsXTotal;
|
||
| 290 | } |
||
| 291 | |||
| 292 | |||
| 293 | /**
|
||
| 294 | * @return error total para la coordenada y
|
||
| 295 | * */
|
||
| 296 | 18611 | nbrodin | public double getRMSy() { |
| 297 | 18530 | nbrodin | return rmsYTotal;
|
| 298 | } |
||
| 299 | |||
| 300 | /**
|
||
| 301 | 18768 | amunoz | * @return error total
|
| 302 | 18530 | nbrodin | * */
|
| 303 | 18611 | nbrodin | public double getRMSTotal() { |
| 304 | 18530 | nbrodin | return rmsTotal;
|
| 305 | } |
||
| 306 | |||
| 307 | 18611 | nbrodin | /*
|
| 308 | * (non-Javadoc)
|
||
| 309 | * @see org.gvsig.gui.beans.incrementabletask.IIncrementable#getTitle()
|
||
| 310 | */
|
||
| 311 | 18530 | nbrodin | public String getTitle() { |
| 312 | 18611 | nbrodin | return RasterToolsUtil.getText(this, "transformacion"); |
| 313 | 18530 | nbrodin | } |
| 314 | |||
| 315 | 18611 | nbrodin | /*
|
| 316 | * (non-Javadoc)
|
||
| 317 | * @see org.gvsig.raster.RasterProcess#getLog()
|
||
| 318 | */
|
||
| 319 | 18530 | nbrodin | public String getLog() { |
| 320 | 18611 | nbrodin | return RasterToolsUtil.getText(this, "calculando_transformacion"); |
| 321 | 18530 | nbrodin | } |
| 322 | |||
| 323 | 18611 | nbrodin | /*
|
| 324 | * (non-Javadoc)
|
||
| 325 | * @see org.gvsig.gui.beans.incrementabletask.IIncrementable#getPercent()
|
||
| 326 | */
|
||
| 327 | 18530 | nbrodin | public int getPercent() { |
| 328 | return percent;
|
||
| 329 | } |
||
| 330 | |||
| 331 | |||
| 332 | 18768 | amunoz | |
| 333 | 18760 | amunoz | /**
|
| 334 | * Funci?n que evalua el polinomio de transformaci?n para obtener las coordenadas
|
||
| 335 | * reales de unas coordenadas pixeles que se pasan como parametros.
|
||
| 336 | * @param x coordenada x del punto
|
||
| 337 | * @param y coordenada y del punto
|
||
| 338 | * @return resultado de la evaluaci?n para x e y
|
||
| 339 | * */
|
||
| 340 | public double[] getCoordMap(double x, double y){ |
||
| 341 | 20006 | amunoz | //setExpPolinomial();
|
| 342 | 18760 | amunoz | double eval[]= new double [2]; |
| 343 | for(int i=0; i<pixelToMapCoefX.length;i++) |
||
| 344 | {
|
||
| 345 | eval[0]+=pixelToMapCoefX[i] * Math.pow(x, exp[i][0]) * Math.pow(y,exp[i][1]); |
||
| 346 | eval[1]+=pixelToMapCoefY[i] * Math.pow(x, exp[i][0]) * Math.pow(y, exp[i][1]); |
||
| 347 | } |
||
| 348 | return eval;
|
||
| 349 | } |
||
| 350 | |||
| 351 | |||
| 352 | /**
|
||
| 353 | * Funci?n que evalua el polinomio de transformaci?n para obtener las coordenadas
|
||
| 354 | * pixeles de unas coordenadas mapa en un proceso de transformacion.
|
||
| 355 | * @param x coordenada x del punto
|
||
| 356 | * @param y coordenada y del punto
|
||
| 357 | * @return resultado de la evaluaci?n para x e y
|
||
| 358 | * */
|
||
| 359 | public double[] getCoordPixel(double x, double y){ |
||
| 360 | 20006 | amunoz | //setExpPolinomial();
|
| 361 | 18760 | amunoz | double eval[]= new double [2]; |
| 362 | for(int i=0; i<pixelToMapCoefX.length;i++) |
||
| 363 | {
|
||
| 364 | eval[0]+=mapToPixelCoefX[i] * Math.pow(x, exp[i][0]) * Math.pow(y,exp[i][1]); |
||
| 365 | eval[1]+=mapToPixelCoefY[i] * Math.pow(x, exp[i][0]) * Math.pow(y, exp[i][1]); |
||
| 366 | } |
||
| 367 | return eval;
|
||
| 368 | } |
||
| 369 | |||
| 370 | |||
| 371 | 20006 | amunoz | /*private void setExpPolinomial(){
|
| 372 | 18760 | amunoz | if(exp==null)
|
| 373 | {
|
||
| 374 | int minGPC=(orden+1)*(orden+2)/2;
|
||
| 375 | exp=new int[minGPC][2];
|
||
| 376 | int k=-1;
|
||
| 377 | for (int i=0; i<=orden; i++){
|
||
| 378 | for(int j=0; j<=i; j++){
|
||
| 379 | k++;
|
||
| 380 | exp[k][0]=i-j;
|
||
| 381 | exp[k][1]=j;
|
||
| 382 | }
|
||
| 383 | }
|
||
| 384 | }
|
||
| 385 | 20006 | amunoz | }*/
|
| 386 | 18760 | amunoz | |
| 387 | 18530 | nbrodin | } |