Examples with InterpolationCase org.geotoolkit.image.interpolation.InterpolationCase used on opensource projects

Example 1 with InterpolationCase

use of org.geotoolkit.image.interpolation.InterpolationCase in project geotoolkit by Geomatys.

the class ComputeVolumeProcess method execute.

/**
 * {@inheritDoc }.
 */
@Override
protected void execute() throws ProcessException {
    ArgumentChecks.ensureNonNull("inputParameters", inputParameters);
    final GridCoverageResource gcReader = inputParameters.getValue(ComputeVolumeDescriptor.IN_GRIDCOVERAGE_READER);
    final Geometry jtsGeom = inputParameters.getValue(ComputeVolumeDescriptor.IN_JTSGEOMETRY);
    CoordinateReferenceSystem geomCRS = inputParameters.getValue(ComputeVolumeDescriptor.IN_GEOMETRY_CRS);
    final Integer bIndex = inputParameters.getValue(ComputeVolumeDescriptor.IN_INDEX_BAND);
    final Double zMinCeil = inputParameters.getValue(ComputeVolumeDescriptor.IN_GEOMETRY_ALTITUDE);
    final double zMaxCeiling = inputParameters.getValue(ComputeVolumeDescriptor.IN_MAX_ALTITUDE_CEILING);
    final int bandIndex = (bIndex == null) ? 0 : (int) bIndex;
    final double zGroundCeiling = (zMinCeil == null) ? 0 : (double) zMinCeil;
    final boolean positiveSens = zGroundCeiling < zMaxCeiling;
    if (zGroundCeiling == zMaxCeiling) {
        outputParameters.getOrCreate(ComputeVolumeDescriptor.OUT_VOLUME_RESULT).setValue(0);
        return;
    }
    try {
        /*
             * geomCRS attribut should be null, we looking for find another way to define geometry CoordinateReferenceSystem.
             * It may be already stipulate in JTS geometry.
             */
        if (geomCRS == null) {
            geomCRS = JTS.findCoordinateReferenceSystem(jtsGeom);
        }
        final GridGeometry covGridGeom = gcReader.getGridGeometry();
        /*
             * If we have no CRS informations from geometry we consider that geometry is defined in same crs as Coverage.
             */
        final CoordinateReferenceSystem covCrs = covGridGeom.getCoordinateReferenceSystem();
        if (geomCRS == null) {
            geomCRS = covCrs;
        }
        final MathTransform covToGeomCRS = CRS.findOperation(covCrs, geomCRS, null).getMathTransform();
        // -- next read only interest area.
        final Envelope envGeom = jtsGeom.getEnvelopeInternal();
        final Envelope2D envGeom2D = new Envelope2D(geomCRS, envGeom.getMinX(), envGeom.getMinY(), envGeom.getWidth(), envGeom.getHeight());
        /**
         ****************************************
         */
        final GridCoverage dem = gcReader.read(gcReader.getGridGeometry().derive().subgrid(envGeom2D).build());
        final SampleDimension gsd = dem.getSampleDimensions().get(bandIndex);
        final MathTransform1D zmt = gsd.getTransferFunction().orElse(null);
        if (zmt == null) {
            throw new ProcessException("you should stipulate MathTransform1D from sampleDimension to geophysic.", this, null);
        }
        final GridGeometry gg2d = dem.getGridGeometry();
        InterpolationCase interpolationChoice;
        // -- adapt interpolation in function of grid extend
        final GridExtent gridEnv2D = gg2d.getExtent();
        final int[] subSpace = gridEnv2D.getSubspaceDimensions(2);
        final long gWidth = gridEnv2D.getSize(subSpace[0]);
        final long gHeight = gridEnv2D.getSize(subSpace[1]);
        if (gWidth < 1 || gHeight < 1) {
            outputParameters.getOrCreate(ComputeVolumeDescriptor.OUT_VOLUME_RESULT).setValue(0);
            return;
        } else if (gWidth < 2 || gHeight < 2) {
            interpolationChoice = InterpolationCase.NEIGHBOR;
        } else if (gWidth < 4 || gHeight < 4) {
            interpolationChoice = InterpolationCase.BILINEAR;
        } else {
            // -- paranoiac assert
            assert gWidth >= 4 && gHeight >= 4;
            interpolationChoice = InterpolationCase.BICUBIC;
        }
        final MathTransform gridToCrs = gg2d.getGridToCRS(PixelInCell.CELL_CENTER);
        final CoordinateSystem destCS = covCrs.getCoordinateSystem();
        final RenderedImage mnt = dem.render(null);
        final Interpolation interpol = Interpolation.create(new PixelIterator.Builder().setIteratorOrder(SequenceType.LINEAR).create(mnt), interpolationChoice, 0, ResampleBorderComportement.EXTRAPOLATION, null);
        final MathTransform gridToGeom = MathTransforms.concatenate(gridToCrs, covToGeomCRS);
        final StepPixelAreaCalculator stePixCalculator;
        if (covCrs instanceof GeographicCRS) {
            stePixCalculator = new GeographicStepPixelAreaCalculator(PIXELSTEP, covCrs, gridToCrs);
        } else {
            if (destCS instanceof CartesianCS) {
                // -- resolution
                final double[] resolution;
                try {
                    resolution = gg2d.getResolution(false);
                } catch (IncompleteGridGeometryException ex) {
                    throw new ProcessException("Cannot estimate resolution", this, ex);
                }
                final int dimDestCS = destCS.getDimension();
                final int destDim = destCS.getDimension();
                final UnitConverter[] unitConverters = new UnitConverter[dimDestCS];
                for (int d = 0; d < destDim; d++) {
                    final CoordinateSystemAxis csA = destCS.getAxis(d);
                    unitConverters[d] = csA.getUnit().getConverterToAny(METER);
                }
                // -- pixel step computing in m²
                stePixCalculator = new CartesianStepPixelAreaCalculator(PIXELSTEP, unitConverters, resolution);
            } else {
                throw new ProcessException("Coordinate reference system configuration not supported. CRS should be instance of geographic crs or has a cartesian coordinate system.", this, null);
            }
        }
        // -- geometry factory to create point at n step to test if it is within geometry
        final GeometryFactory gf = JTS.getFactory();
        // -- coordinate to test if point is within geom
        final Coordinate coords = new Coordinate();
        // -- image attributs
        final double minx = mnt.getMinX() - 0.5;
        final double miny = mnt.getMinY() - 0.5;
        final double maxx = minx + mnt.getWidth();
        final double maxy = miny + mnt.getHeight();
        final double debx = minx + PIXELSTEP / 2.0;
        final double[] pixPoint = new double[] { debx, miny + PIXELSTEP / 2.0 };
        final double[] geomPoint = new double[2];
        double volume = 0;
        final UnitConverter hconverter;
        if (!gsd.getUnits().isPresent() || Units.UNITY.equals(gsd.getUnits().get())) {
            // -- unit unknowed, assume it's meters already
            hconverter = METER.getConverterTo(METER);
        } else {
            hconverter = gsd.getUnits().get().getConverterToAny(METER);
        }
        while (pixPoint[1] < maxy) {
            stopIfDismissed();
            pixPoint[0] = debx;
            while (pixPoint[0] < maxx) {
                stopIfDismissed();
                // -- project point in geomtry CRS
                gridToGeom.transform(pixPoint, 0, geomPoint, 0, 1);
                // -- test if point is within geometry.
                coords.setOrdinate(0, geomPoint[0]);
                coords.setOrdinate(1, geomPoint[1]);
                if (jtsGeom.contains(gf.createPoint(coords))) {
                    // -- get interpolate value
                    double h = interpol.interpolate(pixPoint[0], pixPoint[1], bandIndex);
                    // -- projet h in geophysic value
                    h = zmt.transform(h);
                    // -- convert in meter
                    h = hconverter.convert(h);
                    // -- Verify that h value found is in appropriate interval.
                    if ((positiveSens && h > zGroundCeiling) || (!positiveSens && h < zGroundCeiling)) {
                        // -- add in volum
                        volume += (Math.min(Math.abs(h - zGroundCeiling), Math.abs(zMaxCeiling - zGroundCeiling))) * stePixCalculator.computeStepPixelArea(pixPoint);
                    }
                }
                pixPoint[0] += PIXELSTEP;
            }
            pixPoint[1] += PIXELSTEP;
        }
        outputParameters.getOrCreate(ComputeVolumeDescriptor.OUT_VOLUME_RESULT).setValue(volume);
    } catch (Exception ex) {
        throw new ProcessException(ex.getMessage(), this, ex);
    }
}

