| gvSIG desktop 1

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

 *

 * Copyright (C) 2006 IVER T.I. 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.

 */

package org.gvsig.georeferencing;

import java.awt.geom.AffineTransform;

import java.awt.geom.Point2D;

import java.io.BufferedOutputStream;

import java.io.BufferedReader;

import java.io.BufferedWriter;

import java.io.DataOutputStream;

import java.io.File;

import java.io.FileInputStream;

import java.io.FileNotFoundException;

import java.io.FileOutputStream;

import java.io.IOException;

import java.io.InputStreamReader;

import java.io.OutputStreamWriter;

import java.util.zip.DataFormatException;

import org.cresques.cts.ICoordTrans;

import org.gvsig.georeferencing.utils.AffineT;

import org.gvsig.georeferencing.utils.GeoUtils;

import org.gvsig.georeferencing.utils.PointT;

import org.gvsig.raster.dataset.BandList;

import org.gvsig.raster.dataset.FileNotOpenException;

import org.gvsig.raster.dataset.GeoInfo;

import org.gvsig.raster.dataset.IBuffer;

import org.gvsig.raster.dataset.InvalidSetViewException;

import org.gvsig.raster.dataset.RasterDataset;

import org.gvsig.raster.dataset.RasterDriverException;

import org.gvsig.raster.dataset.RasterMetaFileTags;

import org.gvsig.raster.dataset.io.rmf.RmfBlocksManager;

import org.gvsig.raster.dataset.serializer.GeoInfoRmfSerializer;

import org.gvsig.raster.datastruct.Extent;

import org.xmlpull.v1.XmlSerializer;

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

/** 

 * Operaciones necesarias para la georreferenciaci?n a partir de la capa de puntos. 

 * En base a unos puntos de control de origen y otros de destino se crea una matriz

 * de transformaci?n para asignar la nueva posici?n a la imagen y crear ficheros 

 * de georreferenciaci?n en dos formatos: worldfile y rasterMetaFile

 * 

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

 */

public class GeoOperations{

        //**********************Vars**********************************

        private int                         order = 1;

        private PointT[]                 src = null;

        private PointT[]                 dst = null;

        private AffineT                affine = null;

        private boolean         createWorldFile = false;

        //**********************End Vars******************************

        //**********************Methods*******************************

        /**

         * Constructor

         */

        public GeoOperations(){}

        /**

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

         * @param lyr

         */

        public GeoOperations(FLyrPoints lyr){

                FLyrPoints         lyrPoints = lyr;

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

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

                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.getCountActivePoints() >= 3 * 10)

                        order = 3;

