/*
 * 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.referencing.datum;

import java.util.Date;
import java.util.Arrays;
import java.io.Serializable;
import org.opengis.metadata.extent.Extent;
import org.opengis.referencing.datum.GeodeticDatum;
import org.opengis.referencing.datum.PrimeMeridian;
import org.opengis.referencing.operation.Matrix;
import org.apache.sis.referencing.operation.matrix.Matrix4;
import org.apache.sis.referencing.operation.matrix.Matrices;
import org.apache.sis.referencing.operation.matrix.MatrixSIS;
import org.apache.sis.referencing.IdentifiedObjects;
import org.apache.sis.io.wkt.FormattableObject;
import org.apache.sis.io.wkt.Formatter;
import org.apache.sis.util.Utilities;
import org.apache.sis.util.resources.Errors;
import org.apache.sis.internal.util.DoubleDouble;
import org.apache.sis.internal.referencing.WKTKeywords;
import org.apache.sis.internal.referencing.Resources;

import static java.lang.Math.abs;
import static org.apache.sis.util.ArgumentChecks.*;
import static org.apache.sis.referencing.operation.matrix.Matrix4.SIZE;

// Branch-dependent imports
import java.util.Objects;


/**
 * Parameters for a geographic transformation between two datum having the same prime meridian.
 * Bursa-Wolf parameters are also known as <cite>Helmert transformation parameters</cite>.
 * For an explanation of their purpose, see the <cite>Bursa-Wolf parameters</cite> section
 * of {@link DefaultGeodeticDatum} class javadoc.
 *
 * <p>The Bursa-Wolf parameters shall be applied to geocentric coordinates,
 * where the <var>X</var> axis points towards the Prime Meridian (usually Greenwich),
 * the <var>Y</var> axis points East, and the <var>Z</var> axis points North.</p>
 *
 * <div class="note"><b>Note:</b>
 * The upper case letters are intentional. By convention, (<var>X</var>, <var>Y</var>, <var>Z</var>)
 * stand for <cite>geocentric</cite> coordinates while (<var>x</var>, <var>y</var>, <var>z</var>)
 * stand for <cite>projected</cite> coordinates.</div>
 *
 * The "Bursa-Wolf" formula is expressed with 7 parameters, listed in the table below.
 * The <cite>code</cite>, <cite>name</cite> and <cite>abbreviation</cite> columns list EPSG identifiers,
 * while the <cite>legacy</cite> column lists the identifiers used in the legacy OGC 01-009 specification
 * (still used in some <cite>Well Known Texts</cite>).
 *
 * <div class="horizontal-flow">
 * <div><table class="sis">
 *   <caption>Parameters defined by EPSG</caption>
 *   <tr><th>Code</th> <th>Name</th>               <th>Abbr.</th>       <th>Legacy</th></tr>
 *   <tr><td>8605</td> <td>X-axis translation</td> <td>{@link #tX}</td> <td>{@code dx}</td></tr>
 *   <tr><td>8606</td> <td>Y-axis translation</td> <td>{@link #tY}</td> <td>{@code dy}</td></tr>
 *   <tr><td>8607</td> <td>Z-axis translation</td> <td>{@link #tZ}</td> <td>{@code dz}</td></tr>
 *   <tr><td>8608</td> <td>X-axis rotation</td>    <td>{@link #rX}</td> <td>{@code ex}</td></tr>
 *   <tr><td>8609</td> <td>Y-axis rotation</td>    <td>{@link #rY}</td> <td>{@code ey}</td></tr>
 *   <tr><td>8610</td> <td>Z-axis rotation</td>    <td>{@link #rZ}</td> <td>{@code ez}</td></tr>
 *   <tr><td>8611</td> <td>Scale difference</td>   <td>{@link #dS}</td> <td>{@code ppm}</td></tr>
 * </table></div>
 * <div><p><b>Geocentric coordinates transformation</b></p>
 * <p>from (<var>X</var><sub>s</sub>, <var>Y</var><sub>s</sub>, <var>Z</var><sub>s</sub>)
 *      to (<var>X</var><sub>t</sub>, <var>Y</var><sub>t</sub>, <var>Z</var><sub>t</sub>)
 * <br><span style="font-size:small">(ignoring unit conversions)</span></p>
 *
 * <p>{@include formulas.html#Bursa-Wolf}</p>
 * </div></div>
 *
 * The numerical fields in this {@code BursaWolfParameters} class use the EPSG abbreviations
 * with 4 additional constraints compared to the EPSG definitions:
 *
 * <ul>
 *   <li>Unit of scale difference ({@link #dS}) is fixed to <em>parts per million</em>.</li>
 *   <li>Unit of translation terms ({@link #tX}, {@link #tY}, {@link #tZ}) is fixed to <em>metres</em>.</li>
 *   <li>Unit of rotation terms ({@link #rX}, {@link #rY}, {@link #rZ}) is fixed to <em>arc-seconds</em>.</li>
 *   <li>Sign of rotation terms is fixed to the <em>Position Vector</em> convention (EPSG operation method 9606).
 *       This is the opposite sign than the <cite>Coordinate Frame Rotation</cite> (EPSG operation method 9607).
 *       The Position Vector convention is used by IAG and recommended by ISO 19111.</li>
 * </ul>
 *
 * <h2>Source and target geodetic datum</h2>
 * The <var>source datum</var> in above coordinates transformation is the {@link DefaultGeodeticDatum} instance
 * that contain this {@code BursaWolfParameters}. It can be any datum, including datum that are valid only locally.
 * The <var>{@linkplain #getTargetDatum() target datum}</var> is specified at construction time and is often,
 * but not necessarily, the <cite>World Geodetic System 1984</cite> (WGS 84) datum.
 *
 * <p>If the source and target datum does not have the same {@linkplain DefaultGeodeticDatum#getPrimeMeridian()
 * prime meridian}, then it is user's responsibility to apply longitude rotation before to use the Bursa-Wolf
 * parameters.</p>
 *
 * <h2>When Bursa-Wolf parameters are used</h2>
 * {@code BursaWolfParameters} are used in three contexts:
 * <ol>
 *   <li>Created as a step while creating a {@linkplain org.apache.sis.referencing.operation.AbstractCoordinateOperation
 *       coordinate operation} from the EPSG database.</li>
 *   <li>Associated to a {@link DefaultGeodeticDatum} with the WGS 84 {@linkplain #getTargetDatum() target datum} for
 *       providing the parameter values to display in the {@code TOWGS84[…]} element of <cite>Well Known Text</cite>
 *       (WKT) version 1. Note that WKT version 2 does not have {@code TOWGS84[…]} element anymore.</li>
 *   <li>Specified at {@code DefaultGeodeticDatum} construction time for arbitrary target datum.
 *       Apache SIS will ignore those Bursa-Wolf parameters, except as a fallback if no parameters
 *       can been found in the EPSG database for a given pair of source and target CRS.</li>
 * </ol>
 *
 * <div class="note"><b>Note:</b>
 * In EPSG terminology, Apache SIS gives precedence to the <cite>late-binding</cite> approach
 * (case 1 above) over the <cite>early-binding</cite> approach (case 3 above).</div>
 *
 * @author  Martin Desruisseaux (IRD, Geomatys)
 * @version 0.7
 *
 * @see DefaultGeodeticDatum#getBursaWolfParameters()
 * @see <a href="http://en.wikipedia.org/wiki/Helmert_transformation">Wikipedia: Helmert transformation</a>
 *
 * @since 0.4
 * @module
 */
