Application: gvSIG desktop

package com.iver.cit.gvsig;

import java.awt.geom.Point2D;

import java.io.BufferedOutputStream;

import java.io.DataOutputStream;

import java.io.File;

import java.io.FileOutputStream;

import java.io.IOException;

import java.util.zip.DataFormatException;

import com.iver.cit.gvsig.fmap.layers.FLyrPoints;

import com.iver.cit.gvsig.utils.AffineT;

import com.iver.cit.gvsig.utils.MathUtils;

import com.iver.cit.gvsig.utils.PointT;

/**

 * @author Nacho Brodin (brodin_ign@gva.es)

 */

public class GeoOperations{

        private FLyrPoints         lyrPoints = null;

        private int                         order = 1;

        private PointT[]                 src = null;

        private PointT[]                 dst = null;

        private AffineT                affine = null;

        /**

         * Constructor

         */

        public GeoOperations(){}

        /**

         * Constructor. Crea las estructuras de puntos y calcula la transformaci?n af?n.

         * @param lyr

         */

        public GeoOperations(FLyrPoints lyr){

                this.lyrPoints = lyr;

                src = new PointT[lyr.getCountPoints()];

                dst = new PointT[lyr.getCountPoints()];

                int nPoint = 0;

                for(int i = 0; i<lyr.getCountPoints(); i++){

                        if(lyr.getPoint(i).active == true){

                                src[nPoint] = new PointT();

                                dst[nPoint] = new PointT();

                                src[nPoint].setX(lyr.getPoint(i).pixelPoint.getX());

                                src[nPoint].setY(lyr.getPoint(i).pixelPoint.getY());

                                src[nPoint].setI(lyr.getCountPoints());

                                dst[nPoint].setX(lyr.getPoint(i).mapPoint.getX());

                                dst[nPoint].setY(lyr.getPoint(i).mapPoint.getY());

                                dst[nPoint].setI(lyr.getCountPoints());

                                nPoint++;

                        }

                }

                if(lyr.getCountPoints() >= 3 * 10)

                        order = 3;

                else if(lyr.getCountPoints() >= 3 * 6)

                        order = 2;

                else

                        order = 1;

                affine = new AffineT();

                //Calcular la transformaci?n afin por m?nimos cuadrados

                affine = computeLeastSquaresAffine(affine, src, dst, order);

        }

        /**

         * A partir de la transformaci?n af?n creada en el construtor

         * genera un fichero de georreferenciaci?n para la imagen.

         * @param lyr                Capa de puntos

         * @param order        Orden del sistema de ecuaciones

         * @param widthPx        Ancho en pixeles de la imagen a georreferenciar

         * @param heightPx        Alto en pixeles de la imagen a georreferenciar

         * @param file                Nombre del fichero raster a georreferenciar

         */

        public void createWorldFile(int widthPx, int heightPx, String file){

                try{

                        File f = new File(file);

                        String nameWorldFile = f.getPath().substring(0, f.getPath().lastIndexOf(".")) + getExtensionWorldFile(file);

                        createWorldFile(affine, widthPx, heightPx, nameWorldFile);

                }catch(IOException ex){

                        System.err.println("Can't create WorldFile");

                }

        }

        /**

         * A partir de la transformaci?n af?n calculada en el contructor

         * transforma los puntos pasados como par?metros.

         * @param lyr                Capa de puntos

         * @param list                Lista de puntos a transformar

         * @return                  Lista de puntos transformados

         */

        public Point2D[] transformPoints(Point2D[] list){

                Point2D[] result = new Point2D[list.length];

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

                        double[] pto = transformPoint((int)list[i].getX(), (int)list[i].getY(), affine);

                        result[i] = new Point2D.Double(pto[0], pto[1]);

                }

