| /* |
| * Licensed to the Apache Software Foundation (ASF) under one or more |
| * contributor license agreements. See the NOTICE file distributed with |
| * this work for additional information regarding copyright ownership. |
| * The ASF licenses this file to You under the Apache License, Version 2.0 |
| * (the "License"); you may not use this file except in compliance with |
| * the License. You may obtain a copy of the License at |
| * |
| * http://www.apache.org/licenses/LICENSE-2.0 |
| * |
| * Unless required by applicable law or agreed to in writing, software |
| * distributed under the License is distributed on an "AS IS" BASIS, |
| * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. |
| * See the License for the specific language governing permissions and |
| * limitations under the License. |
| */ |
| package org.apache.sis.math; |
| |
| import java.io.Serializable; |
| import java.util.Objects; |
| import java.util.Iterator; |
| import java.util.function.DoubleBinaryOperator; |
| import static java.lang.Math.abs; |
| import static java.lang.Math.sqrt; |
| import static java.lang.Math.ulp; |
| import org.opengis.geometry.DirectPosition; |
| import org.apache.sis.util.ArgumentChecks; |
| import org.apache.sis.util.privy.DoubleDouble; |
| import org.apache.sis.util.privy.Numerics; |
| import org.apache.sis.util.privy.Strings; |
| import org.apache.sis.util.resources.Errors; |
| |
| // Specific to the geoapi-3.1 and geoapi-4.0 branches: |
| import org.opengis.coordinate.MismatchedDimensionException; |
| |
| |
| /** |
| * Equation of a plane in a three-dimensional space (<var>x</var>,<var>y</var>,<var>z</var>). |
| * The plane equation is expressed by {@linkplain #slopeX() sx}, {@linkplain #slopeY() sy} and |
| * {@linkplain #z0() z₀} coefficients as below: |
| * |
| * <blockquote> |
| * <var>{@linkplain #z(double, double) z}</var>(<var>x</var>,<var>y</var>) = |
| * <var>sx</var>⋅<var>x</var> + <var>sy</var>⋅<var>y</var> + <var>z₀</var> |
| * </blockquote> |
| * |
| * Those coefficients can be set directly, or computed by a linear regression of this plane |
| * through a set of three-dimensional points. |
| * |
| * @author Martin Desruisseaux (MPO, IRD) |
| * @author Howard Freeland (MPO, for algorithmic inspiration) |
| * @version 1.4 |
| * |
| * @see Line |
| * @see org.apache.sis.referencing.operation.builder.LinearTransformBuilder |
| * |
| * @since 0.5 |
| */ |
| public class Plane implements DoubleBinaryOperator, Cloneable, Serializable { |
| /** |
| * Serial number for compatibility with different versions. |
| */ |
| private static final long serialVersionUID = 2956201711131316723L; |
| |
| /** |
| * Number of dimensions. |
| */ |
| private static final int DIMENSION = 3; |
| |
| /** |
| * Threshold value relative to 1 ULP of other terms in the z = sx⋅x + sy⋅y + z₀ equation. |
| * A value of 1 would be theoretically sufficient since adding a value smaller to 1 ULP to |
| * a {@code double} has no effect. Nevertheless we use a smaller value as a safety because: |
| * |
| * <ul> |
| * <li>We perform our checks using the points given by the user, but we don't know how |
| * representative those points are.</li> |
| * <li>We perform our checks using an approximation which may result in slightly different |
| * decisions (false positives) than what we would got by a more robust check.</li> |
| * </ul> |
| * |
| * This arbitrary threshold value may change in any future SIS version according experience gained. |
| * |
| * @see org.apache.sis.referencing.operation.matrix.GeneralMatrix#ZERO_THRESHOLD |
| * @see Numerics#COMPARISON_THRESHOLD |
| */ |
| private static final double ZERO_THRESHOLD = 1E-14; |
| |
| /** |
| * The slope along the <var>x</var> values. This coefficient appears in the plane equation |
| * <var><b><u>sx</u></b></var>⋅<var>x</var> + <var>sy</var>⋅<var>y</var> + <var>z₀</var>. |
| */ |
| private double sx; |
| |
| /** |
| * The slope along the <var>y</var> values. This coefficient appears in the plane equation |
| * <var>sx</var>⋅<var>x</var> + <var><b><u>sy</u></b></var>⋅<var>y</var> + <var>z₀</var>. |
| */ |
| private double sy; |
| |
| /** |
| * The <var>z</var> value at (<var>x</var>,<var>y</var>) = (0,0). This coefficient appears in the plane equation |
| * <var>sx</var>⋅<var>x</var> + <var>sy</var>⋅<var>y</var> + <b><u><var>z₀</var></u></b>. |
| */ |
| private double z0; |
| |
| /** |
| * Constructs a new plane with all coefficients initialized to {@link Double#NaN}. |
| */ |
| public Plane() { |
| sx = sy = z0 = Double.NaN; |
| } |
| |
| /** |
| * Constructs a new plane initialized to the given coefficients. |
| * |
| * @param sx the slope along the <var>x</var> values. |
| * @param sy the slope along the <var>y</var> values. |
| * @param z0 the <var>z</var> value at (<var>x</var>,<var>y</var>) = (0,0). |
| * |
| * @see #setEquation(double, double, double) |
| */ |
| public Plane(final double sx, final double sy, final double z0) { |
| this.sx = sx; |
| this.sy = sy; |
| this.z0 = z0; |
| } |
| |
| /** |
| * Returns the slope along the <var>x</var> values. This coefficient appears in the plane equation |
| * <var><b><u>sx</u></b></var>⋅<var>x</var> + <var>sy</var>⋅<var>y</var> + <var>z₀</var>. |
| * |
| * @return the <var>sx</var> term. |
| */ |
| public final double slopeX() { |
| return sx; |
| } |
| |
| /** |
| * Returns the slope along the <var>y</var> values. This coefficient appears in the plane equation |
| * <var>sx</var>⋅<var>x</var> + <var><b><u>sy</u></b></var>⋅<var>y</var> + <var>z₀</var>. |
| * |
| * @return the <var>sy</var> term. |
| */ |
| public final double slopeY() { |
| return sy; |
| } |
| |
| /** |
| * Returns the <var>z</var> value at (<var>x</var>,<var>y</var>) = (0,0). This coefficient appears in the |
| * plane equation <var>sx</var>⋅<var>x</var> + <var>sy</var>⋅<var>y</var> + <b><var>z₀</var></b>. |
| * |
| * @return the <var>z₀</var> term. |
| * |
| * @see #z(double, double) |
| */ |
| public final double z0() { |
| return z0; |
| } |
| |
| /** |
| * Computes the <var>x</var> value for the specified (<var>y</var>,<var>z</var>) point. |
| * The <var>x</var> value is computed using the following equation: |
| * |
| * <blockquote>x(y,z) = (z - ({@linkplain #z0() z₀} + {@linkplain #slopeY() sy}⋅y)) / {@linkplain #slopeX() sx}</blockquote> |
| * |
| * @param y the <var>y</var> value where to compute <var>x</var>. |
| * @param z the <var>z</var> value where to compute <var>x</var>. |
| * @return the <var>x</var> value. |
| */ |
| public final double x(final double y, final double z) { |
| return (z - (z0 + sy*y)) / sx; |
| } |
| |
| /** |
| * Computes the <var>y</var> value for the specified (<var>x</var>,<var>z</var>) point. |
| * The <var>y</var> value is computed using the following equation: |
| * |
| * <blockquote>y(x,z) = (z - ({@linkplain #z0() z₀} + {@linkplain #slopeX() sx}⋅x)) / {@linkplain #slopeY() sy}</blockquote> |
| * |
| * @param x the <var>x</var> value where to compute <var>y</var>. |
| * @param z the <var>z</var> value where to compute <var>y</var>. |
| * @return the <var>y</var> value. |
| */ |
| public final double y(final double x, final double z) { |
| return (z - (z0 + sx*x)) / sy; |
| } |
| |
| /** |
| * Computes the <var>z</var> value for the specified (<var>x</var>,<var>y</var>) point. |
| * The <var>z</var> value is computed using the following equation: |
| * |
| * <blockquote>z(x,y) = {@linkplain #slopeX() sx}⋅x + {@linkplain #slopeY() sy}⋅y + {@linkplain #z0() z₀}</blockquote> |
| * |
| * @param x the <var>x</var> value where to compute <var>z</var>. |
| * @param y the <var>y</var> value where to compute <var>z</var>. |
| * @return the <var>z</var> value. |
| * |
| * @see #z0() |
| */ |
| public final double z(final double x, final double y) { |
| return z0 + sx*x + sy*y; |
| } |
| |
| /** |
| * Evaluates this equation for the given values. The default implementation delegates to |
| * {@link #z(double,double) z(x,y)}, but subclasses may override with different formulas. |
| * This method is provided for interoperability with libraries making use of {@link java.util.function}. |
| * |
| * @param x the first operand where to evaluate the function. |
| * @param y the second operand where to evaluate the function. |
| * @return the function value for the given operands. |
| * |
| * @since 1.0 |
| */ |
| @Override |
| public double applyAsDouble(double x, double y) { |
| return z(x, y); |
| } |
| |
| /** |
| * Sets the equation of this plane to the given coefficients. |
| * |
| * @param sx the slope along the <var>x</var> values. |
| * @param sy the slope along the <var>y</var> values. |
| * @param z0 the <var>z</var> value at (<var>x</var>,<var>y</var>) = (0,0). |
| */ |
| public void setEquation(final double sx, final double sy, final double z0) { |
| this.sx = sx; |
| this.sy = sy; |
| this.z0 = z0; |
| } |
| |
| /** |
| * Sets this plane from values of arbitrary {@code Number} type. This method is invoked by algorithms that |
| * may produce other kind of numbers (for example with different precision) than the usual {@code double} |
| * primitive type. The default implementation delegates to {@link #setEquation(double, double, double)}, |
| * but subclasses can override this method if they want to process other kind of numbers in a special way. |
| * |
| * @param sx the slope along the <var>x</var> values. |
| * @param sy the slope along the <var>y</var> values. |
| * @param z0 the <var>z</var> value at (<var>x</var>,<var>y</var>) = (0,0). |
| * |
| * @since 0.8 |
| */ |
| public void setEquation(final Number sx, final Number sy, final Number z0) { |
| setEquation(sx.doubleValue(), sy.doubleValue(), z0.doubleValue()); |
| } |
| |
| /** |
| * Computes the plane's coefficients from the given coordinate values. |
| * This method uses a linear regression in the least-square sense, with the assumption that |
| * the (<var>x</var>,<var>y</var>) values are precise and all uncertainty is in <var>z</var>. |
| * {@link Double#NaN} values are ignored. The result is undetermined if all points are colinear. |
| * |
| * <p>The default implementation delegates to {@link #fit(Vector, Vector, Vector)}.</p> |
| * |
| * @param x vector of <var>x</var> coordinates. |
| * @param y vector of <var>y</var> coordinates. |
| * @param z vector of <var>z</var> values. |
| * @return an estimation of the Pearson correlation coefficient. |
| * The closer this coefficient is to +1 or -1, the better the fit. |
| * @throws IllegalArgumentException if <var>x</var>, <var>y</var> and <var>z</var> do not have the same length. |
| */ |
| public double fit(final double[] x, final double[] y, final double[] z) { |
| ArgumentChecks.ensureNonNull("x", x); |
| ArgumentChecks.ensureNonNull("y", y); |
| ArgumentChecks.ensureNonNull("z", z); |
| return fit(new ArrayVector.Doubles(x), new ArrayVector.Doubles(y), new ArrayVector.Doubles(z)); |
| } |
| |
| /** |
| * Computes the plane's coefficients from the given coordinate values. |
| * This method uses a linear regression in the least-square sense, with the assumption that |
| * the (<var>x</var>,<var>y</var>) values are precise and all uncertainty is in <var>z</var>. |
| * {@link Double#NaN} values are ignored. The result is undetermined if all points are colinear. |
| * |
| * <p>The default implementation delegates to {@link #fit(Iterable)}.</p> |
| * |
| * @param x vector of <var>x</var> coordinates. |
| * @param y vector of <var>y</var> coordinates. |
| * @param z vector of <var>z</var> values. |
| * @return an estimation of the Pearson correlation coefficient. |
| * The closer this coefficient is to +1 or -1, the better the fit. |
| * @throws IllegalArgumentException if <var>x</var>, <var>y</var> and <var>z</var> do not have the same length. |
| * |
| * @since 0.8 |
| */ |
| public double fit(final Vector x, final Vector y, final Vector z) { |
| ArgumentChecks.ensureNonNull("x", x); |
| ArgumentChecks.ensureNonNull("y", y); |
| ArgumentChecks.ensureNonNull("z", z); |
| return fit(new CompoundDirectPositions(x, y, z)); |
| } |
| |
| /** |
| * Computes the plane's coefficients from values distributed on a regular grid. Invoking this method |
| * is equivalent (except for NaN handling) to invoking {@link #fit(Vector, Vector, Vector)} where all |
| * vectors have a length of {@code nx} × {@code ny} and the <var>x</var> and <var>y</var> vectors have |
| * the following content: |
| * |
| * <blockquote> |
| * <table class="compact"> |
| * <caption><var>x</var> and <var>y</var> vectors content</caption> |
| * <tr> |
| * <th><var>x</var> vector</th> |
| * <th><var>y</var> vector</th> |
| * </tr><tr> |
| * <td> |
| * 0 1 2 3 4 5 … n<sub>x</sub>-1<br> |
| * 0 1 2 3 4 5 … n<sub>x</sub>-1<br> |
| * 0 1 2 3 4 5 … n<sub>x</sub>-1<br> |
| * …<br> |
| * 0 1 2 3 4 5 … n<sub>x</sub>-1<br> |
| * </td><td> |
| * 0 0 0 0 0 0 … 0<br> |
| * 1 1 1 1 1 1 … 1<br> |
| * 2 2 2 2 2 2 … 2<br> |
| * …<br> |
| * n<sub>y</sub>-1 n<sub>y</sub>-1 n<sub>y</sub>-1 … n<sub>y</sub>-1<br> |
| * </td> |
| * </tr> |
| * </table> |
| * </blockquote> |
| * |
| * This method uses a linear regression in the least-square sense, with the assumption that |
| * the (<var>x</var>,<var>y</var>) values are precise and all uncertainty is in <var>z</var>. |
| * The result is undetermined if all points are colinear. |
| * |
| * @param nx number of columns. |
| * @param ny number of rows. |
| * @param z values of a matrix of {@code nx} columns by {@code ny} rows organized in a row-major fashion. |
| * @return an estimation of the Pearson correlation coefficient. |
| * The closer this coefficient is to +1 or -1, the better the fit. |
| * @throws IllegalArgumentException if <var>z</var> does not have the expected length or if a <var>z</var> |
| * value is {@link Double#NaN}. |
| * |
| * @since 0.8 |
| */ |
| public double fit(final int nx, final int ny, final Vector z) { |
| ArgumentChecks.ensureStrictlyPositive("nx", nx); |
| ArgumentChecks.ensureStrictlyPositive("ny", ny); |
| ArgumentChecks.ensureNonNull("z", z); |
| final int length = Math.multiplyExact(nx, ny); |
| if (z.size() != length) { |
| throw new IllegalArgumentException(Errors.format(Errors.Keys.UnexpectedArrayLength_2, length, z.size())); |
| } |
| final Fit r = new Fit(nx, ny, z); |
| r.resolve(); |
| final double p = r.correlation(nx, length, z, null); |
| setEquation(r.sx, r.sy, r.z0); |
| return p; |
| } |
| |
| /** |
| * Computes the plane's coefficients from the given sequence of points. |
| * This method uses a linear regression in the least-square sense, with the assumption that |
| * the (<var>x</var>,<var>y</var>) values are precise and all uncertainty is in <var>z</var>. |
| * Points shall be three dimensional with coordinate values in the (<var>x</var>,<var>y</var>,<var>z</var>) order. |
| * {@link Double#NaN} values are ignored. The result is undetermined if all points are colinear. |
| * |
| * @param points the three-dimensional points. |
| * @return an estimation of the Pearson correlation coefficient. |
| * The closer this coefficient is to +1 or -1, the better the fit. |
| * @throws MismatchedDimensionException if a point is not three-dimensional. |
| */ |
| public double fit(final Iterable<? extends DirectPosition> points) { |
| final Fit r = new Fit(Objects.requireNonNull(points)); |
| r.resolve(); |
| final double p = r.correlation(0, 0, null, points.iterator()); |
| /* |
| * Store the result only when we are done, so we have a "all or nothing" behavior. |
| * We invoke the setEquation(sx, sy, z₀) method in case the user override it. |
| */ |
| setEquation(r.sx, r.sy, r.z0); |
| return p; |
| } |
| |
| /** |
| * Computes the plane coefficients. This class needs to iterate over the points two times: |
| * one for computing the coefficients, and one for the computing the Pearson coefficient. |
| * The second pass can also opportunistically checks if some small coefficients can be replaced by zero. |
| */ |
| private static final class Fit { |
| private DoubleDouble sum_x = DoubleDouble.ZERO; |
| private DoubleDouble sum_y = DoubleDouble.ZERO; |
| private DoubleDouble sum_z = DoubleDouble.ZERO; |
| private DoubleDouble sum_xx = DoubleDouble.ZERO; |
| private DoubleDouble sum_yy = DoubleDouble.ZERO; |
| private DoubleDouble sum_xy = DoubleDouble.ZERO; |
| private DoubleDouble sum_zx = DoubleDouble.ZERO; |
| private DoubleDouble sum_zy = DoubleDouble.ZERO; |
| private final int n; |
| |
| /** Solution of the plane equation. */ DoubleDouble sx, sy, z0; |
| |
| /** |
| * Computes the values of all {@code sum_*} fields from randomly distributed points. |
| * Value of all other fields are undetermined.. |
| */ |
| Fit(final Iterable<? extends DirectPosition> points) { |
| int i = 0, n = 0; |
| for (final DirectPosition p : points) { |
| final int dimension = p.getDimension(); |
| if (dimension != DIMENSION) { |
| throw new MismatchedDimensionException(Errors.format(Errors.Keys.MismatchedDimension_3, |
| Strings.toIndexed("points", i), DIMENSION, dimension)); |
| } |
| i++; |
| final double x = p.getCoordinate(0); if (Double.isNaN(x)) continue; |
| final double y = p.getCoordinate(1); if (Double.isNaN(y)) continue; |
| final double z = p.getCoordinate(2); if (Double.isNaN(z)) continue; |
| sum_x = sum_x .add(x, false); |
| sum_y = sum_y .add(y, false); |
| sum_z = sum_z .add(z, false); |
| sum_xx = sum_xx.add(DoubleDouble.product(x, x)); |
| sum_yy = sum_yy.add(DoubleDouble.product(y, y)); |
| sum_xy = sum_xy.add(DoubleDouble.product(x, y)); |
| sum_zx = sum_zx.add(DoubleDouble.product(z, x)); |
| sum_zy = sum_zy.add(DoubleDouble.product(z, y)); |
| n++; |
| } |
| this.n = n; |
| } |
| |
| /** |
| * Computes the values of all {@code sum_*} fields from the <var>z</var> values on a regular grid. |
| * Value of all other fields are undetermined.. |
| */ |
| Fit(final int nx, final int ny, final Vector vz) { |
| /* |
| * Computes the sum of x, y and z values. Computes also the sum of x², y², x⋅y, z⋅x and z⋅y values. |
| * When possible, we will avoid to compute the sum inside the loop and use the following identities |
| * instead: |
| * |
| * 1 + 2 + 3 ... + n = n⋅(n+1)/2 (arithmetic series) |
| * 1² + 2² + 3² ... + n² = n⋅(n+0.5)⋅(n+1)/3 |
| * |
| * Note that for exclusive upper bound, we need to replace n by n-1 in above formulas. |
| */ |
| int n = 0; |
| for (int y=0; y<ny; y++) { |
| for (int x=0; x<nx; x++) { |
| final double z = vz.doubleValue(n); |
| if (Double.isNaN(z)) { |
| throw new IllegalArgumentException(Errors.format(Errors.Keys.NotANumber_1, Strings.toIndexed("z", n))); |
| } |
| sum_z = sum_z .add(z, false); |
| sum_zx = sum_zx.add(DoubleDouble.product(z, x)); |
| sum_zy = sum_zy.add(DoubleDouble.product(z, y)); |
| n++; |
| } |
| } |
| sum_x = DoubleDouble.of(n/2d, false).multiply(nx-1); // Division by 2 is exact. |
| sum_y = DoubleDouble.of(n/2d, false).multiply(ny-1); |
| sum_xx = DoubleDouble.of(n) .multiply(nx-0.5, false).multiply(nx-1).divide(3); |
| sum_yy = DoubleDouble.of(n) .multiply(ny-0.5, false).multiply(ny-1).divide(3); |
| sum_xy = DoubleDouble.of(n/4d, false).multiply(ny-1) .multiply(nx-1); |
| this.n = n; |
| } |
| |
| /** |
| * Computes the {@link #sx}, {@link #sy} and {@link #z0} values using the sums computed by the constructor. |
| */ |
| private void resolve() { |
| /* |
| * ( sum_zx - sum_z⋅sum_x ) = sx⋅(sum_xx - sum_x⋅sum_x) + sy⋅(sum_xy - sum_x⋅sum_y) |
| * ( sum_zy - sum_z⋅sum_y ) = sx⋅(sum_xy - sum_x⋅sum_y) + sy⋅(sum_yy - sum_y⋅sum_y) |
| */ |
| var zx = sum_zx.subtract(sum_z.multiply(sum_x).divide(n)); // zx = sum_zx - sum_z⋅sum_x/n |
| var zy = sum_zy.subtract(sum_z.multiply(sum_y).divide(n)); // zy = sum_zy - sum_z⋅sum_y/n |
| var xx = sum_xx.subtract(sum_x.multiply(sum_x).divide(n)); // xx = sum_xx - sum_x⋅sum_x/n |
| var xy = sum_xy.subtract(sum_x.multiply(sum_y).divide(n)); // xy = sum_xy - sum_x⋅sum_y/n |
| var yy = sum_yy.subtract(sum_y.multiply(sum_y).divide(n)); // yy = sum_yy - sum_y⋅sum_y/n |
| /* |
| * den = (xy⋅xy - xx⋅yy) |
| */ |
| final DoubleDouble den = xy.square().subtract(xx.multiply(yy)); |
| /* |
| * sx = (zy⋅xy - zx⋅yy) / den |
| * sy = (zx⋅xy - zy⋅xx) / den |
| * z₀ = (sum_z - (sx⋅sum_x + sy⋅sum_y)) / n |
| */ |
| sx = zy.multiply(xy).subtract(zx.multiply(yy)).divide(den); |
| sy = zx.multiply(xy).subtract(zy.multiply(xx)).divide(den); |
| z0 = sum_z.subtract(sx.multiply(sum_x).add(sy.multiply(sum_y))).divide(n); |
| } |
| |
| /** |
| * Computes an estimation of the Pearson correlation coefficient. We do not use double-double arithmetic |
| * here since the Pearson coefficient is for information purpose (quality estimation). |
| * |
| * <p>Only one of ({@code nx}, {@code length}, {@code z}) tuple or {@code points} argument should be non-null.</p> |
| */ |
| double correlation(final int nx, final int length, final Vector vz, |
| final Iterator<? extends DirectPosition> points) |
| { |
| boolean detectZeroSx = true; |
| boolean detectZeroSy = true; |
| boolean detectZeroZ0 = true; |
| final double sx = this.sx.doubleValue(); |
| final double sy = this.sy.doubleValue(); |
| final double z0 = this.z0.doubleValue(); |
| final double mean_x = sum_x.doubleValue() / n; |
| final double mean_y = sum_y.doubleValue() / n; |
| final double mean_z = sum_z.doubleValue() / n; |
| final double offset = abs((sx * mean_x + sy * mean_y) + z0); // Offsetted z₀ - see comment before usage. |
| int index = 0; |
| double sum_ds2 = 0, sum_dz2 = 0, sum_dsz = 0; |
| for (;;) { |
| double x, y, z; |
| if (vz != null) { |
| if (index >= length) break; |
| x = index % nx; |
| y = index / nx; |
| z = vz.doubleValue(index++); |
| } else { |
| if (!points.hasNext()) break; |
| final DirectPosition p = points.next(); |
| x = p.getCoordinate(0); |
| y = p.getCoordinate(1); |
| z = p.getCoordinate(2); |
| } |
| x = (x - mean_x) * sx; |
| y = (y - mean_y) * sy; |
| z = (z - mean_z); |
| final double s = x + y; |
| if (!Double.isNaN(s) && !Double.isNaN(z)) { |
| sum_ds2 += s * s; |
| sum_dz2 += z * z; |
| sum_dsz += s * z; |
| } |
| /* |
| * Algorithm for detecting if a coefficient should be zero: |
| * If for every points given by the user, adding (sx⋅x) in (sx⋅x + sy⋅y + z₀) does not make any difference |
| * because (sx⋅x) is smaller than 1 ULP of (sy⋅y + z₀), then it is not worth adding it and sx can be set |
| * to zero. The same rational applies to (sy⋅y) and z₀. |
| * |
| * Since we work with differences from the means, the z = sx⋅x + sy⋅y + z₀ equation can be rewritten as: |
| * |
| * Δz = sx⋅Δx + sy⋅Δy + (sx⋅mx + sy⋅my + z₀ - mz) where the term between (…) is close to zero. |
| * |
| * The check for (sx⋅Δx) and (sy⋅Δy) below ignore the (…) term since it is close to zero. |
| * The check for z₀ is derived from an equation without the -mz term. |
| */ |
| if (detectZeroSx && abs(x) >= ulp(y * ZERO_THRESHOLD)) detectZeroSx = false; |
| if (detectZeroSy && abs(y) >= ulp(x * ZERO_THRESHOLD)) detectZeroSy = false; |
| if (detectZeroZ0 && offset >= ulp(s * ZERO_THRESHOLD)) detectZeroZ0 = false; |
| } |
| if (detectZeroSx) this.sx = DoubleDouble.ZERO; |
| if (detectZeroSy) this.sy = DoubleDouble.ZERO; |
| if (detectZeroZ0) this.z0 = DoubleDouble.ZERO; |
| return Math.min(sum_dsz / sqrt(sum_ds2 * sum_dz2), 1); |
| } |
| } |
| |
| /** |
| * Returns a clone of this plane. |
| * |
| * @return a clone of this plane. |
| */ |
| @Override |
| public Plane clone() { |
| try { |
| return (Plane) super.clone(); |
| } catch (CloneNotSupportedException exception) { |
| throw new AssertionError(exception); |
| } |
| } |
| |
| /** |
| * Compares this plane with the specified object for equality. |
| * |
| * @param object the object to compare with this plane for equality. |
| * @return {@code true} if both objects are equal. |
| */ |
| @Override |
| public boolean equals(final Object object) { |
| if (object != null && getClass() == object.getClass()) { |
| final Plane that = (Plane) object; |
| return Numerics.equals(this.z0, that.z0 ) && |
| Numerics.equals(this.sx, that.sx) && |
| Numerics.equals(this.sy, that.sy); |
| } else { |
| return false; |
| } |
| } |
| |
| /** |
| * Returns a hash code value for this plane. |
| */ |
| @Override |
| public int hashCode() { |
| return Long.hashCode(serialVersionUID |
| ^ (Double.doubleToLongBits(z0) |
| + 31 * (Double.doubleToLongBits(sx) |
| + 31 * (Double.doubleToLongBits(sy))))); |
| } |
| |
| /** |
| * Returns a string representation of this plane. |
| * The string will contain the plane's equation, as below: |
| * |
| * <blockquote> |
| * <var>z</var>(<var>x</var>,<var>y</var>) = {@linkplain #slopeX() sx}⋅<var>x</var> |
| * + {@linkplain #slopeY() sy}⋅<var>y</var> + {@linkplain #z0() z₀} |
| * </blockquote> |
| */ |
| @Override |
| public String toString() { |
| final StringBuilder buffer = new StringBuilder(60).append("z(x,y) = "); |
| String separator = ""; |
| if (sx != 0) { |
| buffer.append(sx).append("⋅x"); |
| separator = " + "; |
| } |
| if (sy != 0) { |
| buffer.append(separator).append(sy).append("⋅y"); |
| separator = " + "; |
| } |
| return buffer.append(separator).append(z0).toString(); |
| } |
| } |