public class BursaWolfParameters extends FormattableObject implements Cloneable, Serializable {
    /**
     * Serial number for inter-operability with different versions.
     */
    private static final long serialVersionUID = 754825592343010900L;

    /**
     * The conversion factor from <cite>parts per million</cite> to scale minus one.
     */
    static final double PPM = 1E+6;

    /**
     * X-axis translation in metres (EPSG:8605).
     * The legacy OGC parameter name is {@code "dx"}.
     */
    public double tX;

    /**
     * Y-axis translation in metres (EPSG:8606).
     * The legacy OGC parameter name is {@code "dy"}.
     */
    public double tY;

    /**
     * Z-axis translation in metres (EPSG:8607).
     * The legacy OGC parameter name is {@code "dz"}.
     */
    public double tZ;

    /**
     * X-axis rotation in arc-seconds (EPSG:8608), sign following the <cite>Position Vector</cite> convention.
     * The legacy OGC parameter name is {@code "ex"}.
     */
    public double rX;

    /**
     * Y-axis rotation in arc-seconds (EPSG:8609), sign following the <cite>Position Vector</cite> convention.
     * The legacy OGC parameter name is {@code "ey"}.
     */
    public double rY;

    /**
     * Z-axis rotation in arc-seconds (EPSG:8610), sign following the <cite>Position Vector</cite> convention.
     * The legacy OGC parameter name is {@code "ez"}.
     */
    public double rZ;

