| gvSIG desktop 1

/* gvSIG. Sistema de Informaci?n Geogr?fica de la Generalitat Valenciana

         *

         * Copyright (C) 2006 Instituto de Desarrollo Regional and Generalitat Valenciana.

         *

         * This program is free software; you can redistribute it and/or

         * modify it under the terms of the GNU General Public License

         * as published by the Free Software Foundation; either version 2

         * of the License, or (at your option) any later version.

         *

         * This program is distributed in the hope that it will be useful,

         * but WITHOUT ANY WARRANTY; without even the implied warranty of

         * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the

         * GNU General Public License for more details.

         *

         * You should have received a copy of the GNU General Public License

         * along with this program; if not, write to the Free Software

         * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307,USA.

         *

         * For more information, contact:

         *

         *  Generalitat Valenciana

         *   Conselleria d'Infraestructures i Transport

         *   Av. Blasco Iba?ez, 50

         *   46010 VALENCIA

         *   SPAIN

         *

         *      +34 963862235

         *   gvsig@gva.es

         *      www.gvsig.gva.es

         *

         *    or

         *

         *   Instituto de Desarrollo Regional (Universidad de Castilla La-Mancha)

         *   Campus Universitario s/n

         *   02071 Alabacete

         *   Spain

         *

         *   +34 967 599 200

         */

package org.gvsig.remotesensing.georeferencing.geotransform;

import org.gvsig.rastertools.RasterProcess;

import com.iver.andami.PluginServices;

import Jama.Matrix;

/**

 *  Clase que representa una transformacion polinomial  para la convertir las

 *  coordenadas de una imagen en la imagen rectificada.

 *  

 *  

 *  @author Alejandro Mu?oz Sanchez (alejandro.munoz@uclm.es)

 *         @author Diego Guerrero Sevilla (diego.guerrero@uclm.es)  

 *         @version 20/1/2008

 **/

public class GeoTransformProcess extends RasterProcess {

        // Esto va a cambiar dependiendo del formato de los punto de la interfaz

        // de momento geo_points puntos sobre imagen georeferenciada, image_points

        // puntos sobre imagen a georeferenciar.

        // Utilizar GeoPoint

        protected double geo_points[][]=new double[2][];

        protected double image_points[][]=new double[2][];

        // porcentage del proceso global de la tarea

        private int percent=0;

        // orden del polinomio aproximador

        protected int orden = 0;

        // numero minimo de puntos 

        protected int minGPC=0;

        //coeficientes del polinomio de transformacion para las coordenadas en X

        private double []coefX=null;

    //coeficientes del polinomio de transformacion para las coordenadas en Y

        private double []coefY=null;

        // Lista con los valores de las potencias de x e y  

        private int exp[][] = null;

        // error de la funcion en x

        private double[] xError= null;

        // error de la funcion para las y

        private double [] yError=null;

        // rms

        private double rms[]=null;

        /** Metodo que recoge los parametros del proceso de clasificacion de 

        * m?xima probabilidad 

        * <LI>geo_points: lista de puntos con coordenadas reales</LI>

        * <LI> image_points: lista de puntos coordenadas pixel</LI> 

        * <LI> orden: orden del polinomio de transformacion </LI> 

        */

        public void init() {                

                geo_points =(double[][]) getParam("geoPoints");

                image_points=(double[][]) getParam("imagePoints");

                orden= (int)getIntParam("orden");

                minGPC = (orden+1)*(orden+2)/2;

            exp=new int[minGPC][2];

                // Chequear si el numero de puntos es suficiente para determinar la transformacion de orden n. 

                if(!enoughPoints()){

                        // NOTIFICAR, NO SUFICIENTES PUNTOS PARA ORDEN SELECCIONADO

                        minGPC=0;

                        // Detener ejecucion del calculo.

                }

        }

        /**

         * @return true si se proporciona el numero minimo de puntos de entrada 

         * para la transformaci?n. false en caso contrario.

         * */                

        public boolean enoughPoints() {

                // obligar a que las listas de puntos tengan la misma longitud

                return (geo_points.length>=minGPC && image_points.length==geo_points.length);

        }

