org.csa.rstb.polarimetric.gpf.OrientationAngleCorrectionOp.java Source code

Java tutorial

Introduction

Here is the source code for org.csa.rstb.polarimetric.gpf.OrientationAngleCorrectionOp.java

Source

/*
 * Copyright (C) 2015 by Array Systems Computing Inc. http://www.array.ca
 *
 * 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 3 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, see http://www.gnu.org/licenses/
 */
package org.csa.rstb.polarimetric.gpf;

import com.bc.ceres.core.ProgressMonitor;
import org.apache.commons.math3.util.FastMath;
import org.esa.s1tbx.io.PolBandUtils;
import org.esa.snap.core.datamodel.Band;
import org.esa.snap.core.datamodel.Product;
import org.esa.snap.core.datamodel.ProductData;
import org.esa.snap.core.gpf.Operator;
import org.esa.snap.core.gpf.OperatorException;
import org.esa.snap.core.gpf.OperatorSpi;
import org.esa.snap.core.gpf.Tile;
import org.esa.snap.core.gpf.annotations.OperatorMetadata;
import org.esa.snap.core.gpf.annotations.SourceProduct;
import org.esa.snap.core.gpf.annotations.TargetProduct;
import org.esa.snap.core.util.ProductUtils;
import org.esa.snap.engine_utilities.datamodel.AbstractMetadata;
import org.esa.snap.engine_utilities.eo.Constants;
import org.esa.snap.engine_utilities.gpf.InputProductValidator;
import org.esa.snap.engine_utilities.gpf.OperatorUtils;
import org.esa.snap.engine_utilities.gpf.TileIndex;

import java.awt.*;
import java.util.ArrayList;
import java.util.List;
import java.util.Map;

/**
 * Perform polarization orientation angle correction for given coherency matrix
 */

@OperatorMetadata(alias = "Orientation-Angle-Correction", category = "Radar/Polarimetric", authors = "Jun Lu, Luis Veci", version = "1.0", copyright = "Copyright (C) 2014 by Array Systems Computing Inc.", description = "Perform polarization orientation angle correction for given coherency matrix")
public final class OrientationAngleCorrectionOp extends Operator {

    @SourceProduct(alias = "source")
    private Product sourceProduct;
    @TargetProduct
    private Product targetProduct;

    private PolBandUtils.PolSourceBand[] srcBandList;
    private PolBandUtils.MATRIX sourceProductType;

    private final static double PI4 = Constants.PI / 4.0;
    private final static double PI2 = Constants.PI / 2.0;

    /**
     * Initializes this operator and sets the one and only target product.
     * <p>The target product can be either defined by a field of type {@link Product} annotated with the
     * {@link TargetProduct TargetProduct} annotation or
     * by calling {@link #setTargetProduct} method.</p>
     * <p>The framework calls this method after it has created this operator.
     * Any client code that must be performed before computation of tile data
     * should be placed here.</p>
     *
     * @throws OperatorException If an error occurs during operator initialisation.
     * @see #getTargetProduct()
     */
    @Override
    public void initialize() throws OperatorException {

        try {
            final InputProductValidator validator = new InputProductValidator(sourceProduct);
            validator.checkIfSARProduct();
            validator.checkIfSLC();
            validator.checkIfTOPSARBurstProduct(false);

            sourceProductType = PolBandUtils.getSourceProductType(sourceProduct);

            checkSourceProductType();

            srcBandList = PolBandUtils.getSourceBands(sourceProduct, sourceProductType);

            createTargetProduct();

        } catch (Throwable e) {
            OperatorUtils.catchOperatorException(getId(), e);
        }
    }

    /**
     * Create target product.
     */
    private void createTargetProduct() {

        targetProduct = new Product(sourceProduct.getName(), sourceProduct.getProductType(),
                sourceProduct.getSceneRasterWidth(), sourceProduct.getSceneRasterHeight());

        addSelectedBands();

        ProductUtils.copyProductNodes(sourceProduct, targetProduct);

        AbstractMetadata.getAbstractedMetadata(targetProduct).setAttributeInt(AbstractMetadata.polsarData, 1);
    }

    private void checkSourceProductType() {

        if (sourceProductType != PolBandUtils.MATRIX.FULL && sourceProductType != PolBandUtils.MATRIX.T3
                && sourceProductType != PolBandUtils.MATRIX.T4 && sourceProductType != PolBandUtils.MATRIX.C3
                && sourceProductType != PolBandUtils.MATRIX.C4) {

            throw new OperatorException("OA correction cannot be performed on the input product type");
        }
    }