    /**
     * The scale difference in parts per million (EPSG:8611).
     * The legacy OGC parameter name is {@code "ppm"}.
     *
     * <div class="note"><b>Example:</b>
     * If a distance of 100 km in the source coordinate reference system translates into a distance of 100.001 km
     * in the target coordinate reference system, the scale difference is 1 ppm (the ratio being 1.000001).</div>
     */
    public double dS;

    /**
     * The target datum for this set of parameters, or {@code null} if unknown.
     * This is usually the WGS 84 datum, but other targets are allowed.
     *
     * <p>The source datum is the {@link DefaultGeodeticDatum} that contain this {@code BursaWolfParameters}
     * instance.</p>
     *
     * @see #getTargetDatum()
     */
    private final GeodeticDatum targetDatum;

    /**
     * Region or timeframe in which a coordinate transformation based on those Bursa-Wolf parameters is valid,
     * or {@code null} if unspecified.
     *
     * @see #getDomainOfValidity()
     */
    private final Extent domainOfValidity;

    /**
     * Creates a new instance for the given target datum and domain of validity.
     * All numerical parameters are initialized to 0, which correspond to an identity transform.
     * Callers can assign numerical values to the public fields of interest after construction.
     * For example, many coordinate transformations will provide values only for the translation
     * terms ({@link #tX}, {@link #tY}, {@link #tZ}).
     *
     * <p>Alternatively, numerical fields can also be initialized by a call to
     * {@link #setPositionVectorTransformation(Matrix, double)}.</p>
     *
     * @param targetDatum       the target datum (usually WGS 84) for this set of parameters, or {@code null} if unknown.
     * @param domainOfValidity  area or region in which a coordinate transformation based on those Bursa-Wolf parameters
     *                          is valid, or {@code null} is unspecified.
     */
    public BursaWolfParameters(final GeodeticDatum targetDatum, final Extent domainOfValidity) {
        this.targetDatum = targetDatum;
        this.domainOfValidity = domainOfValidity;
    }

    /**
     * Verifies parameters validity after initialization of {@link DefaultGeodeticDatum}.
     * This method requires that the prime meridian of the target datum is either the same
     * than the enclosing {@code GeodeticDatum}, or Greenwich. We put this restriction for
     * avoiding ambiguity about whether the longitude rotation should be applied before or
     * after the datum shift.
     *
     * <p>If the target prime meridian is Greenwich, then SIS will assume that the datum shift
     * needs to be applied in a coordinate system having Greenwich as the prime meridian.</p>
     *
     * <p><b>Maintenance note:</b>
     * if the above policy regarding prime meridians is modified, then some {@code createOperationStep(…)} method
     * implementations in {@link org.apache.sis.referencing.operation.CoordinateOperationFinder} may need to be
     * revisited. See especially the methods creating a transformation between a pair of {@code GeocentricCRS} or
     * between a pair of {@code GeographicCRS} (tip: search for {@code DefaultGeodeticDatum}).</p>
     *
     * @param  pm  the prime meridian of the enclosing {@code GeodeticDatum}.
     */
    void verify(final PrimeMeridian pm) throws IllegalArgumentException {
        if (targetDatum != null) {
            final PrimeMeridian actual = targetDatum.getPrimeMeridian();
            if (actual.getGreenwichLongitude() != 0 && !Utilities.equalsIgnoreMetadata(pm, actual)) {
                throw new IllegalArgumentException(Resources.format(Resources.Keys.MismatchedPrimeMeridian_2,
                        IdentifiedObjects.getName(pm, null), IdentifiedObjects.getName(actual, null)));
            }
        }
        ensureFinite("tX", tX);
        ensureFinite("tY", tY);
        ensureFinite("tZ", tZ);
        ensureFinite("rX", rX);
        ensureFinite("rY", rY);
        ensureFinite("rZ", rZ);
        ensureBetween("dS", -PPM, PPM, dS);     // For preventing zero or negative value on the matrix diagonal.
    }