        /**

         *  Proceso de calculo de los coeficientes del polinomio. 

         **/

        public void process() throws InterruptedException {

                // Obtencion  polinomio de transformacion en x

                getPolinomyalCoefX();

                // Obtencion del polinomio de transformaci?n en y

                getPolinomialCoefY();

                // calculo de los resultados

                getRMSerror();

        }

        /**

         * @return coeficientes para el polinomio de transformaci?n en x.

         **/

        public double[] getPolinomyalCoefX(){

                if(coefX==null)

                        setDxGeneral();

                return coefX;

        }

        /**

         * @return coeficientes para el polinomio de transformaci?n en y.

         * */

        public double[] getPolinomialCoefY(){

                if(coefY==null)

                        setDyGeneral();

                return coefY;

        }

/////////!!!!!!! Obtiene matriz general a partir de los puntos de control y el orden del polinomio aproximador

        public void setDxGeneral(){

                double matrixDx[][]= new double [minGPC][minGPC]; 

                double result[]= new double[minGPC];

                int exp[][]=new int[minGPC][2];

                int k=-1;

                // Obtencion de la primera fila de la matriz

                for (int filas=0; filas<minGPC;filas++)

                {        k=-1;

                for (int i=0; i<=orden; i++)

                        for(int j=0; j<=i; j++){

                                k++;

                                for(int v=0; v<geo_points.length;v++)

                                        matrixDx[filas][k]+=(Math.pow(geo_points[v][0],i-j)* Math.pow(geo_points[v][0],exp[filas][0]))

                                                        * (Math.pow(geo_points[v][1],j)* Math.pow(geo_points[v][1],exp[filas][1]));                        

                                // Para la fila 0 se guardan los exponentes

                                if(filas==0){

                                        exp[k][0]=i-j;

                                        exp[k][1]=j;

                                        // Se calcula el resultado de !!!!!

                                        for(int v=0; v<geo_points.length;v++)

                                                result[k]+=(Math.pow(geo_points[v][0],i-j)* Math.pow(geo_points[v][0],exp[filas][0]))

                                                * (Math.pow(geo_points[v][1],j)* Math.pow(geo_points[v][1],exp[filas][1]))*image_points[v][0];                                

                                }

                        }

                }

                Matrix matrixResult= new Matrix(matrixDx);

                coefX=solveSistem(matrixResult,result);

        }

        // TO DO: Ver la manera de unificar con setDxGeneral(Parametrizar un metodo general)..... Estudiar optimizaciones!!! 

        // Cambios necesarios para el caolculo de coeficientes para coordenadas y

        public void setDyGeneral(){

                double matrixDy[][]= new double [minGPC][minGPC]; 

                double result[]= new double[minGPC];

                int k=-1;

                // Obtencion de la primera fila de la matriz

                for (int filas=0; filas<minGPC;filas++)

                {        k=-1;

                for (int i=0; i<=orden; i++)

                        for(int j=0; j<=i; j++){

                                k++;

                                for(int v=0; v<geo_points.length;v++)

                                        matrixDy[filas][k]+=(Math.pow(geo_points[v][0],i-j)* Math.pow(geo_points[v][0],exp[filas][0]))

                                                        * (Math.pow(geo_points[v][1],j)* Math.pow(geo_points[v][1],exp[filas][1]));                        

                                // Para la fila 0 se guardan los exponentes

                                if(filas==0){

                                        exp[k][0]=i-j;

                                        exp[k][1]=j;

                                        // Se calcula el resultado de !!!!!

                                        for(int v=0; v<geo_points.length;v++)

                                                result[k]+=(Math.pow(geo_points[v][0],i-j)* Math.pow(geo_points[v][0],exp[filas][0]))

                                                * (Math.pow(geo_points[v][1],j)* Math.pow(geo_points[v][1],exp[filas][1]))*image_points[v][1];                                

                                }

                        }

                }

                Matrix matrixResult= new Matrix(matrixDy);

                coefY=solveSistem(matrixResult,result);        

        }

        /**

         * @return array con la solucion al sistema de ecuadiones.

         * */

