Search in sources :

Example 1 with ComplexDoubleMatrix

use of org.jblas.ComplexDoubleMatrix in project s1tbx by senbox-org.

the class CoherenceOp method computePartialTile.

private void computePartialTile(final int subSwathIndex, final int burstIndex, final int firstLineIdx, final Rectangle targetRectangle, final Map<Band, Tile> targetTileMap) throws OperatorException {
    try {
        final BorderExtender border = BorderExtender.createInstance(BorderExtender.BORDER_ZERO);
        final int y0 = targetRectangle.y;
        final int yN = y0 + targetRectangle.height - 1;
        final int x0 = targetRectangle.x;
        final int xN = x0 + targetRectangle.width - 1;
        final int cohx0 = targetRectangle.x - (cohWinRg - 1) / 2;
        final int cohy0 = targetRectangle.y - (cohWinAz - 1) / 2;
        final int cohw = targetRectangle.width + cohWinRg - 1;
        final int cohh = targetRectangle.height + cohWinAz - 1;
        final Rectangle extRect = new Rectangle(cohx0, cohy0, cohw, cohh);
        final org.jlinda.core.Window tileWindow = new org.jlinda.core.Window(cohy0 - firstLineIdx, cohy0 + cohh - 1 - firstLineIdx, cohx0, cohx0 + cohw - 1);
        final SLCImage mstMeta = targetMap.values().iterator().next().sourceMaster.metaData.clone();
        updateMstMetaData(burstIndex, mstMeta);
        final Orbit mstOrbit = targetMap.values().iterator().next().sourceMaster.orbit;
        DemTile demTile = null;
        if (subtractTopographicPhase) {
            demTile = TopoPhase.getDEMTile(tileWindow, mstMeta, mstOrbit, dem, demNoDataValue, demSamplingLat, demSamplingLon, tileExtensionPercent);
            if (demTile == null) {
                throw new OperatorException("The selected DEM has no overlap with the image or is invalid.");
            }
            if (demTile.getData().length < 3 || demTile.getData()[0].length < 3) {
                throw new OperatorException("The resolution of the selected DEM is too low, " + "please select DEM with higher resolution.");
            }
        }
        final int minLine = 0;
        final int maxLine = subSwath[subSwathIndex - 1].linesPerBurst - 1;
        final int minPixel = 0;
        final int maxPixel = subSwath[subSwathIndex - 1].samplesPerBurst - 1;
        for (String cohKey : targetMap.keySet()) {
            final ProductContainer product = targetMap.get(cohKey);
            final SLCImage slvMeta = product.sourceSlave.metaData.clone();
            updateSlvMetaData(product, burstIndex, slvMeta);
            final Orbit slvOrbit = product.sourceSlave.orbit;
            final Tile mstTileReal = getSourceTile(product.sourceMaster.realBand, extRect, border);
            final Tile mstTileImag = getSourceTile(product.sourceMaster.imagBand, extRect, border);
            final ComplexDoubleMatrix dataMaster = TileUtilsDoris.pullComplexDoubleMatrix(mstTileReal, mstTileImag);
            final Tile slvTileReal = getSourceTile(product.sourceSlave.realBand, extRect, border);
            final Tile slvTileImag = getSourceTile(product.sourceSlave.imagBand, extRect, border);
            final ComplexDoubleMatrix dataSlave = TileUtilsDoris.pullComplexDoubleMatrix(slvTileReal, slvTileImag);
            final String polynomialName = product.sourceSlave.name + '_' + (subSwathIndex - 1) + '_' + burstIndex;
            if (subtractFlatEarthPhase) {
                final DoubleMatrix flatEarthPhase = computeFlatEarthPhase(cohx0, cohx0 + cohw - 1, cohw, cohy0 - firstLineIdx, cohy0 + cohh - 1 - firstLineIdx, cohh, minPixel, maxPixel, minLine, maxLine, polynomialName);
                final ComplexDoubleMatrix complexReferencePhase = new ComplexDoubleMatrix(MatrixFunctions.cos(flatEarthPhase), MatrixFunctions.sin(flatEarthPhase));
                dataSlave.muli(complexReferencePhase);
                if (OUTPUT_PHASE) {
                    saveFlatEarthPhase(x0, xN, y0, yN, flatEarthPhase, product, targetTileMap);
                }
            }
            if (subtractTopographicPhase) {
                TopoPhase topoPhase = TopoPhase.computeTopoPhase(mstMeta, mstOrbit, slvMeta, slvOrbit, tileWindow, demTile, false);
                final ComplexDoubleMatrix ComplexTopoPhase = new ComplexDoubleMatrix(MatrixFunctions.cos(new DoubleMatrix(topoPhase.demPhase)), MatrixFunctions.sin(new DoubleMatrix(topoPhase.demPhase)));
                dataSlave.muli(ComplexTopoPhase);
                if (OUTPUT_PHASE) {
                    saveTopoPhase(x0, xN, y0, yN, topoPhase.demPhase, product, targetTileMap);
                }
            }
            for (int i = 0; i < dataMaster.length; i++) {
                double tmp = norm(dataMaster.get(i));
                dataMaster.put(i, dataMaster.get(i).mul(dataSlave.get(i).conj()));
                dataSlave.put(i, new ComplexDouble(norm(dataSlave.get(i)), tmp));
            }
            DoubleMatrix cohMatrix = SarUtils.coherence2(dataMaster, dataSlave, cohWinAz, cohWinRg);
            saveCoherence(cohMatrix, product, targetTileMap, targetRectangle);
        }
    } catch (Throwable e) {
        OperatorUtils.catchOperatorException(getId(), e);
    }
}
Also used : Orbit(org.jlinda.core.Orbit) BorderExtender(javax.media.jai.BorderExtender) DemTile(org.jlinda.core.geom.DemTile) Tile(org.esa.snap.core.gpf.Tile) Point(org.jlinda.core.Point) GeoPoint(org.jlinda.core.GeoPoint) DoubleMatrix(org.jblas.DoubleMatrix) ComplexDoubleMatrix(org.jblas.ComplexDoubleMatrix) ComplexDouble(org.jblas.ComplexDouble) DemTile(org.jlinda.core.geom.DemTile) TopoPhase(org.jlinda.core.geom.TopoPhase) SLCImage(org.jlinda.core.SLCImage) OperatorException(org.esa.snap.core.gpf.OperatorException) ComplexDoubleMatrix(org.jblas.ComplexDoubleMatrix)