    /**
     * Add bands to the target product.
     *
     * @throws OperatorException The exception.
     */
    private void addSelectedBands() throws OperatorException {

        final List<String> targetBandNameList = new ArrayList<>(10);
        targetBandNameList.add("T11");
        targetBandNameList.add("T12_real");
        targetBandNameList.add("T12_imag");
        targetBandNameList.add("T13_real");
        targetBandNameList.add("T13_imag");
        targetBandNameList.add("T22");
        targetBandNameList.add("T23_real");
        targetBandNameList.add("T23_imag");
        targetBandNameList.add("T33");
        //targetBandNameList.add("Ori_Ang");

        final String[] bandNames = targetBandNameList.toArray(new String[targetBandNameList.size()]);
        for (PolBandUtils.PolSourceBand bandList : srcBandList) {
            final Band[] targetBands = OperatorUtils.addBands(targetProduct, bandNames, bandList.suffix);
            bandList.addTargetBands(targetBands);
        }
    }

    /**
     * Called by the framework in order to compute the stack of tiles for the given target bands.
     * <p>The default implementation throws a runtime exception with the message "not implemented".</p>
     *
     * @param targetTiles     The current tiles to be computed for each target band.
     * @param targetRectangle The area in pixel coordinates to be computed (same for all rasters in <code>targetRasters</code>).
     * @param pm              A progress monitor which should be used to determine computation cancelation requests.
     * @throws OperatorException if an error occurs during computation of the target rasters.
     */
    @Override
    public void computeTileStack(Map<Band, Tile> targetTiles, Rectangle targetRectangle, ProgressMonitor pm)
            throws OperatorException {
        try {

            final int x0 = targetRectangle.x;
            final int y0 = targetRectangle.y;
            final int w = targetRectangle.width;
            final int h = targetRectangle.height;
            final int maxY = y0 + h;
            final int maxX = x0 + w;
            //System.out.println("x0 = " + x0 + ", y0 = " + y0 + ", w = " + w + ", h = " + h);

            final double[][] C3r = new double[3][3];
            final double[][] C3i = new double[3][3];
            final double[][] C4r = new double[4][4];
            final double[][] C4i = new double[4][4];
            final double[][] T3r = new double[3][3];
            final double[][] T3i = new double[3][3];
            final double[][] T4r = new double[4][4];
            final double[][] T4i = new double[4][4];

            final TileIndex tgtIndex = new TileIndex(targetTiles.get(getTargetProduct().getBandAt(0)));
            for (final PolBandUtils.PolSourceBand bandList : srcBandList) {

                final Tile[] sourceTiles = new Tile[bandList.srcBands.length];
                final ProductData[] dataBuffers = new ProductData[bandList.srcBands.length];
                for (int i = 0; i < bandList.srcBands.length; i++) {
                    sourceTiles[i] = getSourceTile(bandList.srcBands[i], targetRectangle);
                    dataBuffers[i] = sourceTiles[i].getDataBuffer();
                }
                final TileIndex srcIndex = new TileIndex(sourceTiles[0]);

                final ProductData[] targetDataBuffers = new ProductData[9];
                for (final Band targetBand : bandList.targetBands) {
                    final String targetBandName = targetBand.getName();
                    final ProductData dataBuffer = targetTiles.get(targetBand).getDataBuffer();
                    if (PolBandUtils.isBandForMatrixElement(targetBandName, "11"))
                        targetDataBuffers[0] = dataBuffer;
                    else if (PolBandUtils.isBandForMatrixElement(targetBandName, "12_real"))
                        targetDataBuffers[1] = dataBuffer;
                    else if (PolBandUtils.isBandForMatrixElement(targetBandName, "12_imag"))
                        targetDataBuffers[2] = dataBuffer;
                    else if (PolBandUtils.isBandForMatrixElement(targetBandName, "13_real"))
                        targetDataBuffers[3] = dataBuffer;
                    else if (PolBandUtils.isBandForMatrixElement(targetBandName, "13_imag"))
                        targetDataBuffers[4] = dataBuffer;
                    else if (PolBandUtils.isBandForMatrixElement(targetBandName, "22"))
                        targetDataBuffers[5] = dataBuffer;
                    else if (PolBandUtils.isBandForMatrixElement(targetBandName, "23_real"))
                        targetDataBuffers[6] = dataBuffer;
                    else if (PolBandUtils.isBandForMatrixElement(targetBandName, "23_imag"))
                        targetDataBuffers[7] = dataBuffer;
                    else if (PolBandUtils.isBandForMatrixElement(targetBandName, "33"))
                        targetDataBuffers[8] = dataBuffer;
                }

                final double[][] Tr = new double[3][3];
                final double[][] Ti = new double[3][3];

                int srcIdx, tgtIdx;
                double theta, c, s, c2, s2, cs;
                for (int y = y0; y < maxY; ++y) {
                    srcIndex.calculateStride(y);
                    tgtIndex.calculateStride(y);

                    for (int x = x0; x < maxX; ++x) {
                        srcIdx = srcIndex.getIndex(x);
                        tgtIdx = tgtIndex.getIndex(x);

                        if (sourceProductType == PolBandUtils.MATRIX.FULL) {

                            PolOpUtils.getT3(srcIdx, sourceProductType, dataBuffers, T3r, T3i);

                        } else if (sourceProductType == PolBandUtils.MATRIX.T3) {

                            PolOpUtils.getCoherencyMatrixT3(srcIdx, dataBuffers, T3r, T3i);

                        } else if (sourceProductType == PolBandUtils.MATRIX.T4) {

                            PolOpUtils.getCoherencyMatrixT4(srcIdx, dataBuffers, T4r, T4i);
                            PolOpUtils.t4ToT3(T4r, T4i, T3r, T3i);

                        } else if (sourceProductType == PolBandUtils.MATRIX.C3) {

                            PolOpUtils.getCovarianceMatrixC3(srcIdx, dataBuffers, C3r, C3i);
                            PolOpUtils.c3ToT3(C3r, C3i, T3r, T3i);

                        } else if (sourceProductType == PolBandUtils.MATRIX.C4) {

                            PolOpUtils.getCovarianceMatrixC4(srcIdx, dataBuffers, C4r, C4i);
                            PolOpUtils.c4ToT4(C4r, C4i, T4r, T4i);
                            PolOpUtils.t4ToT3(T4r, T4i, T3r, T3i);
                        }

                        theta = estimateOrientationAngle(T3r[1][2], T3r[1][1], T3r[2][2]);
                        c = FastMath.cos(2 * theta);
                        s = FastMath.sin(2 * theta);
                        c2 = c * c;
                        s2 = s * s;
                        cs = c * s;

                        Tr[0][0] = T3r[0][0];
                        Tr[0][1] = T3r[0][1] * c - T3r[0][2] * s;
                        Ti[0][1] = T3i[0][1] * c - T3i[0][2] * s;
                        Tr[0][2] = T3r[0][1] * s + T3r[0][2] * c;
                        Ti[0][2] = T3i[0][1] * s + T3i[0][2] * c;
                        Tr[1][1] = T3r[1][1] * c2 + T3r[2][2] * s2 - 2 * T3r[1][2] * cs;
                        Tr[1][2] = T3r[1][2] * (c2 - s2) + (T3r[1][1] - T3r[2][2]) * cs;
                        Ti[1][2] = T3i[1][2];
                        Tr[2][2] = T3r[1][1] * s2 + T3r[2][2] * c2 + 2 * T3r[1][2] * cs;

                        saveT3(Tr, Ti, tgtIdx, targetDataBuffers);
                    }
                }
            }
        } catch (Throwable e) {
            OperatorUtils.catchOperatorException(getId(), e);
        } finally {
            pm.done();
        }
    }

