Application: gvSIG desktop » gvSIG raster

/* gvSIG. Geographic Information System of the Valencian Government

 *

 * Copyright (C) 2007-2008 Infrastructures and Transports Department

 * of the Valencian Government (CIT)

 *

 * 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., 51 Franklin Street, Fifth Floor, Boston,

 * MA  02110-1301, USA.

 *

 */

package org.gvsig.raster.principalcomponents.algorithm;

import java.awt.geom.Rectangle2D;

import java.util.List;

import org.cresques.Messages;

import org.gvsig.fmap.dal.coverage.RasterLocator;

import org.gvsig.fmap.dal.coverage.dataset.Buffer;

import org.gvsig.fmap.dal.coverage.datastruct.Extent;

import org.gvsig.fmap.dal.coverage.datastruct.NoData;

import org.gvsig.fmap.dal.coverage.exception.CloneException;

import org.gvsig.fmap.dal.coverage.exception.ProcessInterruptedException;

import org.gvsig.fmap.dal.coverage.exception.ROIException;

import org.gvsig.fmap.dal.coverage.store.RasterDataStore;

import org.gvsig.fmap.dal.exception.CloseException;

import org.gvsig.raster.algorithm.process.DataProcess;

import org.gvsig.raster.algorithm.process.ProcessException;

import org.gvsig.raster.roi.ROI;

import org.slf4j.Logger;

import org.slf4j.LoggerFactory;

import Jama.EigenvalueDecomposition;

import Jama.Matrix;

/**

 * This Process computes a <code>PCStatsDataStructure</code> Object from bands of two images in a study zone.

 * <code>PCStatsDataStructure</code> is a object that contains the needed statistics to Principal Components decomposition

 * of images.

 * 

 * @author Nacho Brodin (nachobrodin@gmail.com)

 * 

 */

public class PCAStatisticsProcess extends DataProcess {

    public static final String          LAYER_ZONA_ESTUDIO                         = "LAYER_ZONA_ESTUDIO";

    public static String         BANDS                                                         = "BANDS";

    public static String         PATH                                                         = "PATH";

    public static String         ROI_EPSG                       = "ROI_EPSG";

    public static String         WINDOW                         = "WINDOW";

        public static String         RASTER_STORE                   = "RASTER_STORE";

        public static String         STATS_RESULT                   = "STATS_RESULT";

    private boolean[]                          bands                                                         = null;

    private RasterDataStore      store                                                         = null;

    private Buffer               buffer                                                 = null;

    private Buffer               roiBuffer                                                 = null;

    private boolean              roiBufferCalculated            = false;

        private double                                 autovalors[]                                        = null;

        private Matrix                                  coVarMatrix                                        = null;

        private Matrix                                  autoVectorMatrix                                = null;

        private double                                  coVar[][]                                                 = null;

        private static Logger                  logger                                                       = LoggerFactory.getLogger(PCAStatisticsProcess.class.getName());

    private double[][]                          maxminmean                                                = null;

    private        Extent                                  extentResult                                        = null; 

    private String               path                           = null;

    private List<ROI>            rois                           = null;

    private boolean[]            statsByBandCalculated          = null;

        public static void registerParameters() {

                registerInputParameter(RASTER_STORE, RasterDataStore.class, PrincipalComponentsAlgorithmLibrary.PC_STATS_PROCESS_LABEL);

                registerInputParameter(BANDS, Boolean[].class, PrincipalComponentsAlgorithmLibrary.PC_STATS_PROCESS_LABEL);

                registerInputParameter(PATH, String.class, PrincipalComponentsAlgorithmLibrary.PC_STATS_PROCESS_LABEL);

                registerOutputParameter(STATS_RESULT, PCStatsDataStructure.class, PrincipalComponentsAlgorithmLibrary.PC_STATS_PROCESS_LABEL);

                registerOutputParameter(PATH, String.class, PrincipalComponentsAlgorithmLibrary.PC_STATS_PROCESS_LABEL);

                registerOutputParameter(BANDS, Boolean[].class, PrincipalComponentsAlgorithmLibrary.PC_STATS_PROCESS_LABEL);

                registerOutputParameter(ROI_EPSG, String.class, PrincipalComponentsAlgorithmLibrary.PC_STATS_PROCESS_LABEL);

        }