    /**
     * Returns the target datum for this set of parameters, or {@code null} if unknown.
     * This is usually the WGS 84 datum, but other targets are allowed.
     *
     * <p>The source datum is the {@link DefaultGeodeticDatum} that contain this {@code BursaWolfParameters}
     * instance.</p>
     *
     * @return the target datum for this set of parameters, or {@code null} if unknown.
     */
    public GeodeticDatum getTargetDatum() {
        return targetDatum;
    }

    /**
     * Returns the parameter values. The length of the returned array depends on the values:
     *
     * <ul>
     *   <li>If this instance is an {@link TimeDependentBWP}, then the array length will be 14.</li>
     *   <li>Otherwise if this instance contains a non-zero {@link #dS} value, then the array length will be 7 with
     *       {@link #tX}, {@link #tY}, {@link #tZ}, {@link #rX}, {@link #rY}, {@link #rZ} and {@link #dS} values
     *       in that order.</li>
     *   <li>Otherwise if this instance contains non-zero rotation terms,
     *       then this method returns the first 6 of the above-cited values.</li>
     *   <li>Otherwise (i.e. this instance {@linkplain #isTranslation() is a translation}),
     *       this method returns only the first 3 of the above-cited values.</li>
     * </ul>
     *
     * <div class="note"><b>Note:</b>
     * the rules about the arrays of length 3, 6 or 7 are derived from the <cite>Well Known Text</cite> (WKT)
     * version 1 specification. The rule about the array of length 14 is an extension.</div>
     *
     * @return the parameter values as an array of length 3, 6, 7 or 14.
     *
     * @since 0.6
     */
    @SuppressWarnings("fallthrough")
    public double[] getValues() {
        final double[] elements = new double[(dS != 0) ? 7 : (rZ != 0 || rY != 0 || rX != 0) ? 6 : 3];
        switch (elements.length) {
            default: elements[6] = dS;      // Fallthrough everywhere.
            case 6:  elements[5] = rZ;
                     elements[4] = rY;
                     elements[3] = rX;
            case 3:  elements[2] = tZ;
                     elements[1] = tY;
                     elements[0] = tX;
        }
        return elements;
    }

    /**
     * Sets the parameters to the given values. The given array can have any length. The first array elements will be
     * assigned to the {@link #tX}, {@link #tY}, {@link #tZ}, {@link #rX}, {@link #rY}, {@link #rZ} and {@link #dS}
     * fields in that order.
     *
     * <ul>
     *   <li>If the length of the given array is not sufficient for assigning a value to every fields,
     *       then the remaining fields are left unchanged (they are <strong>not</strong> reset to zero,
     *       but this is not a problem if this {@code BursaWolfParameters} is a new instance).</li>
     *   <li>If the length of the given array is greater than necessary, then extra elements are ignored by this base
     *       class. Note however that those extra elements may be used by subclasses like {@link TimeDependentBWP}.</li>
     * </ul>
     *
     * @param  elements  the new parameter values, as an array of any length.
     *
     * @since 0.6
     */
    @SuppressWarnings("fallthrough")
    public void setValues(final double... elements) {
        switch (elements.length) {
            default: dS = elements[6];      // Fallthrough everywhere.
            case 6:  rZ = elements[5];
            case 5:  rY = elements[4];
            case 4:  rX = elements[3];
            case 3:  tZ = elements[2];
            case 2:  tY = elements[1];
            case 1:  tX = elements[0];
            case 0:  break;
         }
    }

