Application: gvSIG desktop

/* 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 Jama.Matrix;

import com.iver.andami.PluginServices;

/**

 *  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.

        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;

        /** 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;

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

                if(!enoughPoints()){

                        // NOTIFICAR, NO SUFICIENTES PUNTOS PARA ORDEN SELECCIONADO

                        return;

                }

        }

        /**

         * @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();

        }

        /**

         * @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);

                try {

                        coefX=solveSistemforCramer(matrixResult,result);

                } catch (BadDeteminantException e) {

                        // Resolver sistema de ecuaciones o por otro metodo

                }

        }

        // 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 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++)

                                        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][0];

                                }

                        }

                }

                Matrix matrixResult= new Matrix(matrixDy);

                try {

                        coefY=solveSistemforCramer(matrixResult,result);

                } catch (BadDeteminantException e) {

                        // Resolver sistema de ecuaciones opor otro metodo

                }

        }

        /**

         *         Resoluci?n del sistema por la regla de Cramer.

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

         * */

        public double[] solveSistemforCramer(Matrix matrix, double columResult[]) throws BadDeteminantException{

                double xCoef []= new double[minGPC];

                double det= matrix.det();

                if(det==0)

                        throw new BadDeteminantException();

                double aux[][]= matrix.getArrayCopy();

                Matrix m=null;

                // Resolucion del sistema por cramer

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

                {

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

                                aux[i][k]= columResult[i];

                        m= new Matrix (aux);

                        xCoef[k]= m.det()/det;

                        aux=matrix.getArrayCopy();;

                }

                return xCoef;

        }

        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;

        }

}