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.sar.gpf.geometric; import com.bc.ceres.core.ProgressMonitor; import org.apache.commons.math3.util.FastMath; import org.esa.snap.core.datamodel.Band; import org.esa.snap.core.datamodel.CrsGeoCoding; import org.esa.snap.core.datamodel.GeoCoding; import org.esa.snap.core.datamodel.GeoPos; import org.esa.snap.core.datamodel.MetadataElement; import org.esa.snap.core.datamodel.PixelPos; import org.esa.snap.core.datamodel.Product; import org.esa.snap.core.datamodel.ProductData; import org.esa.snap.core.datamodel.RasterDataNode; import org.esa.snap.core.datamodel.Stx; import org.esa.snap.core.datamodel.VirtualBand; import org.esa.snap.core.dataop.resamp.Resampling; import org.esa.snap.core.dataop.resamp.ResamplingFactory; 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.Parameter; import org.esa.snap.core.gpf.annotations.SourceProducts; 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.core.util.math.MathUtils; import org.esa.snap.engine_utilities.datamodel.AbstractMetadata; import org.esa.snap.engine_utilities.datamodel.Unit; import org.esa.snap.engine_utilities.eo.Constants; import org.esa.snap.engine_utilities.gpf.OperatorUtils; import org.esa.snap.engine_utilities.gpf.TileGeoreferencing; import org.esa.snap.engine_utilities.gpf.TileIndex; import org.geotools.geometry.jts.ReferencedEnvelope; import org.geotools.referencing.crs.DefaultGeographicCRS; import org.geotools.referencing.wkt.UnformattableObjectException; import org.opengis.referencing.crs.CoordinateReferenceSystem; import java.awt.*; import java.awt.geom.Rectangle2D; import java.text.MessageFormat; import java.util.ArrayList; import java.util.HashMap; import java.util.HashSet; import java.util.List; import java.util.Map; import java.util.Set; /** * The Mosaic operator. */ @OperatorMetadata(alias = "SAR-Mosaic", category = "Radar/Geometric", authors = "Jun Lu, Luis Veci", version = "1.0", copyright = "Copyright (C) 2014 by Array Systems Computing Inc.", description = "Mosaics two or more products based on their geo-codings.") public class MosaicOp extends Operator { @SourceProducts private Product[] sourceProduct; @TargetProduct private Product targetProduct = null; @Parameter(description = "The list of source bands.", alias = "sourceBands", label = "Source Bands") private String[] sourceBandNames = null; @Parameter(valueSet = { ResamplingFactory.NEAREST_NEIGHBOUR_NAME, ResamplingFactory.BILINEAR_INTERPOLATION_NAME, ResamplingFactory.CUBIC_CONVOLUTION_NAME }, defaultValue = ResamplingFactory.NEAREST_NEIGHBOUR_NAME, description = "The method to be used when resampling the slave grid onto the master grid.", label = "Resampling Type") private String resamplingMethod = ResamplingFactory.NEAREST_NEIGHBOUR_NAME; @Parameter(defaultValue = "true", description = "Average the overlapping areas", label = "Average Overlap") private Boolean average = true; @Parameter(defaultValue = "true", description = "Normalize by Mean", label = "Normalize by Mean") private Boolean normalizeByMean = true; @Parameter(defaultValue = "false", description = "Gradient Domain Mosaic", label = "Gradient Domain Mosaic") private Boolean gradientDomainMosaic = false; @Parameter(defaultValue = "0", description = "Pixel Size (m)", label = "Pixel Size (m)") private double pixelSize = 0; @Parameter(defaultValue = "0", description = "Target width", label = "Scene Width (pixels)") private int sceneWidth = 0; @Parameter(defaultValue = "0", description = "Target height", label = "Scene Height (pixels)") private int sceneHeight = 0; @Parameter(defaultValue = "0", description = "Feather amount around source image", label = "Feature (pixels)") private int feather = 0; @Parameter(defaultValue = "5000", description = "Maximum number of iterations", label = "Maximum Iterations") private int maxIterations = 5000; @Parameter(defaultValue = "1e-4", description = "Convergence threshold for Relaxed Gauss-Seidel method", label = "Convergence Threshold") private double convergenceThreshold = 1e-4; private final SceneProperties scnProp = new SceneProperties(); private final Map<Integer, Band> bandIndexSet = new HashMap<>(20); private final Map<Product, Rectangle> srcRectMap = new HashMap<>(10); private Product[] selectedProducts = null; private boolean outputGradientBand = false; @Override public void initialize() throws OperatorException { try { if (gradientDomainMosaic && resamplingMethod.contains(ResamplingFactory.NEAREST_NEIGHBOUR_NAME)) { throw new OperatorException("Nearest neighbour resampling method produces poor result with gradient" + " domain mosaic, please select other method"); } GeoCoding srcGeocoding = null; for (final Product prod : sourceProduct) { if (prod.getSceneGeoCoding() == null) { throw new OperatorException( MessageFormat.format("Product ''{0}'' has no geo-coding.", prod.getName())); } if (srcGeocoding == null) { srcGeocoding = prod.getSceneGeoCoding(); } // todo check that all source products have same geocoding } getSourceBands(); computeImageGeoBoundary(selectedProducts, scnProp); if (sceneWidth == 0 || sceneHeight == 0) { final MetadataElement absRoot = AbstractMetadata.getAbstractedMetadata(sourceProduct[0]); if (pixelSize == 0 && absRoot != null) { final double rangeSpacing = AbstractMetadata.getAttributeDouble(absRoot, AbstractMetadata.range_spacing); final double azimuthSpacing = AbstractMetadata.getAttributeDouble(absRoot, AbstractMetadata.azimuth_spacing); pixelSize = Math.min(rangeSpacing, azimuthSpacing); } getSceneDimensions(pixelSize, scnProp); sceneWidth = scnProp.sceneWidth; sceneHeight = scnProp.sceneHeight; final double ratio = sceneWidth / (double) sceneHeight; long dim = (long) sceneWidth * (long) sceneHeight; while (sceneWidth > 0 && sceneHeight > 0 && dim > Integer.MAX_VALUE) { sceneWidth -= 1000; sceneHeight = (int) (sceneWidth / ratio); dim = (long) sceneWidth * (long) sceneHeight; } } targetProduct = new Product("mosaic", "mosaic", sceneWidth, sceneHeight); targetProduct.setSceneGeoCoding(createCRSGeoCoding(srcGeocoding)); // add all bands for (Map.Entry<Integer, Band> integerBandEntry : bandIndexSet.entrySet()) { final Band srcBand = integerBandEntry.getValue(); int targetBandDataType; if (gradientDomainMosaic) { targetBandDataType = ProductData.TYPE_FLOAT32; } else { targetBandDataType = srcBand.getDataType(); } final Band targetBand = new Band(srcBand.getName(), targetBandDataType, sceneWidth, sceneHeight); targetBand.setUnit(srcBand.getUnit()); targetBand.setDescription(srcBand.getDescription()); targetBand.setNoDataValue(srcBand.getNoDataValue()); targetBand.setNoDataValueUsed(true); targetProduct.addBand(targetBand); if (gradientDomainMosaic && outputGradientBand) { // add gradient band final String targetBandName = srcBand.getName() + "_gradient"; final Band gradientBand = new Band(targetBandName, ProductData.TYPE_FLOAT32, sceneWidth, sceneHeight); targetProduct.addBand(gradientBand); } } // transfer index coding if (sourceProduct[0].getIndexCodingGroup().getNodeCount() > 0 && sourceProduct[0].getIndexCodingGroup().get(0) != null) { try { ProductUtils.copyIndexCodings(sourceProduct[0], targetProduct); } catch (Exception e) { if (!resamplingMethod.equals(Resampling.NEAREST_NEIGHBOUR)) { throw new OperatorException( "Use Nearest Neighbour with Classificaitons: " + e.getMessage()); } } } for (Product srcProduct : selectedProducts) { final Rectangle srcRect = getSrcRect(targetProduct.getSceneGeoCoding(), scnProp.srcCornerLatitudeMap.get(srcProduct), scnProp.srcCornerLongitudeMap.get(srcProduct)); srcRectMap.put(srcProduct, srcRect); } updateTargetProductMetadata(); } catch (Throwable e) { OperatorUtils.catchOperatorException(getId(), e); } } private CrsGeoCoding createCRSGeoCoding(GeoCoding srcGeocoding) throws Exception { final CoordinateReferenceSystem srcCRS = srcGeocoding.getMapCRS(); String wkt; try { wkt = srcCRS.toWKT(); } catch (UnformattableObjectException e) { // if too complex to convert using strict wkt = srcCRS.toString(); } final CoordinateReferenceSystem targetCRS = CRSGeoCodingHandler.getCRS(wkt); final double pixelSpacingInDegree = pixelSize / Constants.semiMajorAxis * Constants.RTOD; double pixelSizeX = pixelSize; double pixelSizeY = pixelSize; if (targetCRS.getName().getCode().equals("WGS84(DD)")) { pixelSizeX = pixelSpacingInDegree; pixelSizeY = pixelSpacingInDegree; } final Rectangle2D bounds = new Rectangle2D.Double(); bounds.setFrameFromDiagonal(scnProp.lonMin, scnProp.latMin, scnProp.lonMax, scnProp.latMax); final ReferencedEnvelope boundsEnvelope = new ReferencedEnvelope(bounds, DefaultGeographicCRS.WGS84); final ReferencedEnvelope targetEnvelope = boundsEnvelope.transform(targetCRS, true); return new CrsGeoCoding(targetCRS, sceneWidth, sceneHeight, targetEnvelope.getMinimum(0), targetEnvelope.getMaximum(1), pixelSizeX, pixelSizeY); } /** * Update metadata in the target product. */ private void updateTargetProductMetadata() { MetadataElement absRoot = AbstractMetadata.getAbstractedMetadata(targetProduct); if (absRoot == null) { absRoot = AbstractMetadata.addAbstractedMetadataHeader(targetProduct.getMetadataRoot()); } AbstractMetadata.setAttribute(absRoot, AbstractMetadata.range_spacing, pixelSize); AbstractMetadata.setAttribute(absRoot, AbstractMetadata.azimuth_spacing, pixelSize); } private Band[] getSourceBands() throws OperatorException { final List<Band> bandList = new ArrayList<>(20); final Set<Product> selectedProductSet = new HashSet<>(sourceProduct.length); if (sourceBandNames == null || sourceBandNames.length == 0) { for (final Product srcProduct : sourceProduct) { for (final Band band : srcProduct.getBands()) { if (band instanceof VirtualBand) continue; bandList.add(band); bandIndexSet.put(srcProduct.getBandIndex(band.getName()), band); } selectedProductSet.add(srcProduct); } } else { for (final String name : sourceBandNames) { final String bandName = getBandName(name); final String productName = getProductName(name, sourceProduct[0].getName()); final Product srcProduct = getProduct(productName); final Band band = srcProduct.getBand(bandName); final String bandUnit = band.getUnit(); if (bandUnit != null) { if (bandUnit.contains(Unit.IMAGINARY) || bandUnit.contains(Unit.REAL)) { throw new OperatorException("Real and imaginary bands not handled"); } } bandList.add(band); bandIndexSet.put(srcProduct.getBandIndex(band.getName()), band); selectedProductSet.add(srcProduct); } } selectedProducts = selectedProductSet.toArray(new Product[selectedProductSet.size()]); return bandList.toArray(new Band[bandList.size()]); } private Product getProduct(final String productName) { for (Product prod : sourceProduct) { if (prod.getName().equals(productName)) { return prod; } } return null; } private static String getBandName(final String name) { if (name.contains("::")) return name.substring(0, name.indexOf("::")); return name; } private static String getProductName(final String name, final String defaultName) { if (name.contains("::")) return name.substring(name.indexOf("::") + 2, name.length()); return defaultName; } private static Rectangle getSrcRect(final GeoCoding destGeoCoding, final double[] lats, final double[] lons) { double srcLatMin = 90.0; double srcLatMax = -90.0; double srcLonMin = 180.0; double srcLonMax = -180.0; for (double lat : lats) { if (lat < srcLatMin) { srcLatMin = lat; } if (lat > srcLatMax) { srcLatMax = lat; } } for (double lon : lons) { if (lon < srcLonMin) { srcLonMin = lon; } if (lon > srcLonMax) { srcLonMax = lon; } } final GeoPos geoPos = new GeoPos(); final PixelPos[] pixelPos = new PixelPos[4]; geoPos.setLocation(srcLatMin, srcLonMin); pixelPos[0] = destGeoCoding.getPixelPos(geoPos, null); geoPos.setLocation(srcLatMin, srcLonMax); pixelPos[1] = destGeoCoding.getPixelPos(geoPos, null); geoPos.setLocation(srcLatMax, srcLonMax); pixelPos[2] = destGeoCoding.getPixelPos(geoPos, null); geoPos.setLocation(srcLatMax, srcLonMin); pixelPos[3] = destGeoCoding.getPixelPos(geoPos, null); return getBoundingBox(pixelPos, 0, 0, Integer.MAX_VALUE, Integer.MAX_VALUE, 4); } private static Rectangle getBoundingBox(final PixelPos[] pixelPositions, final int minOffsetX, final int minOffsetY, final int maxWidth, final int maxHeight, final int margin) { int minX = Integer.MAX_VALUE; int maxX = -Integer.MAX_VALUE; int minY = Integer.MAX_VALUE; int maxY = -Integer.MAX_VALUE; for (final PixelPos pixelsPos : pixelPositions) { if (pixelsPos != null) { final int x = (int) Math.floor(pixelsPos.getX()); final int y = (int) Math.floor(pixelsPos.getY()); if (x < minX) { minX = x; } if (x > maxX) { maxX = x; } if (y < minY) { minY = y; } if (y > maxY) { maxY = y; } } } if (minX > maxX || minY > maxY) { return null; } minX = Math.max(minX - margin, minOffsetX); maxX = Math.min(maxX + margin, maxWidth - 1); minY = Math.max(minY - margin, minOffsetY); maxY = Math.min(maxY + margin, maxHeight - 1); if (minX > maxX || minY > maxY) { return null; } int w = maxX - minX + 1; int h = maxY - minY + 1; return new Rectangle(minX, minY, w, h); } /** * 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 List<Product> validProducts = new ArrayList<>(sourceProduct.length); for (final Product srcProduct : selectedProducts) { final Rectangle srcRect = srcRectMap.get(srcProduct); if (srcRect == null || !srcRect.intersects(targetRectangle)) { continue; } validProducts.add(srcProduct); } if (validProducts.isEmpty()) { return; } final GeoPos geoPos = new GeoPos(); final PixelPos pixelPos = new PixelPos(); final int minX = targetRectangle.x; final int minY = targetRectangle.y; final int maxX = targetRectangle.x + targetRectangle.width - 1; final int maxY = targetRectangle.y + targetRectangle.height - 1; final TileGeoreferencing tileGeoRef = new TileGeoreferencing(targetProduct, minX, minY, maxX - minX + 1, maxY - minY + 1); final List<PixelPos[]> srcPixelCoords = new ArrayList<>(validProducts.size()); final int numPixelPos = targetRectangle.width * targetRectangle.height; for (Product validProduct : validProducts) { srcPixelCoords.add(new PixelPos[numPixelPos]); } int coordIndex = 0; int prodIndex; for (int y = minY; y <= maxY; ++y) { for (int x = minX; x <= maxX; ++x) { tileGeoRef.getGeoPos(x, y, geoPos); prodIndex = 0; for (final Product srcProduct : validProducts) { srcProduct.getSceneGeoCoding().getPixelPos(geoPos, pixelPos); if (pixelPos.x >= feather && pixelPos.y >= feather && pixelPos.x < srcProduct.getSceneRasterWidth() - feather && pixelPos.y < srcProduct.getSceneRasterHeight() - feather) { srcPixelCoords.get(prodIndex)[coordIndex] = new PixelPos(pixelPos.x, pixelPos.y); } else { srcPixelCoords.get(prodIndex)[coordIndex] = null; } ++prodIndex; } ++coordIndex; } } final Resampling resampling = ResamplingFactory.createResampling(resamplingMethod); if (gradientDomainMosaic) { performGradientDomainMosaic(targetTiles, targetRectangle, srcPixelCoords, validProducts, resampling, pm); return; } final List<SourceData> validSourceData = new ArrayList<>(validProducts.size()); for (final Map.Entry<Band, Tile> bandTileEntry : targetTiles.entrySet()) { final String trgBandName = bandTileEntry.getKey().getName(); validSourceData.clear(); prodIndex = 0; for (final Product srcProduct : validProducts) { final Band srcBand = srcProduct.getBand(trgBandName); if (srcBand == null) { continue; } final PixelPos[] pixPos = srcPixelCoords.get(prodIndex); final Rectangle sourceRectangle = getBoundingBox(pixPos, feather, feather, srcProduct.getSceneRasterWidth() - feather, srcProduct.getSceneRasterHeight() - feather, 4); if (sourceRectangle != null) { double min = 0, max = 0, mean = 0, std = 0; if (normalizeByMean) { // get stat values try { final Stx stats = srcBand.getStx(true, ProgressMonitor.NULL); mean = stats.getMean(); min = stats.getMinimum(); max = stats.getMaximum(); std = stats.getStandardDeviation(); } catch (Throwable e) { //OperatorUtils.catchOperatorException(getId(), e); normalizeByMean = false; // temporary disable } } try { final Tile srcTile = getSourceTile(srcBand, sourceRectangle); if (srcTile != null) { validSourceData .add(new SourceData(srcTile, pixPos, resampling, min, max, mean, std)); } } catch (Exception e) { SystemUtils.LOG.severe("Mosaic getSourceTile failed " + e.getMessage()); //continue } } ++prodIndex; } if (!validSourceData.isEmpty()) { collocateSourceBand(validSourceData, resampling, bandTileEntry.getValue()); } } } catch (Throwable e) { OperatorUtils.catchOperatorException(getId(), e); } finally { pm.done(); } } private void collocateSourceBand(final List<SourceData> validSourceData, final Resampling resampling, final Tile targetTile) throws OperatorException { try { final Rectangle targetRectangle = targetTile.getRectangle(); final ProductData trgBuffer = targetTile.getDataBuffer(); double sample; final int maxY = targetRectangle.y + targetRectangle.height; final int maxX = targetRectangle.x + targetRectangle.width; final TileIndex trgIndex = new TileIndex(targetTile); final double[] sampleList = new double[validSourceData.size()]; final int[] sampleDistanceList = new int[validSourceData.size()]; for (int y = targetRectangle.y, index = 0; y < maxY; ++y) { trgIndex.calculateStride(y); for (int x = targetRectangle.x; x < maxX; ++x, ++index) { double targetVal = 0; int numSamples = 0; for (final SourceData srcDat : validSourceData) { final PixelPos sourcePixelPos = srcDat.srcPixPos[index]; if (sourcePixelPos == null) { continue; } resampling.computeIndex(sourcePixelPos.x, sourcePixelPos.y, srcDat.srcRasterWidth - feather, srcDat.srcRasterHeight - feather, srcDat.resamplingIndex); sample = resampling.resample(srcDat.resamplingRaster, srcDat.resamplingIndex); if (!Double.isNaN(sample) && sample != srcDat.nodataValue && !MathUtils.equalValues(sample, 0.0F, 1e-4F)) { if (normalizeByMean) { sample -= srcDat.srcMean; sample /= srcDat.srcStd; } targetVal = sample; if (average) { sampleList[numSamples] = sample; sampleDistanceList[numSamples] = (int) (Math.min(sourcePixelPos.x + 1, srcDat.srcRasterWidth - sourcePixelPos.x) * Math.min(sourcePixelPos.y + 1, srcDat.srcRasterHeight - sourcePixelPos.y)); numSamples++; } } } if (targetVal != 0) { if (average && numSamples > 1) { double sum = 0; int totalWeight = 0; for (int i = 0; i < numSamples; i++) { sum += sampleList[i] * sampleDistanceList[i]; totalWeight += sampleDistanceList[i]; } targetVal = sum / totalWeight; } trgBuffer.setElemDoubleAt(trgIndex.getIndex(x), targetVal); } } } } catch (Throwable e) { OperatorUtils.catchOperatorException(getId(), e); } } private void performGradientDomainMosaic(final Map<Band, Tile> targetTiles, final Rectangle targetRectangle, final List<PixelPos[]> srcPixelCoords, final List<Product> validProducts, final Resampling resampling, ProgressMonitor pm) throws OperatorException { try { final int minX = targetRectangle.x; final int minY = targetRectangle.y; final int maxX = targetRectangle.x + targetRectangle.width - 1; final int maxY = targetRectangle.y + targetRectangle.height - 1; double[][] mosaicedTile = new double[targetRectangle.height][targetRectangle.width]; double[][] gradientTile = new double[targetRectangle.height][targetRectangle.width]; byte[][] mask = new byte[targetRectangle.height][targetRectangle.width]; // -1: no data, 0: used by existing product, 1: used by new product, 2: need mosaic final List<SourceData> validSourceData = new ArrayList<>(validProducts.size()); // loop through all target bands for (final Map.Entry<Band, Tile> bandTileEntry : targetTiles.entrySet()) { final String trgBandName = bandTileEntry.getKey().getName(); if (trgBandName.contains("_gradient")) { continue; } final Tile trgTile = bandTileEntry.getValue(); final ProductData trgBuffer = trgTile.getDataBuffer(); // for each target band, get source data for all related source bands getValidSourceData(validProducts, trgBandName, srcPixelCoords, resampling, validSourceData, pm); // for each target band, find all related source bands and use them in mosaic // for now we assume that source products have been sorted according to time with the oldest first for (int i = 0; i < validSourceData.size(); i++) { if (i == 0) { readFirstProduct(minX, maxX, minY, maxY, validSourceData.get(i), resampling, mosaicedTile, mask); } else { readNextProduct(minX, maxX, minY, maxY, validSourceData.get(i), resampling, mosaicedTile, mask, gradientTile); performMosaic(mask, gradientTile, mosaicedTile); cleanUpMask(mask); } } // save mosaiced image final TileIndex trgIndex = new TileIndex(trgTile); for (int y = minY; y <= maxY; y++) { trgIndex.calculateStride(y); for (int x = minX; x <= maxX; x++) { trgBuffer.setElemDoubleAt(trgIndex.getIndex(x), mosaicedTile[y - minY][x - minX]); } } // save gradient if (outputGradientBand) { final Band gradientBand = targetProduct.getBand(trgBandName + "_gradient"); final ProductData gradientBuffer = targetTiles.get(gradientBand).getDataBuffer(); for (int y = minY; y <= maxY; y++) { trgIndex.calculateStride(y); for (int x = minX; x <= maxX; x++) { gradientBuffer.setElemDoubleAt(trgIndex.getIndex(x), gradientTile[y - minY][x - minX]); } } } } } catch (Throwable e) { OperatorUtils.catchOperatorException(getId(), e); } } private void getValidSourceData(final List<Product> validProducts, final String trgBandName, final List<PixelPos[]> srcPixelCoords, final Resampling resampling, List<SourceData> validSourceData, ProgressMonitor pm) { try { validSourceData.clear(); int prodIndex = 0; for (final Product srcProduct : validProducts) { final Band srcBand = srcProduct.getBand(trgBandName); if (srcBand == null) { continue; } final PixelPos[] pixPos = srcPixelCoords.get(prodIndex); final Rectangle sourceRectangle = getBoundingBox(pixPos, 0, 0, srcProduct.getSceneRasterWidth(), srcProduct.getSceneRasterHeight(), feather); if (sourceRectangle != null) { double mean = 0, min = 0, max = 0, std = 0; if (normalizeByMean) { try { final Stx stats = srcBand.getStx(true, pm); mean = stats.getMean(); min = stats.getMinimum(); max = stats.getMaximum(); std = stats.getStandardDeviation(); } catch (Throwable e) { //OperatorUtils.catchOperatorException(getId(), e); normalizeByMean = false; } } try { final Tile srcTile = getSourceTile(srcBand, sourceRectangle); if (srcTile != null) { validSourceData.add(new SourceData(srcTile, pixPos, resampling, min, max, mean, std)); } } catch (Exception e) { SystemUtils.LOG.severe("Mosaic getSourceTile failed " + e.getMessage()); //continue } } ++prodIndex; } } catch (Throwable e) { OperatorUtils.catchOperatorException(getId(), e); } } private void readFirstProduct(final int minX, final int maxX, final int minY, final int maxY, final SourceData srcDat, final Resampling resampling, double[][] mosaicedTile, byte[][] mask) throws OperatorException { try { double sample; int yy, xx; for (int y = minY, index = 0; y <= maxY; ++y) { yy = y - minY; for (int x = minX; x <= maxX; ++x, ++index) { xx = x - minX; final PixelPos sourcePixelPos = srcDat.srcPixPos[index]; if (sourcePixelPos == null) { mosaicedTile[yy][xx] = srcDat.nodataValue; mask[yy][xx] = -1; continue; } resampling.computeIndex(sourcePixelPos.x, sourcePixelPos.y, srcDat.srcRasterWidth, srcDat.srcRasterHeight, srcDat.resamplingIndex); sample = resampling.resample(srcDat.resamplingRaster, srcDat.resamplingIndex); if (isValidSample(sample, srcDat.nodataValue)) { if (normalizeByMean) { sample -= srcDat.srcMean; sample /= srcDat.srcStd; } mosaicedTile[yy][xx] = sample; mask[yy][xx] = 0; } else { mosaicedTile[yy][xx] = srcDat.nodataValue; mask[yy][xx] = -1; } } } } catch (Throwable e) { OperatorUtils.catchOperatorException(getId(), e); } } private void readNextProduct(final int minX, final int maxX, final int minY, final int maxY, final SourceData srcDat, final Resampling resampling, double[][] mosaicedTile, byte[][] mask, double[][] gradientTile) throws OperatorException { try { final int targetTileWidth = mosaicedTile[0].length; final int targetTileHeight = mosaicedTile.length; double[] adjacentPixels = new double[4]; double sample; int yy, xx; for (int y = minY, index = 0; y <= maxY; ++y) { yy = y - minY; for (int x = minX; x <= maxX; ++x, ++index) { xx = x - minX; final PixelPos sourcePixelPos = srcDat.srcPixPos[index]; if (sourcePixelPos == null) { continue; } resampling.computeIndex(sourcePixelPos.x, sourcePixelPos.y, srcDat.srcRasterWidth, srcDat.srcRasterHeight, srcDat.resamplingIndex); sample = resampling.resample(srcDat.resamplingRaster, srcDat.resamplingIndex); if (isValidSample(sample, srcDat.nodataValue)) { if (normalizeByMean) { sample -= srcDat.srcMean; sample /= srcDat.srcStd; } if (mask[yy][xx] == -1) { mosaicedTile[yy][xx] = sample; mask[yy][xx] = 1; } else if (mask[yy][xx] == 0 && isInnerPoint(index, targetTileWidth, targetTileHeight, srcDat, resampling, adjacentPixels)) { if (isInnerPoint(xx, yy, mask)) { mask[yy][xx] = 2; mosaicedTile[yy][xx] = sample; //gradientTile[yy][xx] = computeGradient(xx, yy, mosaicedTile, sample, adjacentPixels); gradientTile[yy][xx] = adjacentPixels[0] + adjacentPixels[1] + adjacentPixels[2] + adjacentPixels[3] - 4 * sample; } else { mosaicedTile[yy][xx] = sample; } } } } } } catch (Throwable e) { OperatorUtils.catchOperatorException(getId(), e); } } private boolean isInnerPoint(final int index, final int targetTileWidth, final int targetTileHeight, final SourceData srcDat, final Resampling resampling, double[] adjacentPixels) { try { final int indexUp = index - targetTileWidth; final int indexDown = index + targetTileWidth; final int indexLeft = index - 1; final int indexRight = index + 1; if (indexUp >= 0 && indexDown < targetTileWidth * targetTileHeight && index % targetTileWidth != 0 && (index + 1) % targetTileWidth != 0 && srcDat.srcPixPos[indexUp] != null && srcDat.srcPixPos[indexDown] != null && srcDat.srcPixPos[indexLeft] != null && srcDat.srcPixPos[indexRight] != null) { resampling.computeIndex(srcDat.srcPixPos[indexUp].x, srcDat.srcPixPos[indexUp].y, srcDat.srcRasterWidth, srcDat.srcRasterHeight, srcDat.resamplingIndex); final double s1 = resampling.resample(srcDat.resamplingRaster, srcDat.resamplingIndex); resampling.computeIndex(srcDat.srcPixPos[indexDown].x, srcDat.srcPixPos[indexDown].y, srcDat.srcRasterWidth, srcDat.srcRasterHeight, srcDat.resamplingIndex); final double s2 = resampling.resample(srcDat.resamplingRaster, srcDat.resamplingIndex); resampling.computeIndex(srcDat.srcPixPos[indexLeft].x, srcDat.srcPixPos[indexLeft].y, srcDat.srcRasterWidth, srcDat.srcRasterHeight, srcDat.resamplingIndex); final double s3 = resampling.resample(srcDat.resamplingRaster, srcDat.resamplingIndex); resampling.computeIndex(srcDat.srcPixPos[indexRight].x, srcDat.srcPixPos[indexRight].y, srcDat.srcRasterWidth, srcDat.srcRasterHeight, srcDat.resamplingIndex); final double s4 = resampling.resample(srcDat.resamplingRaster, srcDat.resamplingIndex); if (isValidSample((float) s1, srcDat.nodataValue) && isValidSample((float) s2, srcDat.nodataValue) && isValidSample((float) s3, srcDat.nodataValue) && isValidSample((float) s4, srcDat.nodataValue)) { if (normalizeByMean) { adjacentPixels[0] = (s1 - srcDat.srcMean) / srcDat.srcStd; adjacentPixels[1] = (s2 - srcDat.srcMean) / srcDat.srcStd; adjacentPixels[2] = (s3 - srcDat.srcMean) / srcDat.srcStd; adjacentPixels[3] = (s4 - srcDat.srcMean) / srcDat.srcStd; } else { adjacentPixels[0] = s1; adjacentPixels[1] = s2; adjacentPixels[2] = s3; adjacentPixels[3] = s4; } return true; } else { return false; } } return false; } catch (Throwable e) { OperatorUtils.catchOperatorException(getId(), e); } return false; } private static boolean isInnerPoint(final int xx, final int yy, final byte[][] mask) { if (xx == 0 || yy == 0 || xx == mask[0].length - 1 || yy == mask.length - 1) { return false; } else { return (mask[yy - 1][xx] == 0 || mask[yy - 1][xx] == 2) && (mask[yy + 1][xx] == 0 || mask[yy + 1][xx] == 2) && (mask[yy][xx - 1] == 0 || mask[yy][xx - 1] == 2) && (mask[yy][xx + 1] == 0 || mask[yy][xx + 1] == 2); } } private static boolean isValidSample(final double sample, final double noDataValue) { return (!Double.isNaN(sample) && sample != noDataValue && !MathUtils.equalValues(sample, 0.0F, 1e-4F)); } private double computeGradient(final int xx, final int yy, final double[][] mosaicedTile, final double s0, final double[] adjacentPixels) { double g2 = adjacentPixels[0] + adjacentPixels[1] + adjacentPixels[2] + adjacentPixels[3] - 4 * s0; /* double g1 = mosaicedTile[yy-1][xx] + mosaicedTile[yy+1][xx] + mosaicedTile[yy][xx-1] + mosaicedTile[yy][xx+1] - 4*mosaicedTile[yy][xx]; if (Math.abs(g1) > Math.abs(g2)) { return g1; } else { return g2; } */ return g2; } private void performMosaic(final byte[][] mask, final double[][] gradientTile, double[][] mosaicedTile) { final double w = 1.5; final int rows = mask.length; final int cols = mask[0].length; double sigma, update, error = 0.0; int it; for (it = 0; it < maxIterations; it++) { error = 0.0; for (int r = 0; r < rows; r++) { for (int c = 0; c < cols; c++) { if (mask[r][c] == 2) { sigma = gradientTile[r][c] - mosaicedTile[r - 1][c] - mosaicedTile[r + 1][c] - mosaicedTile[r][c - 1] - mosaicedTile[r][c + 1]; update = (1 - w) * mosaicedTile[r][c] - w * sigma / 4.0; error = Math.max(error, Math.abs(mosaicedTile[r][c] - update)); mosaicedTile[r][c] = update; } } } if (error < convergenceThreshold) { break; } } //if(error != 0) // System.out.println("it = " + it + ", error = " + error); } private static void cleanUpMask(byte[][] mask) { final int rows = mask.length; final int cols = mask[0].length; for (int r = 0; r < rows; r++) { for (int c = 0; c < cols; c++) { if (mask[r][c] > 0) { mask[r][c] = 0; } } } } private static class ResamplingRaster implements Resampling.Raster { private final Tile tile; private final boolean usesNoData; private final boolean scalingApplied; private final double noDataValue; private final double geophysicalNoDataValue; private final ProductData dataBuffer; private final int minX, minY, maxX, maxY; public ResamplingRaster(final Tile tile) { this.tile = tile; this.minX = tile.getMinX(); this.minY = tile.getMinY(); this.maxX = tile.getMaxX(); this.maxY = tile.getMaxY(); this.dataBuffer = tile.getDataBuffer(); final RasterDataNode rasterDataNode = tile.getRasterDataNode(); this.usesNoData = rasterDataNode.isNoDataValueUsed(); this.noDataValue = rasterDataNode.getNoDataValue(); this.geophysicalNoDataValue = rasterDataNode.getGeophysicalNoDataValue(); this.scalingApplied = rasterDataNode.isScalingApplied(); } public final int getWidth() { return tile.getWidth(); } public final int getHeight() { return tile.getHeight(); } public boolean getSamples(final int[] x, final int[] y, final double[][] samples) throws Exception { boolean allValid = true; for (int i = 0; i < y.length; i++) { for (int j = 0; j < x.length; j++) { if (x[j] < minX || y[i] < minY || x[j] > maxX || y[i] > maxY) { samples[i][j] = Float.NaN; allValid = false; } try { samples[i][j] = dataBuffer.getElemDoubleAt(tile.getDataBufferIndex(x[j], y[i])); } catch (Exception e) { samples[i][j] = Float.NaN; allValid = false; } if (usesNoData) { if (scalingApplied && geophysicalNoDataValue == samples[i][j] || noDataValue == samples[i][j]) { samples[i][j] = Float.NaN; allValid = false; } } } } return allValid; } } /** * Compute source image geodetic boundary (minimum/maximum latitude/longitude) from the its corner * latitude/longitude. * * @param sourceProducts the list of input products * @param scnProp the output scene properties */ public static void computeImageGeoBoundary(final Product[] sourceProducts, final SceneProperties scnProp) { scnProp.latMin = 90.0f; scnProp.latMax = -90.0f; scnProp.lonMin = 180.0f; scnProp.lonMax = -180.0f; for (final Product srcProd : sourceProducts) { final GeoCoding geoCoding = srcProd.getSceneGeoCoding(); final GeoPos geoPosFirstNear = geoCoding.getGeoPos(new PixelPos(0, 0), null); final GeoPos geoPosFirstFar = geoCoding.getGeoPos(new PixelPos(srcProd.getSceneRasterWidth() - 1, 0), null); final GeoPos geoPosLastNear = geoCoding.getGeoPos(new PixelPos(0, srcProd.getSceneRasterHeight() - 1), null); final GeoPos geoPosLastFar = geoCoding.getGeoPos( new PixelPos(srcProd.getSceneRasterWidth() - 1, srcProd.getSceneRasterHeight() - 1), null); final double[] lats = { geoPosFirstNear.getLat(), geoPosFirstFar.getLat(), geoPosLastNear.getLat(), geoPosLastFar.getLat() }; final double[] lons = { geoPosFirstNear.getLon(), geoPosFirstFar.getLon(), geoPosLastNear.getLon(), geoPosLastFar.getLon() }; scnProp.srcCornerLatitudeMap.put(srcProd, lats); scnProp.srcCornerLongitudeMap.put(srcProd, lons); for (double lat : lats) { if (lat < scnProp.latMin) { scnProp.latMin = (float) lat; } if (lat > scnProp.latMax) { scnProp.latMax = (float) lat; } } for (double lon : lons) { if (lon < scnProp.lonMin) { scnProp.lonMin = (float) lon; } if (lon > scnProp.lonMax) { scnProp.lonMax = (float) lon; } } } } public static void getSceneDimensions(final double minSpacing, final SceneProperties scnProp) { double minAbsLat; if (scnProp.latMin * scnProp.latMax > 0) { minAbsLat = Math.min(Math.abs(scnProp.latMin), Math.abs(scnProp.latMax)) * Constants.DTOR; } else { minAbsLat = 0.0; } double delLat = minSpacing / Constants.MeanEarthRadius * Constants.RTOD; double delLon = minSpacing / (Constants.MeanEarthRadius * FastMath.cos(minAbsLat)) * Constants.RTOD; delLat = Math.min(delLat, delLon); delLon = delLat; scnProp.sceneWidth = (int) ((scnProp.lonMax - scnProp.lonMin) / delLon) + 1; scnProp.sceneHeight = (int) ((scnProp.latMax - scnProp.latMin) / delLat) + 1; } public static class SceneProperties { public int sceneWidth, sceneHeight; public float latMin, lonMin, latMax, lonMax; public final Map<Product, double[]> srcCornerLatitudeMap = new HashMap<>(10); public final Map<Product, double[]> srcCornerLongitudeMap = new HashMap<>(10); } private static class SourceData { final Tile srcTile; final ResamplingRaster resamplingRaster; final Resampling.Index resamplingIndex; final double nodataValue; final PixelPos[] srcPixPos; final int srcRasterHeight; final int srcRasterWidth; final double srcMean; final double srcMax; final double srcMin; final double srcStd; public SourceData(final Tile tile, final PixelPos[] pixPos, final Resampling resampling, final double min, final double max, final double mean, final double std) { srcTile = tile; resamplingRaster = new ResamplingRaster(srcTile); resamplingIndex = resampling.createIndex(); nodataValue = tile.getRasterDataNode().getNoDataValue(); srcPixPos = pixPos; final Product srcProduct = tile.getRasterDataNode().getProduct(); srcRasterHeight = srcProduct.getSceneRasterHeight(); srcRasterWidth = srcProduct.getSceneRasterWidth(); srcMin = min; srcMax = max; srcMean = mean; srcStd = std; } } /** * Operator SPI. */ public static class Spi extends OperatorSpi { public Spi() { super(MosaicOp.class); } } }