Example 2 with ComplexDoubleMatrix

use of org.jblas.ComplexDoubleMatrix in project s1tbx by senbox-org.

the class CoherenceOp method computeTileForNormalProduct.

private void computeTileForNormalProduct(final Map<Band, Tile> targetTileMap, Rectangle targetRectangle, ProgressMonitor pm) throws OperatorException {
    try {
        final BorderExtender border = BorderExtender.createInstance(BorderExtender.BORDER_ZERO);
        final int y0 = targetRectangle.y;
        final int yN = y0 + targetRectangle.height - 1;
        final int x0 = targetRectangle.x;
        final int xN = targetRectangle.x + targetRectangle.width - 1;
        // System.out.println("x0 = " + x0 +", y0 = " + y0 + ", w = " + targetRectangle.width + ", h = " + targetRectangle.height);
        final int cohx0 = targetRectangle.x - (cohWinRg - 1) / 2;
        final int cohy0 = targetRectangle.y - (cohWinAz - 1) / 2;
        final int cohw = targetRectangle.width + cohWinRg - 1;
        final int cohh = targetRectangle.height + cohWinAz - 1;
        final Rectangle extRect = new Rectangle(cohx0, cohy0, cohw, cohh);
        final org.jlinda.core.Window tileWindow = new org.jlinda.core.Window(cohy0, cohy0 + cohh - 1, cohx0, cohx0 + cohw - 1);
        DemTile demTile = null;
        if (subtractTopographicPhase) {
            demTile = TopoPhase.getDEMTile(tileWindow, targetMap, dem, demNoDataValue, demSamplingLat, demSamplingLon, tileExtensionPercent);
            if (demTile.getData().length < 3 || demTile.getData()[0].length < 3) {
                throw new OperatorException("The resolution of the selected DEM is too low, " + "please select DEM with higher resolution.");
            }
        }
        for (String cohKey : targetMap.keySet()) {
            final ProductContainer product = targetMap.get(cohKey);
            final Tile mstTileReal = getSourceTile(product.sourceMaster.realBand, extRect, border);
            final Tile mstTileImag = getSourceTile(product.sourceMaster.imagBand, extRect, border);
            final ComplexDoubleMatrix dataMaster = TileUtilsDoris.pullComplexDoubleMatrix(mstTileReal, mstTileImag);
            final Tile slvTileReal = getSourceTile(product.sourceSlave.realBand, extRect, border);
            final Tile slvTileImag = getSourceTile(product.sourceSlave.imagBand, extRect, border);
            final ComplexDoubleMatrix dataSlave = TileUtilsDoris.pullComplexDoubleMatrix(slvTileReal, slvTileImag);
            if (subtractFlatEarthPhase) {
                final DoubleMatrix flatEarthPhase = computeFlatEarthPhase(cohx0, cohx0 + cohw - 1, cohw, cohy0, cohy0 + cohh - 1, cohh, 0, sourceImageWidth - 1, 0, sourceImageHeight - 1, product.sourceSlave.name);
                final ComplexDoubleMatrix complexReferencePhase = new ComplexDoubleMatrix(MatrixFunctions.cos(flatEarthPhase), MatrixFunctions.sin(flatEarthPhase));
                dataSlave.muli(complexReferencePhase);
                if (OUTPUT_PHASE) {
                    saveFlatEarthPhase(x0, xN, y0, yN, flatEarthPhase, product, targetTileMap);
                }
            }
            if (subtractTopographicPhase) {
                final TopoPhase topoPhase = TopoPhase.computeTopoPhase(product, tileWindow, demTile, false);
                final ComplexDoubleMatrix ComplexTopoPhase = new ComplexDoubleMatrix(MatrixFunctions.cos(new DoubleMatrix(topoPhase.demPhase)), MatrixFunctions.sin(new DoubleMatrix(topoPhase.demPhase)));
                dataSlave.muli(ComplexTopoPhase);
                if (OUTPUT_PHASE) {
                    saveTopoPhase(x0, xN, y0, yN, topoPhase.demPhase, product, targetTileMap);
                }
            }
            for (int i = 0; i < dataMaster.length; i++) {
                double tmp = norm(dataMaster.get(i));
                dataMaster.put(i, dataMaster.get(i).mul(dataSlave.get(i).conj()));
                dataSlave.put(i, new ComplexDouble(norm(dataSlave.get(i)), tmp));
            }
            DoubleMatrix cohMatrix = SarUtils.coherence2(dataMaster, dataSlave, cohWinAz, cohWinRg);
            saveCoherence(cohMatrix, product, targetTileMap, targetRectangle);
        }
    } catch (Throwable e) {
        OperatorUtils.catchOperatorException(getId(), e);
    } finally {
        pm.done();
    }
}
Also used : BorderExtender(javax.media.jai.BorderExtender) DemTile(org.jlinda.core.geom.DemTile) Tile(org.esa.snap.core.gpf.Tile) Point(org.jlinda.core.Point) GeoPoint(org.jlinda.core.GeoPoint) DoubleMatrix(org.jblas.DoubleMatrix) ComplexDoubleMatrix(org.jblas.ComplexDoubleMatrix) ComplexDouble(org.jblas.ComplexDouble) DemTile(org.jlinda.core.geom.DemTile) TopoPhase(org.jlinda.core.geom.TopoPhase) OperatorException(org.esa.snap.core.gpf.OperatorException) ComplexDoubleMatrix(org.jblas.ComplexDoubleMatrix)

