/*
* 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.calibration.gpf;
import com.bc.ceres.core.ProgressMonitor;
import org.esa.s1tbx.calibration.gpf.calibrators.Sentinel1Calibrator;
import org.esa.s1tbx.insar.gpf.support.Sentinel1Utils;
import org.esa.snap.core.datamodel.Band;
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.Parameter;
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.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.gpf.InputProductValidator;
import org.esa.snap.engine_utilities.gpf.OperatorUtils;
import org.esa.snap.engine_utilities.gpf.ThreadManager;
import org.esa.snap.engine_utilities.gpf.TileIndex;
import org.esa.snap.engine_utilities.util.Maths;
import java.awt.*;
import java.io.IOException;
import java.util.ArrayList;
import java.util.Arrays;
import java.util.List;
import java.util.Map;
import java.util.Set;
/**
* Mask "No-value" pixels for near-range region in a Sentinel-1 Level-1 GRD product.
* <p>
* [1] "Masking No-value Pixels on GRD Products generated by the Sentinel-1 ESA IPF",
* OI-MPC-OTH-0243-S1-Border-Masking-MPCS-916.
*/
@OperatorMetadata(alias = "Remove-GRD-Border-Noise",
category = "Radar/Sentinel-1 TOPS",
authors = "Jun Lu, Luis Veci",
copyright = "Copyright (C) 2014 by Array Systems Computing Inc.",
version = "1.0",
description = "Mask no-value pixels for GRD product")
public final class RemoveGRDBorderNoiseOp extends Operator {
@SourceProduct
private Product sourceProduct;
@TargetProduct
private Product targetProduct;
@Parameter(description = "The list of polarisations", label = "Polarisations")
private String[] selectedPolarisations;
@Parameter(description = "The border margin limit", defaultValue = "500", label = "Border margin limit[pixels]")
private int borderLimit = 500;
@Parameter(description = "The trim threshold", defaultValue = "0.5", label = "Threshold")
private double trimThreshold = 0.5;
private MetadataElement absRoot = null;
private MetadataElement origMetadataRoot = null;
private int sourceImageWidth = 0;
private int sourceImageHeight = 0;
private double version = 0.0f;
private double scalingFactor = 0.0;
private Double noDataValue = 0.0;
private String coPolarization = null;
private Sentinel1Utils.NoiseVector noiseVector = null;
private double[] noiseLUT = null;
private Band coPolBand = null;
private boolean thermalNoiseCorrectionPerformed = false;
private boolean useBorderDetection = true;
private boolean borderDetected = false;
private int topBorder = 0;
private int bottomBorder = 0;
private int leftBorder = 0;
private int rightBorder = 0;
/**
* Default constructor. The graph processing framework
* requires that an operator has a default constructor.
*/
public RemoveGRDBorderNoiseOp() {
}
/**
* 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.checkIfSentinel1Product();
validator.checkIfGRD();
validator.checkIfCalibrated(false);
absRoot = AbstractMetadata.getAbstractedMetadata(sourceProduct);
origMetadataRoot = AbstractMetadata.getOriginalProductMetadata(sourceProduct);
sourceImageWidth = sourceProduct.getSceneRasterWidth();
sourceImageHeight = sourceProduct.getSceneRasterHeight();
topBorder = borderLimit;
bottomBorder = sourceImageHeight - borderLimit;
leftBorder = borderLimit;
rightBorder = sourceImageWidth - borderLimit;
getIPFVersion();
getProductCoPolarization();
getThermalNoiseCorrectionFlag();
if (!thermalNoiseCorrectionPerformed) {
getThermalNoiseVector();
}
computeNoiseScalingFactor();
computeNoiseLUT();
createTargetProduct();
} catch (Throwable e) {
OperatorUtils.catchOperatorException(getId(), e);
}
}
/**
* Get IPF version.
*/
private void getIPFVersion() {
final String procSysId = absRoot.getAttributeString(AbstractMetadata.ProcessingSystemIdentifier);
version = Double.valueOf(procSysId.substring(procSysId.lastIndexOf(" ")));
}
/**
* Get product co-polarization.
*/
private void getProductCoPolarization() {
final String[] sourceBandNames = sourceProduct.getBandNames();
for (String bandName : sourceBandNames) {
if (bandName.contains("HH")) {
coPolarization = "HH";
coPolBand = sourceProduct.getBand(bandName);
break;
} else if (bandName.contains("VV")) {
coPolarization = "VV";
coPolBand = sourceProduct.getBand(bandName);
break;
}
}
if (coPolarization == null) {
throw new OperatorException("Input product does not contain band with HH or VV polarization");
}
noDataValue = coPolBand.getNoDataValue();
}
/**
* Get thermal noise correction flag from the original product metadata.
*/
private void getThermalNoiseCorrectionFlag() {
final MetadataElement annotationElem = origMetadataRoot.getElement("annotation");
final MetadataElement[] annotationDataSetListElem = annotationElem.getElements();
final MetadataElement productElem = annotationDataSetListElem[0].getElement("product");
final MetadataElement imageAnnotationElem = productElem.getElement("imageAnnotation");
final MetadataElement processingInformationElem = imageAnnotationElem.getElement("processingInformation");
thermalNoiseCorrectionPerformed = Boolean.parseBoolean(
processingInformationElem.getAttribute("thermalNoiseCorrectionPerformed").getData().getElemString());
}
/**
* Get the middle thermal noise vector from the original product metadata.
*/
private void getThermalNoiseVector() {
final MetadataElement noiseElem = origMetadataRoot.getElement("noise");
final MetadataElement[] noiseDataSetListElem = noiseElem.getElements();
for (MetadataElement dataSetListElem : noiseDataSetListElem) {
final MetadataElement noiElem = dataSetListElem.getElement("noise");
final MetadataElement adsHeaderElem = noiElem.getElement("adsHeader");
final String pol = adsHeaderElem.getAttributeString("polarisation");
if (coPolarization.contains(pol)) {
final MetadataElement noiseVectorListElem = noiElem.getElement("noiseVectorList");
final int count = Integer.parseInt(noiseVectorListElem.getAttributeString("count"));
Sentinel1Utils.NoiseVector[] noiseVectorList = Sentinel1Utils.getNoiseVector(noiseVectorListElem);
noiseVector = noiseVectorList[count / 2];
break;
}
}
if (noiseVector == null) {
throw new OperatorException("Input product does not have noise vector for HH or VV band");
}
}
/**
* Compute noise LUTs for pixels of a whole range line for given noise vector.
*/
private void computeNoiseLUT() {
try {
noiseLUT = new double[sourceImageWidth];
if (!thermalNoiseCorrectionPerformed) {
int pixelIdx = getPixelIndex(0, noiseVector);
final int maxLength = noiseVector.pixels.length - 2;
for (int x = 0; x < sourceImageWidth; x++) {
if (x > noiseVector.pixels[pixelIdx + 1] && pixelIdx < maxLength) {
pixelIdx++;
}
final int xx0 = noiseVector.pixels[pixelIdx];
final int xx1 = noiseVector.pixels[pixelIdx + 1];
final double muX = (double) (x - xx0) / (double) (xx1 - xx0);
noiseLUT[x] = Maths.interpolationLinear(
noiseVector.noiseLUT[pixelIdx], noiseVector.noiseLUT[pixelIdx + 1], muX) * scalingFactor;
}
} else {
for (int x = 0; x < sourceImageWidth; x++) {
noiseLUT[x] = 0.0;
}
}
} catch (Throwable e) {
OperatorUtils.catchOperatorException("computeNoiseLUT", e);
}
}
/**
* Get pixel index in a given noise vector for a given pixel.
*
* @param x Pixel coordinate.
* @param noiseVector Noise vector.
* @return The pixel index.
*/
private static int getPixelIndex(final int x, final Sentinel1Utils.NoiseVector noiseVector) {
for (int i = 0; i < noiseVector.pixels.length; i++) {
if (x < noiseVector.pixels[i]) {
return i - 1;
}
}
return noiseVector.pixels.length - 2;
}
/**
* Compute noise scaling factor.
*/
private void computeNoiseScalingFactor() throws IOException {
if (version < 2.50) {
final String acquisitionMode = absRoot.getAttributeString(AbstractMetadata.ACQUISITION_MODE);
double knoise;
if (acquisitionMode.contains("IW")) {
knoise = 75088.7;
} else if (acquisitionMode.contains("EW")) {
knoise = 56065.87;
} else {
throw new OperatorException("Cannot apply the operator to the input GRD product.");
}
final double dn0 = getDN0();
if (version < 2.34) {
scalingFactor = knoise * dn0;
} else {
scalingFactor = knoise * dn0 * dn0;
}
} else {
scalingFactor = 1.0;
}
}
/**
* Get the first element in DN vector from the original product metadata.
*/
private double getDN0() throws IOException {
String[] selectedPols = {coPolarization};
Sentinel1Calibrator.CalibrationInfo[] calibration = Sentinel1Calibrator.getCalibrationVectors(
sourceProduct, Arrays.asList(selectedPols), false, false, false, true);
for (Sentinel1Calibrator.CalibrationInfo cal : calibration) {
if (cal.polarization.contains(coPolarization)) {
final float[] dnLUT = Sentinel1Calibrator.getVector(
Sentinel1Calibrator.CALTYPE.DN, cal.getCalibrationVector(0));
return dnLUT[0];
}
}
return 0.0;
}
/**
* Create a target product for output.
*/
private void createTargetProduct() {
targetProduct = new Product(sourceProduct.getName(),
sourceProduct.getProductType(),
sourceImageWidth,
sourceImageHeight);
addSelectedBands();
ProductUtils.copyProductNodes(sourceProduct, targetProduct);
}
/**
* Add user selected bands to target product.
*/
private void addSelectedBands() {
final Band[] sourceBands = sourceProduct.getBands();
for (Band srcBand : sourceBands) {
final String unit = srcBand.getUnit();
if (unit == null) {
throw new OperatorException("band " + srcBand.getName() + " requires a unit");
}
if (!unit.contains(Unit.AMPLITUDE) && !unit.contains(Unit.INTENSITY)) {
continue;
}
final String srcBandName = srcBand.getName();
if (selectedPolarisations != null && selectedPolarisations.length != 0 &&
!containSelectedPolarisations(srcBandName)) {
continue;
}
if (srcBand instanceof VirtualBand) {
ProductUtils.copyVirtualBand(targetProduct, (VirtualBand) srcBand, srcBand.getName());
} else {
final Band targetBand = new Band(
srcBandName,
srcBand.getDataType(),
srcBand.getRasterWidth(),
srcBand.getRasterHeight());
targetBand.setUnit(srcBand.getUnit());
targetBand.setNoDataValue(srcBand.getNoDataValue());
targetBand.setNoDataValueUsed(srcBand.isNoDataValueUsed());
targetBand.setDescription(srcBand.getDescription());
targetProduct.addBand(targetBand);
}
}
}
private boolean containSelectedPolarisations(final String bandName) {
for (String pol : selectedPolarisations) {
if (bandName.contains(pol)) {
return true;
}
}
return false;
}
/**
* 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 {
try {
if (useBorderDetection && !borderDetected) {
detectBorders();
}
final int x0 = targetRectangle.x;
final int y0 = targetRectangle.y;
final int w = targetRectangle.width;
final int h = targetRectangle.height;
final int yMax = y0 + h;
final int xMax = x0 + w;
//System.out.println("x0 = " + x0 + ", y0 = " + y0 + ", w = " + w + ", h = " + h);
final Set<Band> keySet = targetTileMap.keySet();
final int numBands = keySet.size();
final ProductData[] targetData = new ProductData[numBands];
final ProductData[] sourceData = new ProductData[numBands];
final Tile[] sourceTile = new Tile[numBands];
final Tile[] targetTile = new Tile[numBands];
final double[] bandNoDataValues = new double[numBands];
int k = 0;
for (Band targetBand : keySet) {
targetTile[k] = targetTileMap.get(targetBand);
targetData[k] = targetTile[k].getDataBuffer();
final Band srcBand = sourceProduct.getBand(targetBand.getName());
sourceTile[k] = getSourceTile(srcBand, targetRectangle);
sourceData[k] = sourceTile[k].getDataBuffer();
bandNoDataValues[k] = srcBand.getNoDataValue();
k++;
}
final TileIndex srcIndex = new TileIndex(sourceTile[0]);
final TileIndex tgtIndex = new TileIndex(targetTile[0]);
final Tile coPolTile = getSourceTile(coPolBand, targetRectangle);
final ProductData coPolData = coPolTile.getDataBuffer();
double coPolDataValue, deNoisedDataValue;
for (int y = y0; y < yMax; y++) {
srcIndex.calculateStride(y);
tgtIndex.calculateStride(y);
for (int x = x0; x < xMax; x++) {
final int srcIdx = srcIndex.getIndex(x);
boolean testPixel = x < leftBorder || x > rightBorder || y < topBorder || y > bottomBorder;
if (testPixel) {
coPolDataValue = coPolData.getElemDoubleAt(srcIdx);
if (noDataValue.equals(coPolDataValue)) {
continue;
}
deNoisedDataValue =
Math.sqrt(Math.max(coPolDataValue * coPolDataValue - noiseLUT[x], 0.0));
if (deNoisedDataValue < trimThreshold || coPolDataValue < 30) {
final int tgtIdx = tgtIndex.getIndex(x);
for (int i = 0; i < numBands; i++) {
targetData[i].setElemDoubleAt(tgtIdx, bandNoDataValues[i]);
}
} else {
testPixel = false;
}
}
if (!testPixel) {
final int tgtIdx = tgtIndex.getIndex(x);
for (int i = 0; i < numBands; i++) {
targetData[i].setElemDoubleAt(tgtIdx, sourceData[i].getElemDoubleAt(srcIdx));
}
}
}
}
} catch (Throwable e) {
throw new OperatorException(e.getMessage());
}
}
private synchronized void detectBorders() throws OperatorException {
if (borderDetected) return;
detectBorder();
SystemUtils.LOG.fine("topBorder = " + topBorder);
SystemUtils.LOG.fine("bottomBorder = " + bottomBorder);
SystemUtils.LOG.fine("leftBorder = " + leftBorder);
SystemUtils.LOG.fine("rightBorder = " + rightBorder);
borderDetected = true;
}
private enum SIDE {TOP, BOTTOM, LEFT, RIGHT}
private void detectBorder() throws OperatorException {
final Dimension tileSize = new Dimension(borderLimit, borderLimit);
final Rectangle[] topRectangles = getAllTileRectangles(0, 0, sourceImageWidth, borderLimit, tileSize);
final Rectangle[] bottomRectangles = getAllTileRectangles(
0, sourceImageHeight - borderLimit, sourceImageWidth, borderLimit, tileSize);
final Rectangle[] leftRectangles = getAllTileRectangles(0, 0, borderLimit, sourceImageHeight, tileSize);
final Rectangle[] rightRectangles = getAllTileRectangles(
sourceImageWidth - borderLimit, 0, borderLimit, sourceImageHeight, tileSize);
final List<Border> borderList = new ArrayList<>(4);
borderList.add(new Border(SIDE.TOP, borderLimit, topRectangles));
borderList.add(new Border(SIDE.BOTTOM, borderLimit, bottomRectangles));
borderList.add(new Border(SIDE.LEFT, borderLimit, leftRectangles));
borderList.add(new Border(SIDE.RIGHT, borderLimit, rightRectangles));
final int totalRects = topRectangles.length + bottomRectangles.length + leftRectangles.length + rightRectangles.length;
final ThreadManager threadManager = new ThreadManager();
final StatusProgressMonitor status = new StatusProgressMonitor(StatusProgressMonitor.TYPE.SUBTASK);
status.beginTask("Detecting border... ", totalRects);
try {
for (Border border : borderList) {
for (final Rectangle rectangle : border.tileRectangles) {
final Thread worker = new Thread() {
final int xMax = rectangle.x + rectangle.width;
final int yMax = rectangle.y + rectangle.height;
@Override
public void run() {
final Tile coPolTile = getSourceTile(coPolBand, rectangle);
final ProductData coPolData = coPolTile.getDataBuffer();
final TileIndex srcIndex = new TileIndex(coPolTile);
if (border.side.equals(SIDE.TOP) || border.side.equals(SIDE.BOTTOM)) {
final double[] colSum = new double[borderLimit];
for (int y = rectangle.y; y < yMax; ++y) {
srcIndex.calculateStride(y);
final int yy = y - rectangle.y;
double sum = 0.0;
for (int x = rectangle.x; x < xMax; ++x) {
final double v = coPolData.getElemDoubleAt(srcIndex.getIndex(x));
if (v != noDataValue) {
sum += v;
}
}
colSum[yy] += sum;
}
synchronized (border.avg) {
for (int r = 0; r < borderLimit; r++) {
border.avg[r] += (colSum[r] / rectangle.width);
}
}
} else {
final double[] rowSum = new double[borderLimit];
for (int x = rectangle.x; x < xMax; ++x) {
final int xx = x - rectangle.x;
double sum = 0.0;
for (int y = rectangle.y; y < yMax; ++y) {
final double v = coPolData.getElemDoubleAt(coPolTile.getDataBufferIndex(x, y));
if (v != noDataValue) {
sum += v;
}
}
rowSum[xx] += sum;
}
synchronized (border.avg) {
for (int c = 0; c < borderLimit; c++) {
border.avg[c] += (rowSum[c] / rectangle.height);
}
}
}
}
};
threadManager.add(worker);
status.worked(1);
}
}
threadManager.finish();
} catch (Exception e) {
OperatorUtils.catchOperatorException(getId() + " detectBorder ", e);
} finally {
status.done();
}
for (Border border : borderList) {
for (int r = 0; r < borderLimit; r++) {
border.avg[r] /= border.tileRectangles.length;
}
final int peakPos = getPeakPosition(border.avg);
switch (border.side) {
case TOP: {
if (peakPos == -1) {
topBorder = 0;
} else {
topBorder = peakPos;
}
break;
}
case BOTTOM: {
if (peakPos == -1) {
bottomBorder = sourceImageHeight - 1;
} else {
bottomBorder = sourceImageHeight - borderLimit + peakPos;
}
break;
}
case LEFT: {
if (peakPos == -1) {
leftBorder = 0;
} else {
leftBorder = peakPos;
}
break;
}
case RIGHT: {
if (peakPos == -1) {
rightBorder = sourceImageWidth - 1;
} else {
rightBorder = sourceImageWidth - borderLimit + (int) (peakPos * 0.8);
}
break;
}
}
}
}
private static class Border {
public final SIDE side;
public final double[] avg;
public final Rectangle[] tileRectangles;
public Border(final SIDE side, final int borderLimit, final Rectangle[] rectangles) {
this.side = side;
this.avg = new double[borderLimit];
this.tileRectangles = rectangles;
}
}
public static Rectangle[] getAllTileRectangles(
final int x0, final int y0, final int width, final int height, final Dimension tileSize) {
final Rectangle boundary = new Rectangle(x0, y0, width, height);
final int tileCountX = MathUtils.ceilInt(boundary.width / (double) tileSize.width);
final int tileCountY = MathUtils.ceilInt(boundary.height / (double) tileSize.height);
final Rectangle[] rectangles = new Rectangle[tileCountX * tileCountY];
int index = 0;
for (int y = y0; y < y0 + height; y += tileSize.height) {
for (int x = x0; x < x0 + width; x += tileSize.width) {
final Rectangle tileRectangle = new Rectangle(x, y, tileSize.width, tileSize.height);
final Rectangle intersection = boundary.intersection(tileRectangle);
rectangles[index] = intersection;
index++;
}
}
return rectangles;
}
private static int getPeakPosition(final double[] array) {
final double[] derivative = new double[array.length - 1];
for (int i = 0; i < derivative.length; i++) {
derivative[i] = Math.abs(array[i + 1] - array[i]);
}
int maxIdx = -1;
double max = -Double.MAX_VALUE;
double mean = 0.0;
for (int i = 0; i < derivative.length; i++) {
mean += derivative[i];
if (max < derivative[i]) {
max = derivative[i];
maxIdx = i;
}
}
mean /= derivative.length;
if (mean > 0.0 && max / mean > 10.0) {
return maxIdx;
} else {
return -1;
}
}
/**
* 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(RemoveGRDBorderNoiseOp.class);
}
}
}