| 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.georeferencing.process.geotransform;

import org.gvsig.raster.RasterProcess;

import org.gvsig.raster.datastruct.GeoPoint;

import org.gvsig.raster.util.RasterToolsUtil;

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)

 *         @version 20/1/2008

 **/

public class GeoTransformProcess extends RasterProcess {

        // Lista de puntos de control

        private GeoPoint gpcs[]= null;

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

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

        private int exp[][] = null;

        // rms total en las x

        private double rmsXTotal=0;

        // rms total en las y

        private double rmsYTotal=0;

        // rms total

        private double rmsTotal=0;

        GeoTransformDataResult resultData= null;

        //coeficientes polinomios conversion coordenadas reales a coordenadas pixel

        private double mapToPixelCoefX[]= null;

        private double mapToPixelCoefY[]=null;

        //coeficientes polinomio conversion coordenadas pixel a coordenadas reales

        private double pixelToMapCoefX[]= null;

        private double pixelToMapCoefY[]=null;

        /** 

        * Metodo que recoge los parametros del proceso de transformaci?n 

        * <LI>gpcs: lista de puntos de control</LI>> 

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

        */

        public void init() {                

                gpcs = (GeoPoint[])getParam("gpcs");

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

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

            exp = new int[minGPC][2];

            resetErrors();

            resultData= new GeoTransformDataResult();

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

                }

        }

        /**

         * Inicializa los valores de los errores a 0

         */

        private void resetErrors() {

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

                        gpcs[i].resetErrors();

        }

        /**

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

         * para la transformaci?n de orden seleccionado. false en caso contrario.

         * */                

        public boolean enoughPoints() {

                return (gpcs.length>=minGPC);

        }

        /**

         * Obtiene el resultado del proceso.

         * @return

         */

        public Object getResult() {

                return resultData;

        }

        /**

         *  Proceso

         **/

        public void process() {

                if(minGPC > 0) {

                        // Obtencion  polinomio de transformacion en x

                        calculatePolinomialCoefX();

                        // Obtencion del polinomio de transformaci?n en y

                        calculatePolinomialCoefY();

                        // calculo de los resultados

                        calculateRMSerror();

                        // Se almacenan los resultados en dataResult

                        resultData.setGpcs(gpcs);

                        resultData.setPixelToMapCoefX(pixelToMapCoefX);

                        resultData.setPixelToMapCoefY(pixelToMapCoefY);

                        resultData.setMapToPixelCoefX(mapToPixelCoefX);

                        resultData.setMapToPixelCoefY(mapToPixelCoefY);

                        resultData.setRmsXTotal(rmsXTotal);

                        resultData.setRmsYTotal(rmsYTotal);

                        resultData.setRmsTotal(rmsTotal);

                        resultData.setPolynomialOrden(orden);

                        if(externalActions!=null)

                                externalActions.end(resultData);

                        return;

                }

                resultData = null;

        }

    /**

     *  Calculo de los coeficientes del polinimio aproximador.

     *  @return array con los coeficientes para las x. 

     * 

     * */

        public void calculatePolinomialCoefX(){

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

                double matrixDx2[][]= new double [minGPC][minGPC];

                double result[]= new double[minGPC];

                double result2[]= 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<gpcs.length;v++){

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

                                                        * (Math.pow(gpcs[v].pixelPoint.getY(),j)* Math.pow(gpcs[v].pixelPoint.getY(),exp[filas][1]));

                                        matrixDx2[filas][k]+=(Math.pow(gpcs[v].mapPoint.getX(),i-j)* Math.pow(gpcs[v].mapPoint.getX(),exp[filas][0]))

                                        * (Math.pow(gpcs[v].mapPoint.getY(),j)* Math.pow(gpcs[v].mapPoint.getY(),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<gpcs.length;v++){

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

                                                * (Math.pow(gpcs[v].pixelPoint.getY(),j)* Math.pow(gpcs[v].pixelPoint.getY(),exp[filas][1]))*gpcs[v].mapPoint.getX();

                                                result2[k]+=(Math.pow(gpcs[v].mapPoint.getX(),i-j)* Math.pow(gpcs[v].mapPoint.getX(),exp[filas][0]))

                                                * (Math.pow(gpcs[v].mapPoint.getY(),j)* Math.pow(gpcs[v].mapPoint.getY(),exp[filas][1]))*gpcs[v].pixelPoint.getX();        

                                        }        

                                }

                        }

                }

                Matrix matrixResult= new Matrix(matrixDx);

                Matrix matrixResult2= new Matrix(matrixDx2);

                pixelToMapCoefX=solveSistem(matrixResult,result);        

                mapToPixelCoefX=solveSistem(matrixResult2,result2);

        }

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

         /**

     *  Calculo de los coeficientes del polinimio aproximador.

     *  @return array con los coeficientes para las x. 

     * 

     * */

        public void calculatePolinomialCoefY(){

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

                double matrixDy2[][]= new double [minGPC][minGPC]; 

                double result[]= new double[minGPC];

                double result2[]= 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<gpcs.length;v++){

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

                                                        * (Math.pow(gpcs[v].pixelPoint.getY(),j)* Math.pow(gpcs[v].pixelPoint.getY(),exp[filas][1]));        

                                        matrixDy2[filas][k]+=(Math.pow(gpcs[v].mapPoint.getX(),i-j)* Math.pow(gpcs[v].mapPoint.getX(),exp[filas][0]))

                                        * (Math.pow(gpcs[v].mapPoint.getY(),j)* Math.pow(gpcs[v].mapPoint.getY(),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<gpcs.length;v++){

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

                                                                * (Math.pow(gpcs[v].pixelPoint.getY(),j)* Math.pow(gpcs[v].pixelPoint.getY(),exp[filas][1]))*gpcs[v].mapPoint.getY();        

                                                result2[k]+=(Math.pow(gpcs[v].mapPoint.getX(),i-j)* Math.pow(gpcs[v].mapPoint.getX(),exp[filas][0]))

                                                                * (Math.pow(gpcs[v].mapPoint.getY(),j)* Math.pow(gpcs[v].mapPoint.getY(),exp[filas][1]))*gpcs[v].pixelPoint.getY();

                                        }        

                                }

                        }

                }

                Matrix matrixResult= new Matrix(matrixDy);

                Matrix matrixResult2= new Matrix(matrixDy2);

                pixelToMapCoefY=solveSistem(matrixResult,result);        

                mapToPixelCoefY=  solveSistem(matrixResult2,result2);

        }

        /**

         * @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 void calculateRMSerror() {

                int numgpcs= gpcs.length;

                double rmsxTotal=0;

                double rmsyTotal=0;

                rmsTotal=0;

                for(int im_point = 0; im_point < numgpcs; im_point++) {

                        double tmp[]= getCoordMap(gpcs[im_point].pixelPoint.getX(),gpcs[im_point].pixelPoint.getY());

                        double tmp2[]= getCoordPixel(tmp[0],tmp[1]);

                        gpcs[im_point].setEvaluateX(tmp2[0]);

                        gpcs[im_point].setEvaluateY(tmp2[1]);

                        gpcs[im_point].setErrorX(Math.pow(gpcs[im_point].getEvaluateX() - gpcs[im_point].pixelPoint.getX(), 2));

                        rmsxTotal += gpcs[im_point].getErrorX();

                        gpcs[im_point].setErrorY(Math.pow(gpcs[im_point].getEvaluateY() - gpcs[im_point].pixelPoint.getY(), 2));

                        rmsyTotal += gpcs[im_point].getErrorY();

                        gpcs[im_point].setRms(Math.sqrt(gpcs[im_point].getErrorX() + gpcs[im_point].getErrorY()));

                        rmsTotal += gpcs[im_point].getRms();

                }

                this.rmsTotal = rmsTotal / numgpcs; 

                this.rmsXTotal = rmsxTotal / numgpcs;

                this.rmsYTotal = rmsyTotal / numgpcs;

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

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

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

                System.out.print("Pix Y\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");

                // Escribir resultados

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

                        {

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

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

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

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

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

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

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

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

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

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

                        }        

                System.out.println(new Double (rmsTotal).toString());

        }

        /**

         *  @return error total para la coordenada X

         * */

        public double getRMSx(){

                return rmsXTotal;        

        }

        /**

         *  @return error total para la coordenada y

         * */

        public double getRMSy() {

                return rmsYTotal;        

        }

        /**

         *  @return error total

         * */

        public double getRMSTotal() {

                return rmsTotal;        

        }

        /*

         * (non-Javadoc)

         * @see org.gvsig.gui.beans.incrementabletask.IIncrementable#getTitle()

         */

        public String getTitle() {

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

        }

        /*

         * (non-Javadoc)

         * @see org.gvsig.raster.RasterProcess#getLog()

         */

        public String getLog() {

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

        }

        /*

         * (non-Javadoc)

         * @see org.gvsig.gui.beans.incrementabletask.IIncrementable#getPercent()

         */

        public int getPercent() {

                return percent;

        }

        /**

         *  Funci?n que evalua el polinomio de transformaci?n para obtener las coordenadas

         *  reales de unas coordenadas pixeles que se pasan como parametros.

         * @param x coordenada x del punto

         * @param y coordenada y del punto

         * @return resultado de la evaluaci?n para x e y

         * */

        public double[] getCoordMap(double x, double y){

                setExpPolinomial();

                double eval[]= new double [2];        

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

                {

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

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

                }        

                return eval;

        }

        /**

         *  Funci?n que evalua el polinomio de transformaci?n para obtener las coordenadas

         *  pixeles de unas coordenadas mapa en un proceso de transformacion.

         * @param x coordenada x del punto

         * @param y coordenada y del punto

         * @return resultado de la evaluaci?n para x e y

         * */

        public double[] getCoordPixel(double x, double y){

                setExpPolinomial();

                double eval[]= new double [2];        

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

                {

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

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

                }        

                return eval;

        }

        private void setExpPolinomial(){

                if(exp==null)

                        {

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

                                exp=new int[minGPC][2];

                                int k=-1;

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

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

                                                k++;

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

                                                exp[k][1]=j;

                                        }

                                }

                        }        

        }

}