Example 2 with InterpolationCase

use of org.geotoolkit.image.interpolation.InterpolationCase in project geotoolkit by Geomatys.

the class AbstractCoverageSymbolizerRenderer method getObjectiveCoverage.

/**
 * Returns expected {@linkplain GridCoverage elevation coverage} or {@linkplain GridCoverage coverage}
 * from given {@link ProjectedCoverage}.
 *
 * TODO: add a margin or interpolation parameter. To properly interpolate border on output canvas, we need extra
 *  lines/columns on source image. Their number depends on applied interpolation (bilinear, bicubic, etc.).
 *
 * @param ref coverage resource
 * @param canvasGrid Rendering canvas grid geometry
 * @param isElevation {@code true} if we want elevation coverage, else ({@code false}) for features coverage.
 * @param sourceBands coverage source bands to features
 * @return expected {@linkplain GridCoverage elevation coverage} or {@linkplain GridCoverage coverage}
 * @throws org.geotoolkit.coverage.io.CoverageStoreException if problem during coverage reading.
 * @throws org.opengis.referencing.operation.TransformException if problem during {@link Envelope} transformation.
 * @throws org.opengis.util.FactoryException if problem during {@link Envelope} study.
 * @throws org.geotoolkit.process.ProcessException if problem during resampling processing.
 * @see ProjectedCoverage#getElevationCoverage(org.geotoolkit.coverage.io.GridCoverageReadParam)
 * @see ProjectedCoverage#getCoverage(org.geotoolkit.coverage.io.GridCoverageReadParam)
 */