                return result;

        }

        /**

         * Crea un fichero de georreferenciaci?n (worldfile) a partir de la transformaci?n af?n. Para

         * esto necesita obtener las coordenadas reales de la coordenada en pixels (0,0) y calcular

         * el tama?o de pixel.

         * @param affine                Transformaci?n

         * @param widthPx                Ancho en pixeles de la imagen a georreferenciar

         * @param heightPx                Alto en pixeles de la imagen a georreferenciar

         * @param nameWordlFile        Nombre del fichero de georreferenciaci?n

         * @throws IOException

         */

        private void createWorldFile(AffineT affine, int widthPx, int heightPx, String nameWordlFile) throws IOException {

                StringBuffer data = new StringBuffer();

                double[] min = transformPoint(0, 0, affine);

                double[] max = transformPoint(widthPx, heightPx, affine);

                File f = new File(nameWordlFile);

                DataOutputStream dos = new DataOutputStream( new BufferedOutputStream(new FileOutputStream(f)) );

            data.append((max[0] - min[0])/(widthPx - 1)+"\n");

            data.append("0.0\n");

            data.append("0.0\n");

            data.append((max[1] - min[1])/(heightPx - 1)+"\n");

            data.append(""+min[0]+"\n");

            data.append(""+min[1]+"\n");

            dos.writeBytes(data.toString());

                dos.close();

        }

        /**

         * Obtiene la extensi?n del fichero de georreferenciaci?n

         * @return String con la extensi?n del fichero de georreferenciaci?n dependiendo

         * del valor del formato obtenido del servidor. Por defecto asignaremos un .wld

         */

        private String getExtensionWorldFile(String file){

                String ext = file.substring(file.lastIndexOf(".") + 1);

                String extWorldFile = ".wld";

            if(ext.equals("tif") || ext.equals("tiff"))

                    extWorldFile = ".tfw";

            if(ext.equals("jpeg") || ext.equals("jpg"))

                    extWorldFile = ".jpgw";

            return extWorldFile;

        }

        /**

         * Calcula la transformaci?n af?n a partir de los puntos introducidos por m?nimos cuadrados.

         * @param af        Transformaci?n af?n

         * @param src        Puntos de entrada en coordenadas pixel

         * @param dst        Puntos de destino en coordenadas del mundo

         * @param order        Orden del sistema de ecuaciones

         */

        public AffineT computeLeastSquaresAffine(AffineT af, PointT[] src, PointT[] dst, int order){

                double[] b = new double[dst.length];

                int i,cofs;

                af.setOrder(order);

                cofs = order <= 2 ? order*3 : 10;

                af.mallocCofs(cofs);

                //First compute the X cofs

                for(i = 0; i < dst.length; i++) {

                        b[i] = dst[i].getX();

                }

                af.setCofX(singleLeastSquaresAffine(af.getXcofs(), src, b, order));

                //Now compute the Y cofs

                for(i = 0; i < dst.length; i++) {

                        b[i] = dst[i].getY();

                }

                af.setCofY(singleLeastSquaresAffine(af.getYcofs(), src, b, order));

                return af;

        }

        private double[] singleLeastSquaresAffine(double[] a, PointT[] src, double[] dst, int order){

                int i,j;

                int n = dst.length;

                int points = order <= 2 ? order * 3 : 10;

                double[][] combined = new double[points][points + 1];

                double[][] A = new double[n + 1][points];

                double[][] b = new double[n + 1][1];

                for(i = 0; i < n; i++) {

                        A[i][0] = 1.0;

                        A[i][1] = src[i].getX();

                        A[i][2] = src[i].getY();

                        if(order > 1) {

                                A[i][3] = src[i].getX() * src[i].getX();

                                A[i][4] = src[i].getX() * src[i].getY();

                                A[i][5] = src[i].getY() * src[i].getY();

                        }

                        if(order > 2) {

                                A[i][6] = src[i].getX() * src[i].getX() * src[i].getX();

                                A[i][7] = src[i].getX() * src[i].getX() * src[i].getY();

                                A[i][8] = src[i].getX() * src[i].getY() * src[i].getY();

                                A[i][9] = src[i].getY() * src[i].getY() * src[i].getY();

                        }

                        b[i][0] = dst[i];

                }

                double[][] Atrans = MathUtils.transpose(A, n, points);

                double[][] left = MathUtils.multmatrix(Atrans,A,points,n,points);

                double[][] right = MathUtils.multmatrix(Atrans,b,points,n,1);

                for(i = 0; i < points; i++) {

                        combined[i][0] = right[i][0];

                        for(j = 0; j < points; j++) {

                                combined[i][j+1] = left[i][j];

                        }

                }

                try{

                        combined = solve(combined, points, points+1);

                }catch(DataFormatException ex){

                        System.err.println("Can't solve matrix");

                }

                for(i = 0; i < points; i++) {

                        a[i] = combined[i][0];

                }

                return a;

        }

        private double[][] solve(double[][] mat,int rows,int cols)throws DataFormatException{

                int i,j,k;

                double d,tmp;

                int big;

                for(i = 0; i < rows; i++) {

                        // Find largest row

                                big = i;

                                for(j = i; j < rows; j++) {

                                        if(Math.abs(mat[j][i+1]) > Math.abs(mat[big][i+1])) {

                                                big = j;

                                        }

                                }

                                // swap row i and row big

                                for(k = 0; k < cols ; k++) {

                                        tmp = mat[i][k];

                                        mat[i][k] = mat[big][k];

                                        mat[big][k] = tmp;

                                }

                        if(mat[i][i+1] == 0)

                                throw new DataFormatException();

                        d = 1.0 / mat[i][i+1];

                        for(j = 0; j < cols ; j++) {

                                mat[i][j] *= d;

                                //assert(!isnan(mat[i][j]));

                        }

                        for(k = 0; k < rows; k++) {

                                if(k == i)

                                        continue;

                                if(mat[k][i+1] != 0.0) {

                                d = mat[k][i+1] / mat[i][i+1];

                                        for(j = 0; j < cols ; j++) {

                                                mat[k][j] -= d*mat[i][j];

                                                //assert(!isnan(mat[k][j]));

                                        }

                                }

                        }

                }

                return mat;

        }

        /**

         * Obtiene las coordenadas de un punto en coordenadas del mundo a partir de una transformaci?n

         * y el punto de entrada en coordenadas de la imagen.

         * @param sx        Coordenada X del punto de entrada

         * @param syz        Coordenada Y del punto de entrada

         * @param af        Transformaci?n

         * @return                Punto transformado

         */

        public double[] transformPoint(int sx, int sy, AffineT af){

                double[] p = new double[2];

                if(af.getOrder() == 1) {

                        p[0] = af.getCofX(0) + af.getCofX(1) * sx + af.getCofX(2) * sy;

                        p[1] = af.getCofY(0) + af.getCofY(1) * sx + af.getCofY(2) * sy;

                }

                else if(af.getOrder() == 2) {

                        p = quadTransformPoint(sx, sy, af);

                }

                else {

                        p = cubicTransformPoint(sx, sy, af);

                 }

                return p;

        }

        private static double[] quadTransformPoint(int sx, int sy, AffineT af){

                double[] p = new double[2];

                p[0] = af.getCofX(0) + af.getCofX(1) * sx + af.getCofX(2) * sy +

                        af.getCofX(3) * sx * sx + af.getCofX(4) * sx * sy + af.getCofX(5) * sy * sy;

                p[1] = af.getCofY(0) + af.getCofY(1) * sx + af.getCofY(2) * sy +

                        af.getCofY(3) * sx * sx + af.getCofY(4) * sx * sy + af.getCofY(5) * sy * sy;

                return p;

        }

        private static double[] cubicTransformPoint(int sx, int sy, AffineT af){

                double[] p = new double[2];

                p[0] = af.getCofX(0) + af.getCofX(1) * sx + af.getCofX(2) * sy +

                        af.getCofX(3) * sx * sx + af.getCofX(4) * sx * sy + af.getCofX(5) * sy * sy +

                        af.getCofX(6) * sx * sx * sx + af.getCofX(7) * sx * sx * sy +

                        af.getCofX(8)*sx*sy*sy + af.getCofX(9)*sy*sy*sy;

                p[1] = af.getCofY(0) + af.getCofY(1) * sx + af.getCofY(2) * sy +

                        af.getCofY(3) * sx * sx + af.getCofY(4) * sx * sy + af.getCofY(5) * sy * sy +

                        af.getCofY(6) * sx * sx * sx + af.getCofY(7) * sx * sx * sy +

                        af.getCofY(8) * sx * sy * sy + af.getCofY(9) * sy * sy * sy;

                return p;

        }

}