    /**
     * Returns {@code true} if the {@linkplain #targetDatum target datum} is equals (at least on computation purpose)
     * to the WGS84 datum. If the datum is unspecified, then this method returns {@code true} since WGS84 is the only
     * datum supported by the WKT 1 format, and is what users often mean.
     *
     * @return {@code true} if the given datum is equal to WGS84 for computational purpose.
     */
    final boolean isToWGS84() {
        return (targetDatum != null) && IdentifiedObjects.isHeuristicMatchForName(targetDatum, "WGS84");
    }

    /**
     * Returns {@code true} if a transformation built from this set of parameters would perform no operation.
     * This is true when the value of all parameters is zero.
     *
     * @return {@code true} if the parameters describe no operation.
     */
    public boolean isIdentity() {
        return tX == 0 && tY == 0 && tZ == 0 && isTranslation();
    }

    /**
     * Returns {@code true} if a transformation built from this set of parameters would perform only a translation.
     *
     * @return {@code true} if the parameters describe a translation only.
     */
    public boolean isTranslation() {
        return rX == 0 && rY == 0 && rZ == 0 && dS == 0;
    }

    /**
     * Inverts in-place the sign of rotation terms ({@link #rX}, {@link #rY}, {@link #rZ}).
     * This method can be invoked for converting a <cite>Coordinate Frame Rotation</cite> transformation
     * (EPSG operation method 9607) to a <em>Position Vector</em> transformation (EPSG operation method 9606).
     * The later convention is used by IAG and recommended by ISO 19111.
     */
    public void reverseRotation() {
        rX = -rX;
        rY = -rY;
        rZ = -rZ;
    }

    /**
     * Inverts in-place the transformation by inverting the sign of all numerical parameters.
     * The {@linkplain #getPositionVectorTransformation(Date) position vector transformation} matrix
     * created from inverted Bursa-Wolf parameters will be <strong>approximately</strong> equals
     * to the {@linkplain org.apache.sis.referencing.operation.matrix.MatrixSIS#inverse() inverse}
     * of the matrix created from the original parameters. The equality holds approximately only
     * because the parameter values are very small (parts per millions and arc-seconds).
     */
    public void invert() {
        final double[] values = getValues();
        for (int i=0; i<values.length; i++) {
            values[i] = -values[i];
        }
        setValues(values);
    }

    /**
     * Returns the elapsed time from the {@linkplain TimeDependentBWP#getTimeReference() reference time}
     * to the given date, in millennium. If this {@code BursaWolfParameters} is not time-dependent, then
     * returns {@code null}.
     */
    DoubleDouble period(final Date time) {
        return null;
    }

    /**
     * Returns the parameter at the given index. If this {@code BursaWolfParameters} is time-dependent,
     * then the returned value shall be corrected for the given period.
     *
     * @param  index   0 for {@code tX}, 1 for {@code tY}, <i>etc.</i> in {@code TOWGS84[…]} order.
     * @param  period  the value computed by {@link #period(Date)}, or {@code null}.
     */
    DoubleDouble param(final int index, final DoubleDouble period) {
        final double p;
        switch (index) {
            case 0: p = tX; break;
            case 1: p = tY; break;
            case 2: p = tZ; break;
            case 3: p = rX; break;
            case 4: p = rY; break;
            case 5: p = rZ; break;
            case 6: p = dS; break;
            default: throw new AssertionError(index);
        }
        return DoubleDouble.createAndGuessError(p);
    }