        @Override

        public void init() {

                store = (RasterDataStore)getParam(RASTER_STORE);

                bands = (boolean[])getParam(BANDS);

                path = getStringParam(PATH);

        }

        public double[] getAutoValors() {

                return  autovalors;

        }

        /**

         * @return Autovectors Matrix

         */

        public Matrix getAutoVectorMatrix() {

                return autoVectorMatrix;

        }

        /**

         * @return Varianze-covarianze Matrix

         */

        public Matrix getcoVarMatrix() {

                return coVarMatrix;

        }

        public boolean isInsideOfROI(int x, int y, List<ROI> rois, Extent extentResult) {

                if(roiBufferCalculated)

                        return (roiBuffer.getElemByte(y, x, 0) == 1);

                if(roiBuffer == null)

                        roiBuffer = createOutputBuffer(buffer.getWidth(), buffer.getHeight(), 1, Buffer.TYPE_BYTE);

                boolean a = super.isInsideOfROI(x, y, rois, extentResult);

                roiBuffer.setElem(y, x, 0, a ? (byte)1: (byte)0);

                return a;

        }

        /**

         * Calculate the Varianze-covarianze Matrix from bands of a Grid

         * @throws InterruptedException 

         * @throws ProcessInterruptedException 

         */

        private double[][] covarianceOptimize() throws InterruptedException, ProcessInterruptedException {

        //NoData

        NoData nd = RasterLocator.getManager().getDataStructFactory().createDefaultNoData(1, Buffer.TYPE_DOUBLE);

        double nodata = nd.getValue().doubleValue();

                double dSum = 0;

                int iValues = 0;

                double[][] coV = new double[buffer.getBandCount()][buffer.getBandCount()];

                double valorBandai = 0, valorBandaj = 0;

                maxminmean = new double[buffer.getBandCount()][4];

                statsByBandCalculated = new boolean[buffer.getBandCount()];

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

                        statsByBandCalculated[i] = false;

                }

                loadDataStoreStatistics();

                int percent = 0;

                int maxCounter = 0;

                for (int i = buffer.getBandCount(); i >= 0; i--) {

                        maxCounter += i;

                }

                for (int i = 0; i < buffer.getBandCount(); i++) {

                        double meanBufferToOperatei = getMeanValue(i, buffer);

                        for (int j = i; j < buffer.getBandCount(); j++) {

                                iValues = 0;

                                dSum = 0;

                                double meanBufferToOperatej = getMeanValue(j, buffer);

                                for (int row = 0; row < buffer.getHeight(); row++) {

                                        for (int col = 0; col < buffer.getWidth(); col++) {

                                                if(isInsideOfROI(col, row, rois, extentResult)) {

                                                        valorBandai = getData(buffer, row, col, i) - meanBufferToOperatei;

                                                        valorBandaj = getData(buffer, row, col, j) - meanBufferToOperatej;

                                                        dSum += valorBandai * valorBandaj;

                                                        iValues++;

                                                }

                                        }

                                        updatePercent(percent, maxCounter);

                                }

                                roiBufferCalculated = true;

                                // Asigno el valor a la matriz 

                                if (iValues > 1)

                                        coV[i][j] = dSum / (double)(iValues);        

                                else

                                        coV[i][j] = nodata;

                                percent ++;

                                updatePercent(percent, maxCounter);

                        }

                }

                for (int i = 0; i < buffer.getBandCount(); i++) {

                        for (int j = 0; j < buffer.getBandCount(); j++) {

                                if(j < i)

                                        coV[i][j] = coV[j][i];

                        }

                }

                return coV;

        }