Example 3 with ComplexDoubleMatrix

use of org.jblas.ComplexDoubleMatrix in project s1tbx by senbox-org.

the class RangeShiftOp method getFineOffsets.

private void getFineOffsets(final PixelPos mGCPPixelPos, final PixelPos sGCPPixelPos, final double[] offset) {
    try {
        ComplexDoubleMatrix mI = getComplexDoubleMatrix(mstBandI, mstBandQ, mGCPPixelPos, fineWinWidth, fineWinHeight);
        ComplexDoubleMatrix sI = getComplexDoubleMatrix(slvBandI, slvBandQ, sGCPPixelPos, fineWinWidth, fineWinHeight);
        final double[] fineOffset = { 0, 0 };
        final double coherence = CoregistrationUtils.crossCorrelateFFT(fineOffset, mI, sI, fineWinOvsFactor, fineWinAccY, fineWinAccX);
        if (coherence < xCorrThreshold) {
            offset[0] = noDataValue;
            offset[1] = noDataValue;
        } else {
            offset[0] = -fineOffset[0];
            offset[1] = -fineOffset[1];
        }
    } catch (Throwable e) {
        OperatorUtils.catchOperatorException(getId() + " getFineOffsets ", e);
    }
}
Also used : ComplexDoubleMatrix(org.jblas.ComplexDoubleMatrix)

