org.esa.s1tbx.sentinel1.gpf.AzimuthShiftOp.java Source code

Java tutorial

Introduction

Here is the source code for org.esa.s1tbx.sentinel1.gpf.AzimuthShiftOp.java

Source

/*
 * Copyright (C) 2014 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.esa.s1tbx.sentinel1.gpf;

import com.bc.ceres.core.ProgressMonitor;
import edu.emory.mathcs.jtransforms.fft.DoubleFFT_1D;
import org.apache.commons.math3.util.FastMath;
import org.esa.s1tbx.insar.gpf.support.Sentinel1Utils;
import org.esa.snap.core.datamodel.Band;
import org.esa.snap.core.datamodel.MetadataAttribute;
import org.esa.snap.core.datamodel.MetadataElement;
import org.esa.snap.core.datamodel.Product;
import org.esa.snap.core.datamodel.ProductData;
import org.esa.snap.core.datamodel.VirtualBand;
import org.esa.snap.core.dataop.downloadable.StatusProgressMonitor;
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.core.util.SystemUtils;
import org.esa.snap.engine_utilities.datamodel.AbstractMetadata;
import org.esa.snap.engine_utilities.gpf.InputProductValidator;
import org.esa.snap.engine_utilities.gpf.OperatorUtils;
import org.esa.snap.engine_utilities.gpf.ReaderUtils;
import org.esa.snap.engine_utilities.gpf.StackUtils;
import org.esa.snap.engine_utilities.gpf.ThreadManager;
import org.esa.snap.engine_utilities.gpf.TileIndex;

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

/**
 * Estimate global azimuth offset using Enhanced Spectral Diversity (ESD) approach.
 * Perform azimuth shift for all bursts in a sub-swath with the azimuth offset above
 * using a frequency domain method.
 */

@OperatorMetadata(alias = "Azimuth-Shift", category = "Radar/Coregistration/S-1 TOPS Coregistration", authors = "Jun Lu, Luis Veci", version = "1.0", copyright = "Copyright (C) 2014 by Array Systems Computing Inc.", description = "Estimate global azimuth offset for the whole image")
public class AzimuthShiftOp extends Operator {

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

    @TargetProduct(description = "The target product which will use the master's grid.")
    private Product targetProduct = null;

    private boolean isOffsetAvailable = false;
    private double azOffset = 0.0;
    private Sentinel1Utils.SubSwathInfo[] subSwath = null;
    private int subSwathIndex = 0;
    private String swathIndexStr = null;
    private String[] subSwathNames = null;
    private String[] polarizations = null;

    private static final String DerampDemodPhase = "derampDemodPhase";

    /**
     * Default constructor. The graph processing framework
     * requires that an operator has a default constructor.
     */
    public AzimuthShiftOp() {
    }

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

            checkDerampDemodPhaseBand();

            final Sentinel1Utils su = new Sentinel1Utils(sourceProduct);
            su.computeDopplerRate();
            subSwath = su.getSubSwath();

            subSwathNames = su.getSubSwathNames();
            if (subSwathNames.length != 1) {
                throw new OperatorException("Split product is expected.");
            } else {
                subSwathIndex = 1;//Integer.parseInt(subSwathNames[0].substring(subSwathNames[0].length()-1));
                swathIndexStr = subSwathNames[0].substring(2);
            }

            polarizations = su.getPolarizations();

            createTargetProduct();

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

    private void checkDerampDemodPhaseBand() {

        boolean hasDerampDemodPhaseBand = false;
        final Band[] sourceBands = sourceProduct.getBands();
        for (Band band : sourceBands) {
            if (band.getName().contains(DerampDemodPhase)) {
                hasDerampDemodPhaseBand = true;
                break;
            }
        }

        if (!hasDerampDemodPhaseBand) {
            throw new OperatorException("Cannot find derampDemodPhase band in source product. "
                    + "Please run Backgeocoding and select \"Output Deramp and Demod Phase\".");
        }
    }

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

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

