This is an automated email from the ASF dual-hosted git repository. desruisseaux pushed a commit to branch geoapi-4.0 in repository https://gitbox.apache.org/repos/asf/sis.git
commit 411009be20b5789cf28c4717b4d6e0e02a67d44b Author: Martin Desruisseaux <[email protected]> AuthorDate: Thu Sep 16 16:07:34 2021 +0200 Modify Lanczos implementation for: - reducing the amount of `Math.sin(x)` computation - create only one temporary array instead of two - make unnecessary the call to `Arrays.fill(…)`. --- .../org/apache/sis/image/LanczosInterpolation.java | 91 +++++++++++++++------- 1 file changed, 61 insertions(+), 30 deletions(-) diff --git a/core/sis-feature/src/main/java/org/apache/sis/image/LanczosInterpolation.java b/core/sis-feature/src/main/java/org/apache/sis/image/LanczosInterpolation.java index 897788b..3a0bdca 100644 --- a/core/sis-feature/src/main/java/org/apache/sis/image/LanczosInterpolation.java +++ b/core/sis-feature/src/main/java/org/apache/sis/image/LanczosInterpolation.java @@ -16,7 +16,6 @@ */ package org.apache.sis.image; -import java.util.Arrays; import java.awt.Dimension; import java.nio.DoubleBuffer; @@ -87,27 +86,61 @@ final class LanczosInterpolation implements Interpolation { public void interpolate(final DoubleBuffer source, final int numBands, final double xfrac, final double yfrac, final double[] writeTo, final int writeToOffset) { - Arrays.fill(writeTo, writeToOffset, writeToOffset + numBands, 0); + /* + * Lanczos kernel coefficients, pre-computed because reused many times. We store them in a temporary array + * created in this method; we do not reuse a cached array because this method will be invoked concurrently + * in many threads. In the particular case of the last dimension (y in our case), the coefficients will be + * used only once, so we can compute them on-the-fly and avoid the creation of one array. + */ final double[] kx = new double[span]; - final double[] ky = new double[span]; - for (int i=0; i<span; i++) { - /* - * kernel(x) must be invoked with −a ≤ x ≤ +a. - * - * Given x = xfrac − (i+1 - a) - * with i ∈ [0 … 2a-1] - * we get x ∈ [xfrac − a … xfrac + a-1] - * - * The condition is met for xfrac ∈ [0 … 1]. - */ - final double offset = i+1 - a; - kx[i] = kernel(xfrac - offset); - ky[i] = kernel(yfrac - offset); + double ky; + /* + * The kernel values could be computed in two lines such as: + * + * x = (xfrac - (i+1 - a)) * PI; + * kx[i] = sin(x)/x * sin(x/a)/(x/a) for x ≠ 0 + * + * But we can take advantage that, for integer values of i, + * sin(x) will always have the same value except for the sign. + */ + double siny; + { // Block for keeping variables in local scope. + double sinx; + double xp = (xfrac - (1 - a)) * Math.PI; + double yp = (yfrac - (1 - a)) * Math.PI; + sinx = Math.sin(xp); + siny = Math.sin(yp); + kx[0] = kernel(xp, sinx); + ky = kernel(yp, siny); + for (int i=1; i<span;) { + /* + * kernel(x) must be invoked with −a ≤ x ≤ +a. + * + * Given x = xfrac − (i+1 - a) + * with i ∈ [0 … 2a-1] + * we get x ∈ [xfrac − a … xfrac + a-1] + * + * The condition is met for xfrac ∈ [0 … 1]. + */ + kx[i] = kernel((xfrac - (++i - a)) * Math.PI, sinx = -sinx); + } } + /* + * Now multiply pixel values by Lanczos kernel coefficients. + * The first loop over `i` is unrolled for using the existing `ky` value + * and for overwritten existing values in the destination array. + */ source.mark(); - for (int y=0; y<ky.length; y++) { - for (int x=0; x<kx.length; x++) { - final double k = kx[x] * ky[y]; + for (int i=0; i<span; i++) { + final double k = kx[i] * ky; + for (int b=0; b<numBands; b++) { + writeTo[writeToOffset + b] = k * source.get(); + } + } + for (int j=1; j<span;) { + ky = kernel((yfrac - (++j - a)) * Math.PI, siny = -siny); + for (int i=0; i<span; i++) { + final double k = kx[i] * ky; for (int b=0; b<numBands; b++) { writeTo[writeToOffset + b] += k * source.get(); } @@ -118,19 +151,17 @@ final class LanczosInterpolation implements Interpolation { /** * Computes a value of the Lanczos reconstruction kernel L(x). - * This is a component of Lanczos filter's kernel in two dimensions, - * which is L(x,y) = L(x)L(y). + * This is a component of Lanczos filter's kernel in two dimensions, which is L(x,y) = L(x)L(y). + * In this implementation, the given <var>x</var> value must be pre-multiplied by {@link Math#PI}. * - * @param x must be between −{@link #a} and +{@link #a} inclusive. + * <div class="note"><b>Note:</b> the multiplication by π is a normalization that causes + * the definite integral of the function over the real numbers to equal 1.</div> + * + * @param x must be between −{@link #a}⋅π and +{@link #a}⋅π inclusive. + * @param sinx value of {@code Math.sin(x)}. */ - private double kernel(double x) { - /* - * Formula is sin(x)/x * sin(x/a)/(x/a) rewritten below for one less division. - * The multiplication by π is a normalization that causes the definite integral - * of the function over the real numbers to equal 1. - */ - x *= Math.PI; - final double y = Math.sin(x) * Math.sin(x/a)*a / (x*x); + private double kernel(double x, final double sinx) { + final double y = sinx * Math.sin(x/a)*a / (x*x); return (y <= 1) ? y : 1; // Do not use Math.min(…) because we want to replace NaN by 1. } }