        public double[] solveSistem(Matrix matrix, double columResult[]){

                double xCoef []= new double[minGPC];

                double [][]a= new double[columResult.length][1]; 

                for (int i=0; i<columResult.length;i++)

                        a[i][0]=columResult[i];

                Matrix c= matrix.solve(new Matrix(a));

                for (int i=0; i<columResult.length;i++)

                        xCoef[i]=c.get(i,0);

                return xCoef;

        }

        /**

         * @return vector con los RMS 

         * */

        public double[] getRMSerror(){

                double xEvaluate[]= new double [image_points.length];

                double yEvaluate[]= new double [image_points.length];

                rms = new double [image_points.length];

                xError= new double [image_points.length];

                yError= new double[image_points.length];

                for(int im_point=0; im_point<image_points.length;im_point++){

                        for(int i=0; i<minGPC;i++)

                        {

                                xEvaluate[im_point]+=coefX[i] * Math.pow(geo_points[im_point][0], exp[i][0]) * Math.pow(geo_points[im_point][1], exp[i][1]);

                                yEvaluate[im_point]+=coefY[i] * Math.pow(geo_points[im_point][0], exp[i][0]) * Math.pow(geo_points[im_point][1], exp[i][1]);

                        }        

                        xError[im_point]= Math.pow(xEvaluate[im_point] -image_points[im_point][0], 2);

                        yError[im_point]= Math.pow(yEvaluate[im_point] -image_points[im_point][1], 2);

                        rms[im_point]=Math.sqrt

                        ( 

                                        xError[im_point]+ yError[im_point]                

                        );

                }

                System.out.print("Base X\t\t");

                System.out.print("Base Y\t\t");

                System.out.print("WarpX\t\t");

                System.out.print("WarpY\t\t");

                System.out.print("PredicX\t\t\t\t");

                System.out.print("PredicY\t\t\t\t");

                System.out.print("ErrorX\t\t\t\t");

                System.out.print("ErrorY\t\t\t\t");

                System.out.print("RMS");

                // Escriir resultados

                for(int i=0; i<image_points.length;i++)

                        {

                                System.out.print("\n");

                                System.out.print((new Double(geo_points[i][0]).toString()+"\t\t"));

                                System.out.print((new Double(geo_points[i][1]).toString()+"\t\t"));

                                System.out.print((new Double(image_points[i][0]).toString()+"\t\t"));

                                System.out.print((new Double(image_points[i][1]).toString()+"\t\t"));

                                System.out.print((new Double(xEvaluate[i]).toString()+"\t\t"));

                                System.out.print((new Double(yEvaluate[i]).toString()+"\t\t"));

                                System.out.print((new Double(xError[i]).toString()+"\t\t"));

                                System.out.print((new Double(yError[i]).toString()+"\t\t"));

                                System.out.print((new Double(rms[i]).toString()+"\t\t"));

                        }

                return rms;        

        }

        /**

         * @return array con el error en la coordenada x para los puntos de entrada 

         * 

         * */

        public double[] getxError(){

                return xError;

        }

        /**

         * @return array con el error en la coordenada y para los puntos de entrada 

         * 

         * */

        public double[] getyError(){

                return xError;

        }

        /**

         *  Funci?n que evalua el polinomio de transformaci?n para un punto especifico

         * @param x coordenada x del punto

         * @param y coordenada y del punto

         * @return resultado de la evaluaci?n

         * */

        double[] evaluate(double x, double y){

                double eval[]= new double [2];        

                for(int i=0; i<minGPC;i++)

                {

                        eval[0]+=coefX[i] * Math.pow(x, exp[i][0]) * Math.pow(y, exp[i][1]);

                        eval[1]+=coefY[i] * Math.pow(x, exp[i][0]) * Math.pow(y, exp[i][1]);

                }        

                return eval;

        }

        public String getTitle() {

                return PluginServices.getText(this,"transformacion");

        }

        public String getLog() {

                return PluginServices.getText(this,"calculando_transformacion");

        }

        public int getPercent() {

                // TODO Auto-generated method stub

                return percent;

        }

}