        ProductUtils.copyProductNodes(sourceProduct, targetProduct);

        final String[] bandNames = sourceProduct.getBandNames();
        for (String srcBandName : bandNames) {
            final Band band = sourceProduct.getBand(srcBandName);
            if (band instanceof VirtualBand) {
                continue;
            }

            Band targetBand;
            if (srcBandName.contains(StackUtils.MST) || srcBandName.contains("derampDemod")) {
                targetBand = ProductUtils.copyBand(srcBandName, sourceProduct, srcBandName, targetProduct, true);
            } else if (srcBandName.contains("azOffset") || srcBandName.contains("rgOffset")) {
                continue;
            } else {
                targetBand = new Band(srcBandName, band.getDataType(), band.getRasterWidth(),
                        band.getRasterHeight());

                targetBand.setUnit(band.getUnit());
                targetProduct.addBand(targetBand);
            }

            if (targetBand != null && srcBandName.startsWith("q_")) {
                final String suffix = srcBandName.substring(1);
                ReaderUtils.createVirtualIntensityBand(targetProduct, targetProduct.getBand("i" + suffix),
                        targetBand, suffix);
            }
        }

        targetProduct.setPreferredTileSize(512, subSwath[subSwathIndex - 1].linesPerBurst);
    }

    /**
     * Called by the framework in order to compute a tile for the given target band.
     * <p>The default implementation throws a runtime exception with the message "not implemented".</p>
     *
     * @param targetTileMap   The target tiles associated with all target bands to be computed.
     * @param targetRectangle The rectangle of target tile.
     * @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 raster.
     */
    @Override
    public void computeTileStack(Map<Band, Tile> targetTileMap, Rectangle targetRectangle, ProgressMonitor pm)
            throws OperatorException {

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

        try {
            if (!isOffsetAvailable) {
                estimateOffset();
            }

            Band slvBandI = null, slvBandQ = null;
            Band tgtBandI = null, tgtBandQ = null;
            Band derampDemodPhaseBand = null;
            final Band[] sourceBands = sourceProduct.getBands();
            for (Band band : sourceBands) {
                final String bandName = band.getName();
                if (bandName.contains("i_") && bandName.contains(StackUtils.SLV)) {
                    slvBandI = band;
                    tgtBandI = targetProduct.getBand(bandName);
                }

                if (bandName.contains("q_") && bandName.contains(StackUtils.SLV)) {
                    slvBandQ = band;
                    tgtBandQ = targetProduct.getBand(bandName);
                }

                if (bandName.contains(DerampDemodPhase)) {
                    derampDemodPhaseBand = band;
                }
            }

            // get deramp/demodulation phase
            final Tile derampDemodPhaseTile = getSourceTile(derampDemodPhaseBand, targetRectangle);
            final ProductData derampDemodPhaseData = derampDemodPhaseTile.getDataBuffer();
            final TileIndex index = new TileIndex(derampDemodPhaseTile);
            final double[][] derampDemodPhase = new double[h][w];
            for (int y = y0; y < yMax; y++) {
                index.calculateStride(y);
                final int yy = y - y0;
                for (int x = x0; x < xMax; x++) {
                    final int idx = index.getIndex(x);
                    derampDemodPhase[yy][x - x0] = derampDemodPhaseData.getElemDoubleAt(idx);
                }
            }

            // perform deramp and demodulation
            final Tile slvTileI = getSourceTile(slvBandI, targetRectangle);
            final Tile slvTileQ = getSourceTile(slvBandQ, targetRectangle);
            final double[][] derampDemodI = new double[h][w];
            final double[][] derampDemodQ = new double[h][w];
            BackGeocodingOp.performDerampDemod(slvTileI, slvTileQ, targetRectangle, derampDemodPhase, derampDemodI,
                    derampDemodQ);

            // compute shift phase
            final double[] phase = new double[2 * h];
            computeShiftPhaseArray(azOffset, h, phase);

            // perform azimuth shift using FFT, and perform reramp and remodulation
            final Tile tgtTileI = targetTileMap.get(tgtBandI);
            final Tile tgtTileQ = targetTileMap.get(tgtBandQ);
            final ProductData tgtDataI = tgtTileI.getDataBuffer();
            final ProductData tgtDataQ = tgtTileQ.getDataBuffer();

            final double[] col1 = new double[2 * h];
            final double[] col2 = new double[2 * h];
            final DoubleFFT_1D col_fft = new DoubleFFT_1D(h);
            for (int c = 0; c < w; c++) {
                final int x = x0 + c;
                for (int r = 0; r < h; r++) {
                    col1[2 * r] = derampDemodI[r][c];
                    col1[2 * r + 1] = derampDemodQ[r][c];

                    col2[2 * r] = derampDemodPhase[r][c];
                    col2[2 * r + 1] = 0.0;
                }

                col_fft.complexForward(col1);
                col_fft.complexForward(col2);

                multiplySpectrumByShiftFactor(col1, phase);
                multiplySpectrumByShiftFactor(col2, phase);

                col_fft.complexInverse(col1, true);
                col_fft.complexInverse(col2, true);

                for (int r = 0; r < h; r++) {
                    final int y = y0 + r;
                    final double cosPhase = FastMath.cos(col2[2 * r]);
                    final double sinPhase = FastMath.sin(col2[2 * r]);
                    final int idx = tgtTileI.getDataBufferIndex(x, y);
                    tgtDataI.setElemDoubleAt(idx, (float) (col1[2 * r] * cosPhase + col1[2 * r + 1] * sinPhase));
                    tgtDataQ.setElemDoubleAt(idx, (float) (-col1[2 * r] * sinPhase + col1[2 * r + 1] * cosPhase));
                }
            }

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

    /**
     * Estimate azimuth offset using ESD approach.
     * @throws Exception The exception.
     */
    private synchronized void estimateOffset() throws Exception {

        if (isOffsetAvailable) {
            return;
        }

        final int[] overlapSizeArray = computeBurstOverlapSize();
        final int numOverlaps = overlapSizeArray.length;
        final List<Double> azOffsetArray = new ArrayList<>(numOverlaps);
        final List<Integer> overlapIndexArray = new ArrayList<>(numOverlaps);

        final StatusProgressMonitor status = new StatusProgressMonitor(StatusProgressMonitor.TYPE.SUBTASK);
        status.beginTask("Estimating azimuth offset... ", numOverlaps);

        final ThreadManager threadManager = new ThreadManager();
        try {
            for (int i = 0; i < numOverlaps; i++) {
                checkForCancellation();
                final int x0 = 0;
                final int y0 = subSwath[subSwathIndex - 1].linesPerBurst * (i + 1);
                final int w = subSwath[subSwathIndex - 1].samplesPerBurst;
                final int h = overlapSizeArray[i];
                final int overlapIndex = i;

                final double tCycle = subSwath[subSwathIndex - 1].linesPerBurst
                        * subSwath[subSwathIndex - 1].azimuthTimeInterval;

                double sumSpectralSeparation = 0.0;
                for (int b = 0; b < subSwath[subSwathIndex - 1].numOfBursts; b++) {
                    for (int p = 0; p < subSwath[subSwathIndex - 1].samplesPerBurst; p++) {
                        sumSpectralSeparation += subSwath[subSwathIndex - 1].dopplerRate[b][p] * tCycle;
                    }
                }
                final double spectralSeparation = sumSpectralSeparation
                        / (subSwath[subSwathIndex - 1].numOfBursts * subSwath[subSwathIndex - 1].samplesPerBurst);

                final Thread worker = new Thread() {
                    @Override
                    public void run() {
                        try {
                            final Rectangle backwardRectangle = new Rectangle(x0, y0, w, h);
                            final Rectangle forwardRectangle = new Rectangle(x0, y0 - h, w, h);
                            final Band mBandI = getBand(StackUtils.MST, "i_", swathIndexStr, polarizations[0]);
                            final Band mBandQ = getBand(StackUtils.MST, "q_", swathIndexStr, polarizations[0]);
                            final Band sBandI = getBand(StackUtils.SLV, "i_", swathIndexStr, polarizations[0]);
                            final Band sBandQ = getBand(StackUtils.SLV, "q_", swathIndexStr, polarizations[0]);

                            final double azOffset = estimateAzOffsets(mBandI, mBandQ, sBandI, sBandQ,
                                    backwardRectangle, forwardRectangle, spectralSeparation);
                            //System.out.println("azOffset = " + azOffset);

                            synchronized (azOffsetArray) {
                                azOffsetArray.add(azOffset);
                                overlapIndexArray.add(overlapIndex);
                            }
                        } catch (Throwable e) {
                            OperatorUtils.catchOperatorException("estimateOffset", e);
                        }
                    }
                };
                threadManager.add(worker);
                status.worked(1);
            }
            status.done();
            threadManager.finish();

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

        // todo The following simple average should be replaced by weighted average using coherence as weight
        double sumAzOffset = 0.0;
        for (int i = 0; i < azOffsetArray.size(); i++) {
            final double anAzOffset = azOffsetArray.get(i);
            sumAzOffset += anAzOffset;
            SystemUtils.LOG.info("AzimuthShiftOp: overlap area = " + overlapIndexArray.get(i)
                    + ", azimuth offset = " + anAzOffset);
        }
        azOffset = sumAzOffset / numOverlaps;
        SystemUtils.LOG.info("AzimuthShiftOp: whole image azimuth offset = " + azOffset);

        saveAzimuthOffsetToMetadata(overlapIndexArray, azOffsetArray);

        isOffsetAvailable = true;
    }

    private void saveAzimuthOffsetToMetadata(final List<Integer> overlapIndexArray,
            final List<Double> azOffsetArray) {

        final MetadataElement absTgt = AbstractMetadata.getAbstractedMetadata(targetProduct);
        if (absTgt == null) {
            return;
        }

        final double azimuthPixelSpacing = absTgt.getAttributeDouble(AbstractMetadata.azimuth_spacing);
        final MetadataElement ESDMeasurement = new MetadataElement("ESD_Measurement");
        final MetadataElement swathElem = new MetadataElement(subSwathNames[0]);
        swathElem.addAttribute(new MetadataAttribute("count", ProductData.TYPE_INT16));
        swathElem.setAttributeInt("count", overlapIndexArray.size());

        for (int i = 0; i < azOffsetArray.size(); i++) {
            final MetadataElement overlapListElem = new MetadataElement("OverlapList." + i);
            overlapListElem.addAttribute(new MetadataAttribute("azimuthShift", ProductData.TYPE_FLOAT32));
            overlapListElem.setAttributeDouble("azimuthShift", azOffsetArray.get(i) * azimuthPixelSpacing * 100.0);
            overlapListElem.addAttribute(new MetadataAttribute("overlapIndex", ProductData.TYPE_INT16));
            overlapListElem.setAttributeInt("overlapIndex", overlapIndexArray.get(i));
            swathElem.addElement(overlapListElem);
        }

        ESDMeasurement.addElement(swathElem);
        absTgt.addElement(ESDMeasurement);
    }

    /**
     * Compute burst overlap size for all bursts in given sub-swath.
     * @return The burst overlap size array.
     */
    private int[] computeBurstOverlapSize() {

        final int numBursts = subSwath[subSwathIndex - 1].numOfBursts;
        int[] sizeArray = new int[numBursts - 1];

        for (int i = 0; i < numBursts - 1; i++) {
            final double endTime = subSwath[subSwathIndex - 1].burstLastLineTime[i];
            final double startTime = subSwath[subSwathIndex - 1].burstFirstLineTime[i + 1];
            sizeArray[i] = (int) ((endTime - startTime) / subSwath[subSwathIndex - 1].azimuthTimeInterval);
        }

        return sizeArray;
    }

    private double estimateAzOffsets(final Band mBandI, final Band mBandQ, final Band sBandI, final Band sBandQ,
            final Rectangle backwardRectangle, final Rectangle forwardRectangle, final double spectralSeparation) {

        final int mDataType = mBandI.getDataType();
        final int sDataType = sBandI.getDataType();

        final Tile mTileIBack = getSourceTile(mBandI, backwardRectangle);
        final Tile mTileQBack = getSourceTile(mBandQ, backwardRectangle);
        final Tile sTileIBack = getSourceTile(sBandI, backwardRectangle);
        final Tile sTileQBack = getSourceTile(sBandQ, backwardRectangle);

        double[] mIBackArray = null;
        double[] mQBackArray = null;
        if (mDataType == ProductData.TYPE_INT16) {
            final short[] mIBackArrayShort = (short[]) mTileIBack.getDataBuffer().getElems();
            final short[] mQBackArrayShort = (short[]) mTileQBack.getDataBuffer().getElems();
            mIBackArray = new double[mIBackArrayShort.length];
            mQBackArray = new double[mQBackArrayShort.length];
            for (int i = 0; i < mIBackArrayShort.length; i++) {
                mIBackArray[i] = (double) mIBackArrayShort[i];
                mQBackArray[i] = (double) mQBackArrayShort[i];
            }
        } else {
            mIBackArray = (double[]) mTileIBack.getDataBuffer().getElems();
            mQBackArray = (double[]) mTileQBack.getDataBuffer().getElems();
        }

        double[] sIBackArray, sQBackArray;
        if (sDataType == ProductData.TYPE_FLOAT32) {
            final float[] sIBackArrayFloat = (float[]) sTileIBack.getDataBuffer().getElems();
            final float[] sQBackArrayFloat = (float[]) sTileQBack.getDataBuffer().getElems();
            sIBackArray = new double[sIBackArrayFloat.length];
            sQBackArray = new double[sQBackArrayFloat.length];
            for (int i = 0; i < sIBackArrayFloat.length; i++) {
                sIBackArray[i] = (double) sIBackArrayFloat[i];
                sQBackArray[i] = (double) sQBackArrayFloat[i];
            }
        } else {
            sIBackArray = (double[]) sTileIBack.getDataBuffer().getElems();
            sQBackArray = (double[]) sTileQBack.getDataBuffer().getElems();
        }

        final Tile mTileIFor = getSourceTile(mBandI, forwardRectangle);
        final Tile mTileQFor = getSourceTile(mBandQ, forwardRectangle);
        final Tile sTileIFor = getSourceTile(sBandI, forwardRectangle);
        final Tile sTileQFor = getSourceTile(sBandQ, forwardRectangle);

        double[] mIForArray = null;
        double[] mQForArray = null;
        if (mDataType == ProductData.TYPE_INT16) {
            final short[] mIForArrayShort = (short[]) mTileIFor.getDataBuffer().getElems();
            final short[] mQForArrayShort = (short[]) mTileQFor.getDataBuffer().getElems();
            mIForArray = new double[mIForArrayShort.length];
            mQForArray = new double[mQForArrayShort.length];
            for (int i = 0; i < mIForArrayShort.length; i++) {
                mIForArray[i] = (double) mIForArrayShort[i];
                mQForArray[i] = (double) mQForArrayShort[i];
            }
        } else {
            mIForArray = (double[]) mTileIFor.getDataBuffer().getElems();
            mQForArray = (double[]) mTileQFor.getDataBuffer().getElems();
        }

        double[] sIForArray = null;
        double[] sQForArray = null;
        if (sDataType == ProductData.TYPE_FLOAT32) {
            final float[] sIForArrayFloat = (float[]) sTileIFor.getDataBuffer().getElems();
            final float[] sQForArrayFloat = (float[]) sTileQFor.getDataBuffer().getElems();
            sIForArray = new double[sIForArrayFloat.length];
            sQForArray = new double[sQForArrayFloat.length];
            for (int i = 0; i < sIForArrayFloat.length; i++) {
                sIForArray[i] = (double) sIForArrayFloat[i];
                sQForArray[i] = (double) sQForArrayFloat[i];
            }
        } else {
            sIForArray = (double[]) sTileIFor.getDataBuffer().getElems();
            sQForArray = (double[]) sTileQFor.getDataBuffer().getElems();
        }

        final int arrayLength = mIBackArray.length;
        final double[] backIntReal = new double[arrayLength];
        final double[] backIntImag = new double[arrayLength];
        complexArrayMultiplication(mIBackArray, mQBackArray, sIBackArray, sQBackArray, backIntReal, backIntImag);

        final double[] forIntReal = new double[arrayLength];
        final double[] forIntImag = new double[arrayLength];
        complexArrayMultiplication(mIForArray, mQForArray, sIForArray, sQForArray, forIntReal, forIntImag);

        final double[] diffIntReal = new double[arrayLength];
        final double[] diffIntImag = new double[arrayLength];
        complexArrayMultiplication(forIntReal, forIntImag, backIntReal, backIntImag, diffIntReal, diffIntImag);

        double sumReal = 0.0;
        double sumImag = 0.0;
        for (int i = 0; i < arrayLength; i++) {
            final double theta = Math.atan2(diffIntImag[i], diffIntReal[i]);
            sumReal += FastMath.cos(theta);
            sumImag += FastMath.sin(theta);
        }

        final double phase = Math.atan2(sumImag, sumReal);
        return phase / (2 * Math.PI * spectralSeparation * subSwath[subSwathIndex - 1].azimuthTimeInterval);
    }

    private void complexArrayMultiplication(final double[] realArray1, final double[] imagArray1,
            final double[] realArray2, final double[] imagArray2, final double[] realOutput,
            final double[] imagOutput) {

        final int arrayLength = realArray1.length;
        if (imagArray1.length != arrayLength || realArray2.length != arrayLength || imagArray2.length != arrayLength
                || realOutput.length != arrayLength || imagOutput.length != arrayLength) {
            throw new OperatorException("Arrays of the same length are expected.");
        }

        for (int i = 0; i < arrayLength; i++) {
            realOutput[i] = realArray1[i] * realArray2[i] + imagArray1[i] * imagArray2[i];
            imagOutput[i] = imagArray1[i] * realArray2[i] - realArray1[i] * imagArray2[i];
        }
    }

    private Band getBand(final String suffix, final String prefix, final String swathIndexStr,
            final String polarization) {

        final String[] bandNames = sourceProduct.getBandNames();
        for (String bandName : bandNames) {
            if (bandName.contains(suffix) && bandName.contains(prefix) && bandName.contains(swathIndexStr)
                    && bandName.contains(polarization)) {
                return sourceProduct.getBand(bandName);
            }
        }
        return null;
    }

    private static void computeShiftPhaseArray(final double shift, final int signalLength,
            final double[] phaseArray) {

        int k2;
        double phaseK;
        final double phase = -2.0 * Math.PI * shift / signalLength;
        final int halfSignalLength = (int) (signalLength * 0.5 + 0.5);

        for (int k = 0; k < signalLength; ++k) {
            if (k < halfSignalLength) {
                phaseK = phase * k;
            } else {
                phaseK = phase * (k - signalLength);
            }
            k2 = k * 2;
            phaseArray[k2] = FastMath.cos(phaseK);
            phaseArray[k2 + 1] = FastMath.sin(phaseK);
        }
    }

    private static void multiplySpectrumByShiftFactor(final double[] array, final double[] phaseArray) {

        int k2;
        double c, s;
        double real, imag;
        final int signalLength = array.length / 2;
        for (int k = 0; k < signalLength; ++k) {
            k2 = k * 2;
            c = phaseArray[k2];
            s = phaseArray[k2 + 1];
            real = array[k2];
            imag = array[k2 + 1];
            array[k2] = real * c - imag * s;
            array[k2 + 1] = real * s + imag * c;
        }
    }

    /**
     * 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(AzimuthShiftOp.class);
        }
    }

}