    /**
     * Returns the position vector transformation (geocentric domain) as an affine transform.
     * For transformations that do not depend on time, the formula is as below where {@code R}
     * is a conversion factor from arc-seconds to radians:
     *
     * <blockquote><pre> R = toRadians(1″)
     * S = 1 + {@linkplain #dS}/1000000
     * ┌    ┐    ┌                               ┐  ┌   ┐
     * │ X' │    │      S   -{@linkplain #rZ}*RS   +{@linkplain #rY}*RS   {@linkplain #tX} │  │ X │
     * │ Y' │  = │ +{@linkplain #rZ}*RS        S   -{@linkplain #rX}*RS   {@linkplain #tY} │  │ Y │
     * │ Z' │    │ -{@linkplain #rY}*RS   +{@linkplain #rX}*RS        S   {@linkplain #tZ} │  │ Z │
     * │ 1  │    │      0        0        0    1 │  │ 1 │
     * └    ┘    └                               ┘  └   ┘</pre></blockquote>
     *
     * This affine transform can be applied on <strong>geocentric</strong> coordinates.
     * This is identified as operation method 1033 in the EPSG database.
     * Those geocentric coordinates are typically converted from geographic coordinates
     * in the region or timeframe given by {@link #getDomainOfValidity()}.
     *
     * <p>If the source datum and the {@linkplain #getTargetDatum() target datum} do not use the same
     * {@linkplain DefaultGeodeticDatum#getPrimeMeridian() prime meridian}, then it is caller's responsibility
     * to apply longitude rotation before to use the matrix returned by this method.</p>
     *
     * <h4>Time-dependent transformation</h4>
     * Some transformations use parameters that vary with time (e.g. operation method EPSG:1053).
     * Users can optionally specify a date for which the transformation is desired.
     * For transformations that do not depends on time, this date is ignored and can be null.
     * For time-dependent transformations, {@code null} values default to the transformation's
     * {@linkplain TimeDependentBWP#getTimeReference() reference time}.
     *
     * <h4>Inverse transformation</h4>
     * The inverse transformation can be approximated by reversing the sign of the 7 parameters before to use them
     * in the above matrix. This is often considered sufficient since <cite>position vector transformations</cite>
     * are themselves approximations. However Apache SIS will rather use
     * {@link org.apache.sis.referencing.operation.matrix.MatrixSIS#inverse()} in order to increase the chances
     * that concatenation of transformations <var>A</var> → <var>B</var> followed by <var>B</var> → <var>A</var>
     * gives back the identity transform.
     *
     * @param  time  date for which the transformation is desired, or {@code null} for the transformation's reference time.
     * @return an affine transform in geocentric space created from this Bursa-Wolf parameters and the given time.
     *
     * @see DefaultGeodeticDatum#getPositionVectorTransformation(GeodeticDatum, Extent)
     */
    public Matrix getPositionVectorTransformation(final Date time) {
        final DoubleDouble period = period(time);
        if (period == null && isTranslation()) {
            final Matrix4 matrix = new Matrix4();
            matrix.m03 = tX;
            matrix.m13 = tY;
            matrix.m23 = tZ;
            return matrix;
        }
        /*
         * Above was an optimization for the common case where the Bursa-Wolf parameters contain only
         * translation terms. If we have rotation or scale terms, then use double-double arithmetic.
         */
        final DoubleDouble RS = DoubleDouble.createSecondsToRadians();
        final DoubleDouble S = param(6, period);
        S.divide(PPM);
        S.add(1);                                                   // S = 1 + dS / PPM;
        RS.multiply(S);                                             // RS = toRadians(1″) * S;
        final DoubleDouble  X = param(3, period); X.multiply(RS);
        final DoubleDouble  Y = param(4, period); Y.multiply(RS);
        final DoubleDouble  Z = param(5, period); Z.multiply(RS);
        final DoubleDouble mX = new DoubleDouble(X); mX.negate();
        final DoubleDouble mY = new DoubleDouble(Y); mY.negate();
        final DoubleDouble mZ = new DoubleDouble(Z); mZ.negate();
        final Integer       O = 0;                                  // Fetch Integer instance only once.
        return Matrices.create(4, 4, new Number[] {
                 S,  mZ,   Y,  param(0, period),
                 Z,   S,  mX,  param(1, period),
                mY,   X,   S,  param(2, period),
                 O,   O,   O,  1});
    }