//         * Devuelve la media de los valores de los pixeles de un buffer incluidos en la zona de estudio

        /**

         * Get the mean of values of pixels from a buffer in a study zone

         * @param bufferToOperate

         * @return

         */

        private double getMeanValue(int band, Buffer bufferToOperate) {

                if(statsByBandCalculated[band])

                        return maxminmean[band][2];

                double mean = 0;

                double variance = 0;

                double min = Double.MAX_VALUE;

                double max = Double.MIN_VALUE;

                int valuesCounter = 0;

                for (int j = 0; j < bufferToOperate.getHeight(); j++) {

                        for (int i = 0; i < bufferToOperate.getWidth(); i++) {

                                if (isInsideOfROI(i, j, rois, extentResult)) {

                                        double value = getData(bufferToOperate, j, i, band);

                                        if(value < min) {

                                                min = value;

                                        }

                                        if(value > max) {

                                                max = value;

                                        }

                                        mean += value;

                                        variance += value * value;

                                        valuesCounter++;

                                }

                        }

                }

                roiBufferCalculated = true;

                mean = (mean / valuesCounter);

                variance = variance / (double) valuesCounter - mean * mean;

                maxminmean[band][0] = max;

                maxminmean[band][1] = min;

                maxminmean[band][2] = mean;

                maxminmean[band][3] = variance;

                statsByBandCalculated[band] = true;

                return mean;

        }

        private void loadDataStoreStatistics() {

                if(isAnalizedEntireLayer(getOutputWindow(), rois, store) && store.getStatistics().isCalculated()) {

                        double[] max = store.getStatistics().getMax();

                        double[] min = store.getStatistics().getMin();

                        double[] mean = store.getStatistics().getMean();

                        double[] variance = store.getStatistics().getVariance();

                        int iBand = 0;

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

                                if(bands[i]) {

                                        maxminmean[iBand][0] = max[i];

                                        maxminmean[iBand][1] = min[i];

                                        maxminmean[iBand][2] = mean[i];

                                        maxminmean[iBand][3] = variance[i];

                                        statsByBandCalculated[iBand] = true;

                                        iBand ++;

                                }

                        }

                }

        }

        @Override

        public void process() throws ProcessInterruptedException, ProcessException {

                insertLineLog(Messages.getText("preparing_buffers"));

                if(getROIEPSG() != null) {

                        try {

                                rois = store.getRois(getROIEPSG());

                        } catch (ROIException e2) {

                                logger.error(Messages.getText("error_getting_rois"), e2);

                        }

                }

                try {

                        store = store.cloneDataStore();

                } catch (CloneException e1) {

                        throw new ProcessException("Error cloning the sources", e1);

                }

                //Calculo de datos de entrada

                extentResult = getExtentResult(getOutputWindow(), rois, store);

                Rectangle2D sourcePxBBox = getSourcePxBox(extentResult, store);

                buffer = createSourceBuffer(store, sourcePxBBox, bands);

                insertLineLog(Messages.getText("calc_stats"));

                try {

                        coVar = covarianceOptimize();

                } catch (InterruptedException e) {

                        logger.error(Messages.getText("error_covarianza_optimize"), e);

                }

                // Calculo de autovectores:

                coVarMatrix = new Matrix(coVar);

                EigenvalueDecomposition eigenvalueDecomp = new EigenvalueDecomposition(coVarMatrix); 

                //Autovectores

                autoVectorMatrix = eigenvalueDecomp.getV();

                // Autovalores

                autovalors = eigenvalueDecomp.getRealEigenvalues();        

                PCStatsDataStructure resultStatistics = new PCStatsDataStructure(autoVectorMatrix, autovalors, coVarMatrix);

                resultStatistics.setMaxMinMean(maxminmean);

                roiBuffer.dispose();

                //Object[] dataROI = new Object[]{vectROI, (int)shift[0].getX(), (int)shift[0].getY()};        

                addOutputValue(STATS_RESULT, resultStatistics);

                addOutputValue(PATH, path);

                addOutputValue(BANDS, bands);

                addOutputValue(ROI_EPSG, getROIEPSG());

                try {

                        store.close();

                } catch (CloseException e) {

                        throw new ProcessException("Error closing the sources", e);

                }

        }

        public String getTitle() {

                return null;

        }

}