                else if(lyr.getCountActivePoints() >= 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 createGeorefFile(int widthPx, int heightPx, String file){

                try{

                        File f = new File(file);

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

                        if(createWorldFile)

                                createWorldFile(affine, widthPx, heightPx, nameWorldFile);

                        createRasterMetaFile(affine, widthPx, heightPx, nameWorldFile.substring(0, nameWorldFile.lastIndexOf("."))+".rmf");

                }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[] begin = transformPoint(0, 0, affine);

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

                //double pixelSizeX = (end[0] - begin[0])/(widthPx - 1);

                //double pixelSizeY = (end[1] - begin[1])/(heightPx - 1);

                File f = new File(nameWordlFile);

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

            data.append(affine.getCofX(1)+"\n");//pixelSizeX+"\n");

            data.append(affine.getCofX(2)+"\n");

            data.append(affine.getCofY(1)+"\n");

            data.append(affine.getCofY(2)+"\n");//(-1  * pixelSizeY)+"\n");

            data.append(affine.getCofX(0)+"\n");//""+begin[0]+"\n");

            data.append(affine.getCofY(0)+"\n");//""+end[1]+"\n");

            dos.writeBytes(data.toString());

                dos.close();

        }

        /**

         * Si no existe crea un fichero .rmf con la georreferenciaci?n de la imagen. Si existe este 

         * deber? a?adir o modificar la informaci?n de georreferenciaci?n en el fichero.

         * 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 createRasterMetaFile(AffineT affine, int widthPx, int heightPx, String nameWorldFile) throws IOException {

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

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

                RmfBlocksManager manager = new RmfBlocksManager(nameWorldFile);

                DatasetSer dataset = new DatasetSer();

                dataset.setExtent(new Extent(min[0], max[1], max[0], min[1]));

                AffineTransform at = new AffineTransform(affine.getCofX(1), affine.getCofX(2), affine.getCofY(1), affine.getCofY(2), affine.getCofX(0), affine.getCofY(0));

                dataset.setAffineTransform(at);

                dataset.w = widthPx;

                dataset.h = heightPx;

                GeoInfoRmfSerializer ser = new GeoInfoRmfSerializer(dataset);

                manager.addClient(ser);

                try {

                        manager.write();

                } catch (FileNotFoundException e) {

                        e.printStackTrace();

                } catch (IOException e) {

                        e.printStackTrace();

                }

        }

        /**

         * Calcula la transformaci?n af?n a partir de los puntos introducidos por el m?todo de 

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

                if(src.length == 2){

                        af.setCofX(0, 2);

                        af.setCofY(0, 1);

                        double distWCX = dst[1].getX() - dst[0].getX();

                        double distWCY = dst[1].getY() - dst[0].getY();

                        double distPxX = src[1].getX() - src[0].getX();

                        double distPxY = src[1].getY() - src[0].getY();

                        double psX = distWCX / distPxX;

                        double psY = distWCY / distPxY;

                        af.setCofX(psX, 1);

                        af.setCofY(psY, 2);

                        double initWCX = dst[0].getX() - ((src[0].getX() + 1) * psX);

                        double initWCY = dst[0].getY() - ((src[0].getY() + 1) * psY);

                        af.setCofX(initWCX, 0);

                        af.setCofY(initWCY, 0);

                        return af;

                }

                //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 = GeoUtils.transpose(A, n, points);

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

                double[][] right = GeoUtils.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;

        }        

        //**********************End Methods***************************

        //**********************Setters & Getters*********************

        /**

         * M?todo para notificar si se desea crear o no el fichero de coordenadas .tfw, .wld .jpgw ,... 

         * @param createWorldFile

         */

        public void setCreateWorldFile(boolean createWorldFile) {

                this.createWorldFile = createWorldFile;

        }

        /**

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

                String extWorldFile = ".wld";

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

                    extWorldFile = ".tfw";

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

                    extWorldFile = ".jgw";

            if(ext.equals("gif"))

                    extWorldFile = ".gfw";

            if(ext.equals("png"))

                    extWorldFile = ".pgw";

            return extWorldFile;

        }*/

        /**

         * Obtiene la matriz de transformaci?n

         * @return AffineT

         */

        public AffineT getAffine() {

                return affine;

        }

        //**********************End Setters & Getters*****************

        public class DatasetSer extends RasterDataset {

                public int w, h;

                AffineTransform at = new AffineTransform();

                public void close() {}

                public int getBlockSize() {return 0;}

                public Object getData(int arg0, int arg1, int arg2) throws InvalidSetViewException, FileNotOpenException, RasterDriverException {return null;}

                public int getHeight() {return h;}

                public Extent getView() {return null;}

                public int getWidth() {return w;}

                public IBuffer getWindowRaster(int arg0, int arg1, int arg2, int arg3, BandList arg4, IBuffer arg5) {return null;}

                public IBuffer getWindowRaster(double arg0, double arg1, double arg2, double arg3, BandList arg4, IBuffer arg5, boolean arg6) {return null;}

                public IBuffer getWindowRaster(int arg0, int arg1, int arg2, int arg3, int arg4, int arg5, BandList arg6, IBuffer arg7) {return null;}

                public IBuffer getWindowRaster(double arg0, double arg1, double arg2, double arg3, int arg4, int arg5, BandList arg6, IBuffer arg7, boolean arg8) {return null;}

                public GeoInfo load() {return null;}

                public Point2D rasterToWorld(Point2D arg0) {return null;}

                public void reProject(ICoordTrans arg0) {}

                public Object readBlock(int arg0, int arg1) throws InvalidSetViewException, FileNotOpenException, RasterDriverException {return null;}

                public Object readCompleteLine(int arg0, int arg1) throws InvalidSetViewException, FileNotOpenException, RasterDriverException {return null;}

                public void setExtentTransform(double arg0, double arg1, double arg2, double arg3) {}

                public void setView(Extent arg0) {}

                public Point2D worldToRaster(Point2D arg0) { return null;}

                public void setAffineTransform(AffineTransform at) {this.at = at;}

                public AffineTransform getAffineTransform() { return at;}

        }

}