protected GridCoverage getObjectiveCoverage(final GridCoverageResource ref, GridGeometry canvasGrid, final boolean isElevation, int[] sourceBands) throws DataStoreException, TransformException, FactoryException, ProcessException {
    ArgumentChecks.ensureNonNull("projectedCoverage", ref);
    final InterpolationCase interpolation = InterpolationCase.BILINEAR;
    final GridGeometry refGG = ref.getGridGeometry();
    // fast envelope intersection in 2D
    if (refGG.isDefined(GridGeometry.ENVELOPE)) {
        Envelope bbox = renderingContext.getCanvasObjectiveBounds2D();
        Envelope refEnv = Envelopes.transform(refGG.getEnvelope(), bbox.getCoordinateReferenceSystem());
        if (!AbstractEnvelope.castOrCopy(bbox).intersects(refEnv, true)) {
            throw new DisjointExtentException("Coverage resource envelope do not intersect canvas");
        }
    }
    final GridGeometry baseGG = trySubGrid(refGG, canvasGrid);
    final GridGeometry slice = extractSlice(baseGG, canvasGrid, computeMargin2D(interpolation), true);
    if (sourceBands != null && sourceBands.length < 1)
        sourceBands = null;
    GridCoverage coverage = ref.read(slice, sourceBands);
    if (coverage instanceof GridCoverageStack) {
        Logger.getLogger("org.geotoolkit.display2d.primitive").log(Level.WARNING, "Coverage reader return more than one slice.");
    }
    while (coverage instanceof GridCoverageStack) {
        // pick the first slice
        coverage = ((GridCoverageStack) coverage).coverageAtIndex(0);
    }
    // we remove all other dimension to simplify any following operation
    if (coverage.getCoordinateReferenceSystem().getCoordinateSystem().getDimension() > 2) {
        coverage = new ReducedGridCoverage(coverage, 0, 1);
    }
    final CoordinateReferenceSystem crs2d = CRS.getHorizontalComponent(canvasGrid.getCoordinateReferenceSystem());
    if (Utilities.equalsIgnoreMetadata(crs2d, coverage.getCoordinateReferenceSystem())) {
        return coverage;
    } else {
        coverage = prepareCoverageToResampling(coverage, symbol);
        // resample
        final double[] fill = new double[coverage.getSampleDimensions().size()];
        Arrays.fill(fill, Double.NaN);
        // ///// HACK FOR 0/360 /////////////////////////////////////////
        GeneralEnvelope ge = new GeneralEnvelope(coverage.getGridGeometry().getEnvelope());
        try {
            GeneralEnvelope cdt = GeneralEnvelope.castOrCopy(Envelopes.transform(coverage.getGridGeometry().getEnvelope(), CommonCRS.WGS84.normalizedGeographic()));
            cdt.normalize();
            if (!cdt.isEmpty()) {
                ge = cdt;
            }
        } catch (ProjectionException ex) {
            LOGGER.log(Level.INFO, ex.getMessage(), ex);
        }
        GridGeometry resampleGrid = canvasGrid;
        try {
            resampleGrid = resampleGrid.derive().rounding(GridRoundingMode.ENCLOSING).subgrid(ge).build().reduce(0, 1);
        } catch (DisjointExtentException ex) {
        // don't log, still continue
        } catch (IllegalGridGeometryException ex) {
            LOGGER.log(Level.INFO, ex.getMessage(), ex);
        }
        resampleGrid = CoverageUtilities.forceLowerToZero(resampleGrid);
        // ///// HACK FOR 0/360 /////////////////////////////////////////
        final MathTransform gridToCRS = coverage.getGridGeometry().getGridToCRS(PixelInCell.CELL_CENTER);
        if (isNonLinear(gridToCRS)) {
            final GridGeometry slice2 = extractSlice(refGG, canvasGrid, computeMargin2D(interpolation), false);
            coverage = ref.read(slice2, sourceBands);
            // we remove all other dimension to simplify any following operation
            if (coverage.getCoordinateReferenceSystem().getCoordinateSystem().getDimension() > 2) {
                coverage = new ReducedGridCoverage(coverage, 0, 1);
            }
            return forwardResample(coverage, resampleGrid);
        } else {
            ResampleProcess process = new ResampleProcess(coverage, crs2d, resampleGrid, interpolation, fill);
            // do not extrapolate values, can cause large areas of incorrect values
            process.getInput().parameter(ResampleDescriptor.IN_BORDER_COMPORTEMENT_TYPE.getName().getCode()).setValue(ResampleBorderComportement.FILL_VALUE);
            return process.executeNow();
        }
    }
}

