Java tutorial
/* * 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); } } }