Example 4 with ComplexDoubleMatrix

use of org.jblas.ComplexDoubleMatrix in project s1tbx by senbox-org.

the class DInSAROp method computeTileStack.

@Override
public void computeTileStack(Map<Band, Tile> targetTileMap, Rectangle targetRectangle, ProgressMonitor pm) throws OperatorException {
    try {
        int y0 = targetRectangle.y;
        int yN = y0 + targetRectangle.height - 1;
        int x0 = targetRectangle.x;
        int xN = targetRectangle.x + targetRectangle.width - 1;
        final Window tileWindow = new Window(y0, yN, x0, xN);
        Band targetBand_I;
        Band targetBand_Q;
        ComplexDoubleMatrix complexDefoPair = null;
        DoubleMatrix doubleTopoPair = null;
        ProductContainer product = null;
        for (String ifgKey : targetMap.keySet()) {
            product = targetMap.get(ifgKey);
            // / check out results from source ///
            Tile tileReal = getSourceTile(product.sourceMaster.realBand, targetRectangle);
            Tile tileImag = getSourceTile(product.sourceMaster.imagBand, targetRectangle);
            complexDefoPair = TileUtilsDoris.pullComplexDoubleMatrix(tileReal, tileImag);
        }
        // always pull from topoProduct
        for (Band band : topoProduct.getBands()) {
            if (band.getName().contains("unw") || band.getUnit().contains(Unit.ABS_PHASE)) {
                // / check out results from source ///
                Tile tileReal = getSourceTile(band, targetRectangle);
                doubleTopoPair = TileUtilsDoris.pullDoubleMatrix(tileReal);
            }
        }
        dinsar.applyDInSAR(tileWindow, complexDefoPair, doubleTopoPair);
        // / commit complexDefoPair back to target ///
        targetBand_I = targetProduct.getBand(product.targetBandName_I);
        Tile tileOutReal = targetTileMap.get(targetBand_I);
        TileUtilsDoris.pushDoubleMatrix(complexDefoPair.real(), tileOutReal, targetRectangle);
        targetBand_Q = targetProduct.getBand(product.targetBandName_Q);
        Tile tileOutImag = targetTileMap.get(targetBand_Q);
        TileUtilsDoris.pushDoubleMatrix(complexDefoPair.imag(), tileOutImag, targetRectangle);
    } catch (Throwable e) {
        OperatorUtils.catchOperatorException(getId(), e);
    }
}
Also used : Window(org.jlinda.core.Window) DoubleMatrix(org.jblas.DoubleMatrix) ComplexDoubleMatrix(org.jblas.ComplexDoubleMatrix) ProductContainer(org.jlinda.core.utils.ProductContainer) Tile(org.esa.snap.core.gpf.Tile) Band(org.esa.snap.core.datamodel.Band) ComplexDoubleMatrix(org.jblas.ComplexDoubleMatrix)