    /**
     * Compute polarization orientation angle.
     *
     * @param t23Re Real part of the t23 element of coherency matrix
     * @param t22   The t22 element of the coherency matrix
     * @param t33   The t33 element of the coherency matrix
     * @return The polarization orientation angle in radian
     */
    private static double estimateOrientationAngle(final double t23Re, final double t22, final double t33) {
        if (t33 == 0.0) {
            return 0.0;
        }

        double theta = 0.25 * (Math.atan2(2 * t23Re, t33 - t22) + Constants.PI);
        if (theta > PI4) {
            theta -= PI2;
        }
        return theta;
    }

    private static void saveT3(final double[][] Tr, final double[][] Ti, final int idx,
            final ProductData[] targetDataBuffers) {

        targetDataBuffers[0].setElemFloatAt(idx, (float) Tr[0][0]); // T11
        targetDataBuffers[1].setElemFloatAt(idx, (float) Tr[0][1]); // T12_real
        targetDataBuffers[2].setElemFloatAt(idx, (float) Ti[0][1]); // T12_imag
        targetDataBuffers[3].setElemFloatAt(idx, (float) Tr[0][2]); // T13_real
        targetDataBuffers[4].setElemFloatAt(idx, (float) Ti[0][2]); // T13_imag
        targetDataBuffers[5].setElemFloatAt(idx, (float) Tr[1][1]); // T22
        targetDataBuffers[6].setElemFloatAt(idx, (float) Tr[1][2]); // T23_real
        targetDataBuffers[7].setElemFloatAt(idx, (float) Ti[1][2]); // T23_imag
        targetDataBuffers[8].setElemFloatAt(idx, (float) Tr[2][2]); // T33
    }

    /**
     * The SPI is used to register this operator in the graph processing framework
     * via the SPI configuration file
     * {@code META-INF/services/org.esa.snap.core.gpf.OperatorSpi}.
     * This class may also serve as a factory for new operator instances.
     *
     * @see OperatorSpi#createOperator()
     * @see OperatorSpi#createOperator(java.util.Map, java.util.Map)
     */
    public static class Spi extends OperatorSpi {
        public Spi() {
            super(OrientationAngleCorrectionOp.class);
        }
    }
}