Examples with ComplexDouble org.jblas.ComplexDouble used on opensource projects

Example 1 with ComplexDouble

use of org.jblas.ComplexDouble 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);
    }
}

Example 2 with ComplexDouble

use of org.jblas.ComplexDouble 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();
    }
}

Example 3 with ComplexDouble

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

the class DInSARTest method toComplexDoubleMatrix.

private static ComplexDoubleMatrix toComplexDoubleMatrix(final int rows, final int cols, final float[] in) {
    ComplexDoubleMatrix out = new ComplexDoubleMatrix(rows, cols);
    for (int i = 0; i < rows; i++) {
        int cnt = 0;
        for (int j = 0; j < 2 * cols; j = j + 2) {
            double re = in[i * (2 * cols) + j];
            double im = in[i * (2 * cols) + (j + 1)];
            out.put(i, cnt, new ComplexDouble(re, im));
            cnt++;
        }
    }
    return out;
}

Example 4 with ComplexDouble

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

the class LinearAlgebraUtilsTest method setUpTestMatrices.

@BeforeClass
public static void setUpTestMatrices() throws Exception {
    A_PASCAL_22 = DoubleMatrix.ones(2, 2);
    A_PASCAL_22.put(1, 1, 2);
    A_PASCAL_33 = DoubleMatrix.ones(3, 3);
    A_PASCAL_33.put(1, 1, 2);
    A_PASCAL_33.put(1, 2, 3);
    A_PASCAL_33.put(2, 1, 3);
    A_PASCAL_33.put(2, 2, 6);
    A_PASCAL_33_SQUARED = new double[][] { { 1, 1, 1 }, { 1, 8, 27 }, { 1, 27, 216 } };
    A_PASCAL_33_CHOL_EXPECTED = A_PASCAL_33.dup();
    // define lower triangular block
    A_PASCAL_33_CHOL_EXPECTED.put(1, 0, 1);
    A_PASCAL_33_CHOL_EXPECTED.put(1, 1, 1);
    A_PASCAL_33_CHOL_EXPECTED.put(2, 0, 1);
    A_PASCAL_33_CHOL_EXPECTED.put(2, 1, 2);
    A_PASCAL_33_CHOL_EXPECTED.put(2, 2, 1);
    // define inverted matrix
    A_PASCAL_33_CHOL_INV_EXPECTED.put(0, 0, 3);
    A_PASCAL_33_CHOL_INV_EXPECTED.put(0, 1, -3);
    A_PASCAL_33_CHOL_INV_EXPECTED.put(0, 2, 1);
    A_PASCAL_33_CHOL_INV_EXPECTED.put(1, 0, -3);
    A_PASCAL_33_CHOL_INV_EXPECTED.put(1, 1, 5);
    A_PASCAL_33_CHOL_INV_EXPECTED.put(1, 2, -2);
    A_PASCAL_33_CHOL_INV_EXPECTED.put(2, 0, 1);
    A_PASCAL_33_CHOL_INV_EXPECTED.put(2, 1, -2);
    A_PASCAL_33_CHOL_INV_EXPECTED.put(2, 2, 1);
    // define complex PASCAL_22 matrix
    A_PASCAL_22_CPLX = new ComplexDoubleMatrix(A_PASCAL_22, A_PASCAL_22);
    // A_PASCAL_22_CPLX_times2_EXPECTED.put(0, 0, new ComplexDouble(0, 4));
    // A_PASCAL_22_CPLX_times2_EXPECTED.put(0, 1, new ComplexDouble(0, 6));
    // A_PASCAL_22_CPLX_times2_EXPECTED.put(1, 0, new ComplexDouble(0, 6));
    // A_PASCAL_22_CPLX_times2_EXPECTED.put(1, 1, new ComplexDouble(0, 10));
    A_PASCAL_22_CPLX_times2_EXPECTED.put(0, 0, new ComplexDouble(0, 2));
    A_PASCAL_22_CPLX_times2_EXPECTED.put(0, 1, new ComplexDouble(0, 2));
    A_PASCAL_22_CPLX_times2_EXPECTED.put(1, 0, new ComplexDouble(0, 2));
    A_PASCAL_22_CPLX_times2_EXPECTED.put(1, 1, new ComplexDouble(0, 8));
    // A_33
    A_33.put(0, 0, 1);
    A_33.put(0, 1, 2);
    A_33.put(0, 2, 3);
    A_33.put(1, 0, 4);
    A_33.put(1, 1, 5);
    A_33.put(1, 2, 6);
    A_33.put(2, 0, 7);
    A_33.put(2, 1, 8);
    A_33.put(2, 2, 9);
    // A_33
    A_33.put(0, 0, 1);
    A_33.put(0, 1, 2);
    A_33.put(0, 2, 3);
    A_33.put(1, 0, 4);
    A_33.put(1, 1, 5);
    A_33.put(1, 2, 6);
    A_33.put(2, 0, 7);
    A_33.put(2, 1, 8);
    A_33.put(2, 2, 9);
    // ATA_33
    ATA_33_EXPECTED.put(0, 0, 66);
    ATA_33_EXPECTED.put(0, 1, 78);
    ATA_33_EXPECTED.put(0, 2, 90);
    ATA_33_EXPECTED.put(1, 0, 78);
    ATA_33_EXPECTED.put(1, 1, 93);
    ATA_33_EXPECTED.put(1, 2, 108);
    ATA_33_EXPECTED.put(2, 0, 90);
    ATA_33_EXPECTED.put(2, 1, 108);
    ATA_33_EXPECTED.put(2, 2, 126);
    // A_33_FLIPED_EXPECTED
    A_33_FLIP_EXPECTED.putColumn(0, A_33.getColumn(2));
    A_33_FLIP_EXPECTED.putColumn(1, A_33.getColumn(1));
    A_33_FLIP_EXPECTED.putColumn(2, A_33.getColumn(0));
    // A_15_EXPECTED
    A_15 = new DoubleMatrix(MathUtils.increment(5, 0, 1));
    // A_15_SHIFTED_EXPECTED
    A_15_SHIFT_EXPECTED = new DoubleMatrix(new double[] { 3, 4, 0, 1, 2 });
}

Example 5 with ComplexDouble

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

the class PhaseFilter method smoothSpectral.

// Do the same as smoothSpace but faster   -----> for Goldstein filter ??????
// some overhead due to conversion r4<->cr4
private DoubleMatrix smoothSpectral(final DoubleMatrix data, final int blockSize) {
    final int nRows = data.rows;
    final int nCols = data.columns;
    // init to zero...
    final ComplexDoubleMatrix smoothData = new ComplexDoubleMatrix(nRows, nCols);
    // or define fft(R4)
    SpectralUtils.fft2D_inplace(smoothData);
    // init to zeros
    ComplexDoubleMatrix kernel = new ComplexDoubleMatrix(1, nRows);
    // design 1d kernel function of block
    for (int ii = -blockSize; ii <= blockSize; ++ii) {
        kernel.put(0, (ii + nRows) % nRows, new ComplexDouble(1.0 / (2 * blockSize + 1), 0.0));
    }
    ComplexDoubleMatrix kernel2d = LinearAlgebraUtils.matTxmat(kernel, kernel);
    // should be real sinc
    SpectralUtils.fft2D_inplace(kernel2d);
    LinearAlgebraUtils.dotmult(smoothData, kernel2d);
    // convolution, but still complex...
    SpectralUtils.invfft2D_inplace(smoothData);
    return smoothData.real();
}