    /**
     * Sets all Bursa-Wolf parameters from the given <cite>Position Vector transformation</cite> matrix.
     * The matrix shall comply to the following constraints:
     *
     * <ul>
     *   <li>The matrix shall be {@linkplain org.apache.sis.referencing.operation.matrix.MatrixSIS#isAffine() affine}.</li>
     *   <li>The sub-matrix defined by {@code matrix} without the last row and last column shall be
     *       <a href="http://en.wikipedia.org/wiki/Skew-symmetric_matrix">skew-symmetric</a> (a.k.a. antisymmetric).</li>
     * </ul>
     *
     * @param  matrix     the matrix from which to get Bursa-Wolf parameters.
     * @param  tolerance  the tolerance error for the skew-symmetric matrix test, in units of PPM or arc-seconds (e.g. 1E-8).
     * @throws IllegalArgumentException if the specified matrix does not met the conditions.
     *
     * @see #getPositionVectorTransformation(Date)
     */
    public void setPositionVectorTransformation(final Matrix matrix, final double tolerance) throws IllegalArgumentException {
        final int numRow = matrix.getNumRow();
        final int numCol = matrix.getNumCol();
        if (numRow != SIZE || numCol != SIZE) {
            final Integer n = SIZE;
            throw new IllegalArgumentException(Errors.format(Errors.Keys.MismatchedMatrixSize_4, n, n, numRow, numCol));
        }
        if (!Matrices.isAffine(matrix)) {
            throw new IllegalArgumentException(Resources.format(Resources.Keys.NotAnAffineTransform));
        }
        /*
         * Translation terms, taken "as-is".
         * If the matrix contains only translation terms (which is often the case), we are done.
         */
        tX = matrix.getElement(0,3);
        tY = matrix.getElement(1,3);
        tZ = matrix.getElement(2,3);
        if (Matrices.isTranslation(matrix)) {                   // Optimization for a common case.
            return;
        }
        /*
         * Scale factor: take the average of elements on the diagonal. All those
         * elements should have the same value, but we tolerate slight deviation
         * (this will be verified later).
         */
        final DoubleDouble S = DoubleDouble.castOrCopy(getNumber(matrix, 0,0));
        S.addGuessError(getNumber(matrix, 1,1));
        S.addGuessError(getNumber(matrix, 2,2));
        S.divide(3);
        /*
         * Computes: RS = S * toRadians(1″)
         *           dS = (S-1) * PPM
         */
        final DoubleDouble RS = DoubleDouble.createSecondsToRadians();
        RS.multiply(S);
        S.add(-1);
        S.multiply(PPM);
        dS = S.doubleValue();
        /*
         * Rotation terms. Each rotation terms appear twice, with one value being the negative of the other value.
         * We verify this skew symmetric aspect in the loop. We also opportunistically verify that the scale terms
         * are uniform.
         */
        for (int j=0; j < SIZE-1; j++) {
            if (!(abs((matrix.getElement(j,j) - 1)*PPM - dS) <= tolerance)) {
                throw new IllegalArgumentException(Resources.format(Resources.Keys.NonUniformScale));
            }
            for (int i = j+1; i < SIZE-1; i++) {
                S.setFrom(RS);
                S.inverseDivideGuessError(getNumber(matrix, j,i));          // Negative rotation term.
                double value = S.value;
                double error = S.error;
                S.setFrom(RS);
                S.inverseDivideGuessError(getNumber(matrix, i,j));          // Positive rotation term.
                if (!(abs(value + S.value) <= tolerance)) {                 // We expect r1 ≈ -r2
                    throw new IllegalArgumentException(Resources.format(Resources.Keys.NotASkewSymmetricMatrix));
                }
                S.subtract(value, error);
                S.multiply(0.5);
                value = S.doubleValue();                                    // Average of the two rotation terms.
                switch (j*SIZE + i) {
                    case 1: rZ =  value; break;
                    case 2: rY = -value; break;
                    case 6: rX =  value; break;
                }
            }
        }
    }