Example 5 with ComplexDoubleMatrix

use of org.jblas.ComplexDoubleMatrix in project s1tbx by senbox-org.

the class PhaseFilterOp method computeTileStack.

@Override
public void computeTileStack(Map<Band, Tile> targetTileMap, Rectangle targetRectangle, ProgressMonitor pm) throws OperatorException {
    try {
        // String filterType = "goldstein";
        // final float alpha = (float) 0.75;
        // final int overlap = 8;
        // final double[] kernel = new double[]{1. / 5, 1. / 5, 1. / 5, 1. / 5, 1. / 5};
        // final int blockSize = 64;
        Band targetBand_I;
        Band targetBand_Q;
        final Rectangle rectIn = new Rectangle(targetRectangle);
        final Rectangle rectOut = new Rectangle(targetRectangle);
        rectIn.y -= (overlap);
        rectIn.height += (2 * overlap);
        rectIn.x -= (overlap);
        rectIn.width += (2 * overlap);
        for (String ifgKey : targetMap.keySet()) {
            ProductContainer product = targetMap.get(ifgKey);
            Tile tileReal = getSourceTile(product.sourceMaster.realBand, rectIn);
            Tile tileImag = getSourceTile(product.sourceMaster.imagBand, rectIn);
            // put interferogram together
            ComplexDoubleMatrix complexIfg = TileUtilsDoris.pullComplexDoubleMatrix(tileReal, tileImag);
            PhaseFilter phaseFilter = new PhaseFilter(method, complexIfg, blockSize, overlap, kernelArray, alpha);
            phaseFilter.filter();
            complexIfg = phaseFilter.getData();
            // commit to target [note the tile overlap and boundary because of filter]
            // output [x0,y0] is shifted because of the tile overlap
            targetBand_I = targetProduct.getBand(product.targetBandName_I);
            Tile tileOutReal = targetTileMap.get(targetBand_I);
            TileUtilsDoris.pushDoubleMatrix(complexIfg.real(), tileOutReal, rectOut, overlap, overlap);
            targetBand_Q = targetProduct.getBand(product.targetBandName_Q);
            Tile tileOutImag = targetTileMap.get(targetBand_Q);
            TileUtilsDoris.pushDoubleMatrix(complexIfg.imag(), tileOutImag, rectOut, overlap, overlap);
        }
    } catch (Exception e) {
        throw new OperatorException(e);
    }
}
Also used : ProductContainer(org.jlinda.core.utils.ProductContainer) PhaseFilter(org.jlinda.core.filtering.PhaseFilter) Tile(org.esa.snap.core.gpf.Tile) Band(org.esa.snap.core.datamodel.Band) OperatorException(org.esa.snap.core.gpf.OperatorException) OperatorException(org.esa.snap.core.gpf.OperatorException) ComplexDoubleMatrix(org.jblas.ComplexDoubleMatrix)

Aggregations

ComplexDoubleMatrix (org.jblas.ComplexDoubleMatrix)112 DoubleMatrix (org.jblas.DoubleMatrix)40 Test (org.junit.Test)40 Window (org.jlinda.core.Window)24 ComplexDouble (org.jblas.ComplexDouble)20 Tile (org.esa.snap.core.gpf.Tile)8 File (java.io.File)6 OperatorException (org.esa.snap.core.gpf.OperatorException)5 SLCImage (org.jlinda.core.SLCImage)5 BorderExtender (javax.media.jai.BorderExtender)4 Band (org.esa.snap.core.datamodel.Band)4 DemTile (org.jlinda.core.geom.DemTile)4 GeoPoint (org.jlinda.core.GeoPoint)3 Point (org.jlinda.core.Point)3 TopoPhase (org.jlinda.core.geom.TopoPhase)3 ProductContainer (org.jlinda.core.utils.ProductContainer)3 DoubleFFT_2D (edu.emory.mathcs.jtransforms.fft.DoubleFFT_2D)2 IOException (java.io.IOException)2 ProductData (org.esa.snap.core.datamodel.ProductData)2 FloatMatrix (org.jblas.FloatMatrix)2