Example 3 with InterpolationCase

use of org.geotoolkit.image.interpolation.InterpolationCase in project geotoolkit by Geomatys.

the class ResampleProcess method execute.

/**
 * Resamples a grid coverage.
 */
@Override
protected void execute() throws ProcessException {
    final GridCoverage source = inputParameters.getValue(IN_COVERAGE);
    final double[] background = inputParameters.getValue(IN_BACKGROUND);
    InterpolationCase interpolation = inputParameters.getValue(IN_INTERPOLATION_TYPE);
    final ResampleBorderComportement border = inputParameters.getValue(IN_BORDER_COMPORTEMENT_TYPE);
    CoordinateReferenceSystem targetCRS = (CoordinateReferenceSystem) inputParameters.parameter("CoordinateReferenceSystem").getValue();
    final GridGeometry targetGG = inputParameters.getValue(IN_GRID_GEOMETRY);
    final GridCoverage target;
    try {
        target = reproject(source, targetCRS, targetGG, interpolation, border, background);
    } catch (FactoryException exception) {
        throw new CannotReprojectException(Errors.format(Errors.Keys.CantReprojectCoverage_1, CoverageUtilities.getName(source)), exception);
    } catch (TransformException exception) {
        throw new CannotReprojectException(Errors.format(Errors.Keys.CantReprojectCoverage_1, CoverageUtilities.getName(source)), exception);
    }
    outputParameters.getOrCreate(OUT_COVERAGE).setValue(target);
}

Example 4 with InterpolationCase

use of org.geotoolkit.image.interpolation.InterpolationCase in project geotoolkit by Geomatys.

the class RasterPresentation method renderCoverage.