    /**
     * Retrieves the value at the specified row and column of the given matrix, wrapped in a {@code Number}.
     * The {@code Number} type depends on the matrix accuracy.
     *
     * @param  matrix  the matrix from which to get the number.
     * @param  row     the row index, from 0 inclusive to {@link Matrix#getNumRow()} exclusive.
     * @param  column  the column index, from 0 inclusive to {@link Matrix#getNumCol()} exclusive.
     * @return the current value at the given row and column.
     */
    private static Number getNumber(final Matrix matrix, final int row, final int column) {
        if (matrix instanceof MatrixSIS) {
            return ((MatrixSIS) matrix).getNumber(row, column);
        } else {
            return matrix.getElement(row, column);
        }
    }

    /**
     * Returns the region or timeframe in which a coordinate transformation based on those Bursa-Wolf parameters is
     * valid, or {@code null} if unspecified. If an extent was specified at construction time, then that extent is
     * returned. Otherwise the datum domain of validity (which may be {@code null}) is returned.
     *
     * @return area or region or timeframe in which the coordinate transformation is valid, or {@code null}.
     *
     * @see org.apache.sis.metadata.iso.extent.DefaultExtent
     */
    public Extent getDomainOfValidity() {
        if (domainOfValidity == null && targetDatum != null) {
            return targetDatum.getDomainOfValidity();
        }
        return domainOfValidity;
    }

    /**
     * Returns a copy of this object.
     *
     * @return a copy of all parameters.
     */
    @Override
    public BursaWolfParameters clone() {
        try {
            return (BursaWolfParameters) super.clone();
        }  catch (CloneNotSupportedException exception) {
            // Should not happen, since we are cloneable.
            throw new AssertionError(exception);
        }
    }

    /**
     * Compares the specified object with this object for equality.
     *
     * @param  object  the object to compare with the parameters.
     * @return {@code true} if the given object is equal to this {@code BursaWolfParameters}.
     */
    @Override
    public boolean equals(final Object object) {
        if (object != null && object.getClass() == getClass()) {
            final BursaWolfParameters that = (BursaWolfParameters) object;
            return Arrays.equals(this.getValues(),      that.getValues()) &&
                  Objects.equals(this.targetDatum,      that.targetDatum) &&
                  Objects.equals(this.domainOfValidity, that.domainOfValidity);
        }
        return false;
    }

    /**
     * Returns a hash value for this object.
     *
     * @return the hash code value. This value does not need to be the same in past or future versions of this class.
     */
    @Override
    public int hashCode() {
        return Arrays.hashCode(getValues()) ^ (int) serialVersionUID;
    }

    /**
     * Formats this object as a <cite>Well Known Text</cite> {@code ToWGS84[…]} element.
     * The WKT contains the parameters in <var>translation</var>, <var>rotation</var>, <var>scale</var> order,
     * like below:
     *
     * <blockquote><code>TOWGS84[{@linkplain #tX}, {@linkplain #tY}, {@linkplain #tZ}, {@linkplain #rX},
     * {@linkplain #rY}, {@linkplain #rZ}, {@linkplain #dS}]</code></blockquote>
     *
     * <div class="note"><b>Compatibility note:</b>
     * {@code TOWGS84} is defined in the WKT 1 specification only.</div>
     *
     * The element name is {@code "ToWGS84"} in the common case where the {@linkplain #getTargetDatum() target datum}
     * is WGS 84. For other targets, the element name will be derived from the datum name.
     *
     * @param  formatter The formatter where to format the inner content of this WKT element.
     * @return Usually {@code "ToWGS84"}.
     */
    @Override
    protected String formatTo(final Formatter formatter) {
        final double[] values = getValues();
        for (final double value : values) {
            formatter.append(value);
        }
        if (isToWGS84()) {
            if (values.length > 7) {
                formatter.setInvalidWKT(BursaWolfParameters.class, null);
            }
            return WKTKeywords.ToWGS84;
        }
        formatter.setInvalidWKT(BursaWolfParameters.class, null);
        String name = IdentifiedObjects.getSimpleNameOrIdentifier(getTargetDatum());
        if (name == null) {
            name = "Unknown";
        }
        // We may try to build something better here in future SIS versions, if there is a need for that.
        return "To" + name;
    }
}