private boolean renderCoverage(RenderingContext2D renderingContext, org.apache.sis.coverage.grid.GridCoverage coverage, MathTransform trs2D) throws PortrayalException {
    final Graphics2D g2d = renderingContext.getGraphics();
    final CanvasMonitor monitor = renderingContext.getMonitor();
    boolean dataRendered = false;
    RenderedImage img = coverage.render(null);
    /*
         * Try to prefetch image before rendering
         * resampled image or mosaic have deferred tiles
         * java2d render tiles one by one which can be slow when working with
         * computed coverages or distant services like WMTS or TMS
         */
    if ((img.getWidth() * img.getHeight()) < 5000 * 5000) {
        ImageProcessor processor = new ImageProcessor();
        processor.setExecutionMode(ImageProcessor.Mode.PARALLEL);
        img = processor.prefetch(img, null);
    }
    final InterpolationCase interpolationCase = (InterpolationCase) renderingContext.getRenderingHints().get(GO2Hints.KEY_INTERPOLATION);
    if (interpolationCase != null) {
        switch(interpolationCase) {
            case NEIGHBOR:
                g2d.setRenderingHint(RenderingHints.KEY_INTERPOLATION, RenderingHints.VALUE_INTERPOLATION_NEAREST_NEIGHBOR);
                break;
            case BILINEAR:
                g2d.setRenderingHint(RenderingHints.KEY_INTERPOLATION, RenderingHints.VALUE_INTERPOLATION_BILINEAR);
                break;
            case BICUBIC:
                g2d.setRenderingHint(RenderingHints.KEY_INTERPOLATION, RenderingHints.VALUE_INTERPOLATION_BICUBIC);
                break;
            default:
                // resample image ourself
                try {
                    GridCoverage cov = new ResampleProcess(coverage, renderingContext.getGridGeometry().getCoordinateReferenceSystem(), renderingContext.getGridGeometry2D(), interpolationCase, null).executeNow();
                    trs2D = renderingTransform(cov.getGridGeometry());
                    img = cov.render(null);
                } catch (ProcessException ex) {
                    throw new PortrayalException(ex);
                }
                break;
        }
    }
    if (trs2D instanceof AffineTransform) {
        try {
            g2d.drawRenderedImage(RenderingWorkaround.wrap(img), (AffineTransform) trs2D);
            dataRendered = true;
        } catch (Exception ex) {
            final StringWriter sw = new StringWriter();
            final PrintWriter pw = new PrintWriter(sw);
            ex.printStackTrace(pw);
            if (ex instanceof ArrayIndexOutOfBoundsException) {
                // we can recover when it's an inapropriate componentcolormodel
                final StackTraceElement[] eles = ex.getStackTrace();
                if (eles.length > 0 && ComponentColorModel.class.getName().equalsIgnoreCase(eles[0].getClassName())) {
                    try {
                        final Map<String, Object> analyze = StatisticOp.analyze(img);
                        final double[] minArray = (double[]) analyze.get(StatisticOp.MINIMUM);
                        final double[] maxArray = (double[]) analyze.get(StatisticOp.MAXIMUM);
                        final double min = findExtremum(minArray, true);
                        final double max = findExtremum(maxArray, false);
                        final List<InterpolationPoint> values = new ArrayList<>();
                        values.add(new DefaultInterpolationPoint(Double.NaN, GO2Utilities.STYLE_FACTORY.literal(new Color(0, 0, 0, 0))));
                        values.add(new DefaultInterpolationPoint(min, GO2Utilities.STYLE_FACTORY.literal(Color.BLACK)));
                        values.add(new DefaultInterpolationPoint(max, GO2Utilities.STYLE_FACTORY.literal(Color.WHITE)));
                        final Literal lookup = StyleConstants.DEFAULT_CATEGORIZE_LOOKUP;
                        final Literal fallback = StyleConstants.DEFAULT_FALLBACK;
                        final Expression function = GO2Utilities.STYLE_FACTORY.interpolateFunction(lookup, values, Method.COLOR, Mode.LINEAR, fallback);
                        final CompatibleColorModel model = new CompatibleColorModel(img.getColorModel().getPixelSize(), function);
                        final ImageLayout layout = new ImageLayout().setColorModel(model);
                        img = new NullOpImage(img, layout, null, OpImage.OP_COMPUTE_BOUND);
                        g2d.drawRenderedImage(RenderingWorkaround.wrap(img), (AffineTransform) trs2D);
                        dataRendered = true;
                    } catch (Exception e) {
                        // plenty of errors can happen when painting an image
                        monitor.exceptionOccured(e, Level.WARNING);
                        // raise the original error
                        monitor.exceptionOccured(ex, Level.WARNING);
                    }
                } else {
                    // plenty of errors can happen when painting an image
                    monitor.exceptionOccured(ex, Level.WARNING);
                }
            } else {
                // plenty of errors can happen when painting an image
                monitor.exceptionOccured(ex, Level.WARNING);
            }
        }
    } else if (trs2D instanceof LinearTransform) {
        final LinearTransform lt = (LinearTransform) trs2D;
        final int col = lt.getMatrix().getNumCol();
        final int row = lt.getMatrix().getNumRow();
        // TODO using only the first parameters of the linear transform
        throw new PortrayalException("Could not render image, GridToCRS is a not an AffineTransform, found a " + trs2D.getClass());
    } else {
        throw new PortrayalException("Could not render image, GridToCRS is a not an AffineTransform, found a " + trs2D.getClass());
    }
    return dataRendered;
}