/*
 * 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.operation.provider;

import java.util.List;
import jakarta.xml.bind.annotation.XmlTransient;
import org.opengis.util.FactoryException;
import org.opengis.parameter.ParameterValueGroup;
import org.opengis.parameter.ParameterDescriptor;
import org.opengis.parameter.ParameterDescriptorGroup;
import org.opengis.referencing.cs.CartesianCS;
import org.opengis.referencing.cs.EllipsoidalCS;
import org.opengis.referencing.cs.CoordinateSystem;
import org.opengis.referencing.operation.Matrix;
import org.opengis.referencing.operation.MathTransform;
import org.opengis.referencing.operation.MathTransformFactory;
import org.opengis.referencing.operation.Conversion;
import org.opengis.referencing.operation.Transformation;
import org.apache.sis.referencing.IdentifiedObjects;
import org.apache.sis.referencing.util.Formulas;
import org.apache.sis.referencing.util.WKTUtilities;
import org.apache.sis.referencing.util.WKTKeywords;
import org.apache.sis.io.wkt.FormattableObject;
import org.apache.sis.io.wkt.Formatter;
import org.apache.sis.measure.Units;
import org.apache.sis.parameter.Parameters;
import org.apache.sis.parameter.Parameterized;
import org.apache.sis.parameter.ParameterBuilder;
import org.apache.sis.metadata.iso.citation.Citations;
import org.apache.sis.referencing.datum.BursaWolfParameters;
import org.apache.sis.referencing.operation.transform.MathTransforms;
import org.apache.sis.referencing.operation.transform.LinearTransform;
import org.apache.sis.util.logging.Logging;


/**
 * The base class of operation methods performing a translation, rotation and/or scale in geocentric coordinates.
 * Those methods may or may not include a Geographic/Geocentric conversion before the operation in geocentric domain,
 * depending on whether or not implementations extend the {@link GeocentricAffineBetweenGeographic} subclass.
 *
 * <h2>Note on class name</h2>
 * The {@code GeocentricAffine} class name is chosen as a generalization of {@link GeocentricTranslation}.
 * "Geocentric translations" is an operation name defined by EPSG.
 *
 * @author  Martin Desruisseaux (IRD, Geomatys)
 */
@XmlTransient
public abstract class GeocentricAffine extends GeodeticOperation {
    /**
     * The transformation type (translation, frame rotation, <i>etc.</i>).
     *
     * @todo Merge with {@link DatumShiftMethod}.
     *
     * @see #type
     */
    @SuppressWarnings("doclint:missing")
    protected enum Type {TRANSLATION, SEVEN_PARAM, FRAME_ROTATION, MOLODENSKY, CONVERSION};

    /**
     * Serial number for inter-operability with different versions.
     */
    private static final long serialVersionUID = 5597594719123422532L;

    /**
     * The tolerance factor for comparing the {@link BursaWolfParameters} values.
     * We use a tolerance of 1E-6 ({@value Formulas#LINEAR_TOLERANCE} / 10000) based on the knowledge
     * that the translation terms are in metres and the rotation terms have the same order of magnitude.
     * Actually we could use a value of zero, but we add a small tolerance for rounding errors.
     */
    private static final double BURSAWOLF_TOLERANCE = Formulas.LINEAR_TOLERANCE / 10000;

    /**
     * The operation parameter descriptor for the <cite>X-axis translation</cite>
     * ({@linkplain BursaWolfParameters#tX tX}) parameter value. Valid values range
     * from negative to positive infinity. Units are {@linkplain Units#METRE metres}.
     *
     * <!-- Generated by ParameterNameTableGenerator -->
     * <table class="sis">
     *   <caption>Parameter names</caption>
     *   <tr><td> EPSG:    </td><td> X-axis translation </td></tr>
     *   <tr><td> OGC:     </td><td> dx </td></tr>
     * </table>
     */
    public static final ParameterDescriptor<Double> TX;

    /**
     * The operation parameter descriptor for the <cite>Y-axis translation</cite>
     * ({@linkplain BursaWolfParameters#tY tY}) parameter value. Valid values range
     * from negative to positive infinity. Units are {@linkplain Units#METRE metres}.
     *
     * <!-- Generated by ParameterNameTableGenerator -->
     * <table class="sis">
     *   <caption>Parameter names</caption>
     *   <tr><td> EPSG:    </td><td> Y-axis translation </td></tr>
     *   <tr><td> OGC:     </td><td> dy </td></tr>
     * </table>
     */
    public static final ParameterDescriptor<Double> TY;

    /**
     * The operation parameter descriptor for the <cite>Z-axis translation</cite>
     * ({@linkplain BursaWolfParameters#tZ tZ}) parameter value. Valid values range
     * from negative to positive infinity. Units are {@linkplain Units#METRE metres}.
     *
     * <!-- Generated by ParameterNameTableGenerator -->
     * <table class="sis">
     *   <caption>Parameter names</caption>
     *   <tr><td> EPSG:    </td><td> Z-axis translation </td></tr>
     *   <tr><td> OGC:     </td><td> dz </td></tr>
     * </table>
     */
    public static final ParameterDescriptor<Double> TZ;

    /**
     * The operation parameter descriptor for the <cite>X-axis rotation</cite>
     * ({@linkplain BursaWolfParameters#rX rX}) parameter value.
     * Units are {@linkplain Units#ARC_SECOND arc-seconds}.
     *
     * <!-- Generated by ParameterNameTableGenerator -->
     * <table class="sis">
     *   <caption>Parameter names</caption>
     *   <tr><td> EPSG:    </td><td> X-axis rotation </td></tr>
     *   <tr><td> OGC:     </td><td> ex </td></tr>
     * </table>
     * <b>Notes:</b>
     * <ul>
     *   <li>Value domain: [-648000.0 … 648000.0]″</li>
     * </ul>
     */
    static final ParameterDescriptor<Double> RX;

    /**
     * The operation parameter descriptor for the <cite>Y-axis rotation</cite>
     * ({@linkplain BursaWolfParameters#rY rY}) parameter value.
     * Units are {@linkplain Units#ARC_SECOND arc-seconds}.
     *
     * <!-- Generated by ParameterNameTableGenerator -->
     * <table class="sis">
     *   <caption>Parameter names</caption>
     *   <tr><td> EPSG:    </td><td> Y-axis rotation </td></tr>
     *   <tr><td> OGC:     </td><td> ey </td></tr>
     * </table>
     * <b>Notes:</b>
     * <ul>
     *   <li>Value domain: [-648000.0 … 648000.0]″</li>
     * </ul>
     */
    static final ParameterDescriptor<Double> RY;

    /**
     * The operation parameter descriptor for the <cite>Z-axis rotation</cite>
     * ({@linkplain BursaWolfParameters#rZ rZ}) parameter value.
     * Units are {@linkplain Units#ARC_SECOND arc-seconds}.
     *
     * <!-- Generated by ParameterNameTableGenerator -->
     * <table class="sis">
     *   <caption>Parameter names</caption>
     *   <tr><td> EPSG:    </td><td> Z-axis rotation </td></tr>
     *   <tr><td> OGC:     </td><td> ez </td></tr>
     * </table>
     * <b>Notes:</b>
     * <ul>
     *   <li>Value domain: [-648000.0 … 648000.0]″</li>
     * </ul>
     */
    static final ParameterDescriptor<Double> RZ;

    /**
     * The operation parameter descriptor for the <cite>Scale difference</cite>
     * ({@linkplain BursaWolfParameters#dS dS}) parameter value.
     * Valid values range from negative to positive infinity.
     * Units are {@linkplain Units#PPM parts per million}.
     *
     * <!-- Generated by ParameterNameTableGenerator -->
     * <table class="sis">
     *   <caption>Parameter names</caption>
     *   <tr><td> EPSG:    </td><td> Scale difference </td></tr>
     *   <tr><td> OGC:     </td><td> ppm </td></tr>
     * </table>
     */
    static final ParameterDescriptor<Double> DS;
    static {
        final ParameterBuilder builder = builder();
        TX = createShift(builder.addIdentifier("8605").addName("X-axis translation").addName(Citations.OGC, "dx"));
        TY = createShift(builder.addIdentifier("8606").addName("Y-axis translation").addName(Citations.OGC, "dy"));
        TZ = createShift(builder.addIdentifier("8607").addName("Z-axis translation").addName(Citations.OGC, "dz"));
        RX = createRotation(builder.addIdentifier("8608"), "X-axis rotation", "ex");
        RY = createRotation(builder.addIdentifier("8609"), "Y-axis rotation", "ey");
        RZ = createRotation(builder.addIdentifier("8610"), "Z-axis rotation", "ez");
        DS = builder.addIdentifier("8611").addName("Scale difference").addName(Citations.OGC, "ppm").create(0, Units.PPM);
    }

    /**
     * Convenience method for building the rotation parameters.
     */
    private static ParameterDescriptor<Double> createRotation(final ParameterBuilder builder, final String name, final String alias) {
        return builder.addName(name).addName(Citations.OGC, alias).createBounded(-180*60*60, 180*60*60, 0, Units.ARC_SECOND);
    }

    /**
     * The transformation type (translation, frame rotation, <i>etc.</i>).
     */
    private final Type type;

    /**
     * Creates a copy of this provider.
     *
     * @deprecated This is a temporary constructor before replacement by a {@code provider()} method with JDK9.
     */
    @Deprecated
    GeocentricAffine(final GeocentricAffine copy) {
        super(copy);
        type = copy.type;
    }

    /**
     * Constructs a provider with the specified parameters.
     *
     * @param type               the operation type as an enumeration value.
     * @param parameters         description of parameters expected by this operation.
     * @param indexOfDim         number of dimensions as the index in {@link #redimensioned} array.
     * @param sourceCSType       base interface of the coordinate system of source coordinates.
     * @param sourceOnEllipsoid  whether the operation needs source ellipsoid axis lengths.
     * @param targetCSType       base interface of the coordinate system of target coordinates.
     * @param targetOnEllipsoid  whether the operation needs target ellipsoid axis lengths.
     */
    GeocentricAffine(Type operationType, ParameterDescriptorGroup parameters, int indexOfDim,
                     Class<? extends CoordinateSystem> sourceCSType, boolean sourceOnEllipsoid,
                     Class<? extends CoordinateSystem> targetCSType, boolean targetOnEllipsoid)
    {
        super((operationType == Type.CONVERSION) ? Conversion.class : Transformation.class, parameters, indexOfDim,
              sourceCSType, sourceOnEllipsoid,
              targetCSType, targetOnEllipsoid);
        type = operationType;
    }

    /**
     * Constructs a provider with the specified parameters for an operation in Cartesian space.
     */
    GeocentricAffine(Type operationType, ParameterDescriptorGroup parameters) {
        this(operationType, parameters, INDEX_OF_3D,
             CartesianCS.class, false,
             CartesianCS.class, false);
    }

    /**
     * Creates a math transform from the specified group of parameter values.
     * The default implementation creates an affine transform, but some subclasses
     * will wrap that affine operation into Geographic/Geocentric conversions.
     *
     * @param  factory  the factory to use for creating concatenated transforms.
     * @param  values   the group of parameter values.
     * @return the created math transform.
     * @throws FactoryException if a transform cannot be created.
     */
    @Override
    @SuppressWarnings("fallthrough")
    public MathTransform createMathTransform(final MathTransformFactory factory, final ParameterValueGroup values)
            throws FactoryException
    {
        final BursaWolfParameters parameters = new BursaWolfParameters(null, null);
        final Parameters pv = Parameters.castOrWrap(values);
        boolean reverseRotation = false;
        switch (type) {
            default:             throw new AssertionError(type);
            case FRAME_ROTATION: reverseRotation = true;                    // Fall through
            case SEVEN_PARAM:    parameters.rX = pv.doubleValue(RX);
                                 parameters.rY = pv.doubleValue(RY);
                                 parameters.rZ = pv.doubleValue(RZ);
                                 parameters.dS = pv.doubleValue(DS);
            case TRANSLATION:    parameters.tX = pv.doubleValue(TX);        // Fall through
                                 parameters.tY = pv.doubleValue(TY);
                                 parameters.tZ = pv.doubleValue(TZ);
        }
        if (reverseRotation) {
            parameters.reverseRotation();
        }
        return MathTransforms.linear(parameters.getPositionVectorTransformation(null));
    }

    /**
     * Creates parameter values for a Molodensky, Geocentric Translation or Position Vector transformation.
     *
     * @param  descriptor     the {@code PARAMETERS} constant of the subclass describing the operation to create.
     * @param  parameters     Bursa-Wolf parameters from which to get the values.
     * @param  isTranslation  {@code true} if the operation contains only translation terms.
     * @return the operation parameters with their values initialized.
     */
    private static Parameters createParameters(final ParameterDescriptorGroup descriptor,
            final BursaWolfParameters parameters, final boolean isTranslation)
    {
        final Parameters values = Parameters.castOrWrap(descriptor.createValue());
        values.getOrCreate(TX).setValue(parameters.tX);
        values.getOrCreate(TY).setValue(parameters.tY);
        values.getOrCreate(TZ).setValue(parameters.tZ);
        if (!isTranslation) {
            values.getOrCreate(RX).setValue(parameters.rX);
            values.getOrCreate(RY).setValue(parameters.rY);
            values.getOrCreate(RZ).setValue(parameters.rZ);
            values.getOrCreate(DS).setValue(parameters.dS);
        }
        return values;
    }

    /**
     * Returns the parameters for creating a datum shift operation.
     * The operation method will be one of the {@code GeocentricAffine} subclasses,
     * unless the specified {@code method} argument is {@link DatumShiftMethod#NONE}.
     * If no single operation method can be used, then this method returns {@code null}.
     *
     * <p>This method does <strong>not</strong> change the coordinate system type.
     * The source and target coordinate systems can be both {@code EllipsoidalCS} or both {@code CartesianCS}.
     * Any other type or mix of types (e.g. a {@code EllipsoidalCS} source and {@code CartesianCS} target)
     * will cause this method to return {@code null}. In such case, it is caller's responsibility to apply
     * the datum shift itself in Cartesian geocentric coordinates.</p>
     *
     * @param  sourceCS    the source coordinate system. Only the type and number of dimensions is checked.
     * @param  targetCS    the target coordinate system. Only the type and number of dimensions is checked.
     * @param  datumShift  the datum shift as a matrix, or {@code null} if there is no datum shift information.
     * @param  method      the preferred datum shift method. Note that {@code createParameters(…)} may overwrite.
     * @return the parameter values, or {@code null} if no single operation method can be found.
     */
    public static ParameterValueGroup createParameters(final CoordinateSystem sourceCS,
            final CoordinateSystem targetCS, final Matrix datumShift, DatumShiftMethod method)
    {
        final boolean isEllipsoidal = (sourceCS instanceof EllipsoidalCS);
        if (!(isEllipsoidal ? (targetCS instanceof EllipsoidalCS)
                            : (targetCS instanceof CartesianCS && sourceCS instanceof CartesianCS)))
        {
            return null;        // Coordinate systems are not two EllipsoidalCS or two CartesianCS.
        }
        int dimension  = sourceCS.getDimension();
        if (dimension != targetCS.getDimension()) {
            dimension  = 4;     // Any value greater than 3 means "mismatched dimensions" for this method.
        }
        if (method == DatumShiftMethod.NONE) {
            if (dimension <= 3) {
                return Affine.identity(dimension);
            } else if (isEllipsoidal) {
                final ParameterDescriptorGroup descriptor;
                switch (sourceCS.getDimension()) {
                    case 2: descriptor = Geographic2Dto3D.PARAMETERS; break;
                    case 3: descriptor = Geographic3Dto2D.PARAMETERS; break;
                    default: return null;
                }
                return descriptor.createValue();
            } else {
                return null;
            }
        }
        /*
         * Try to convert the matrix into (tX, tY, tZ, rX, rY, rZ, dS) parameters.
         * The matrix may not be convertible, in which case we will let the caller
         * uses the matrix directly in Cartesian geocentric coordinates.
         */
        final BursaWolfParameters parameters = new BursaWolfParameters(null, null);
        if (datumShift != null) try {
            parameters.setPositionVectorTransformation(datumShift, BURSAWOLF_TOLERANCE);
        } catch (IllegalArgumentException e) {
            recoverableException(GeocentricAffine.class, e);
            return null;
        } else {
            /*
             * If there is no datum shift parameters (not to be confused with identity), then those parameters
             * are assumed unknown. Using the most accurate methods would give a false impression of accuracy,
             * so we use the fastest method instead. Since all parameter values are zero, Apache SIS should use
             * the AbridgedMolodenskyTransform2D optimization.
             */
            method = DatumShiftMethod.ABRIDGED_MOLODENSKY;
        }
        final boolean isTranslation = parameters.isTranslation();
        final ParameterDescriptorGroup descriptor;
        /*
         * Following "if" blocks are ordered from most accurate to less accurate datum shift method
         * supported by GeocentricAffine subclasses (except NONE which has already been handled).
         * Special cases:
         *
         *   - If the datum shift is applied between geocentric CRS, then the Molodensky approximations do not apply
         *     as they are designed for transformations between geographic CRS only. User preference is then ignored.
         *
         *   - Molodensky methods are approximations for datum shifts having only translation terms in their Bursa-Wolf
         *     parameters. If there is also a scale or rotation terms, then we cannot use Molodensky methods. The user
         *     preference is then ignored.
         */
        if (!isEllipsoidal) {
            method = DatumShiftMethod.GEOCENTRIC_DOMAIN;
            descriptor = isTranslation ? GeocentricTranslation.PARAMETERS
                                       : PositionVector7Param .PARAMETERS;
        } else if (!isTranslation) {
            method = DatumShiftMethod.GEOCENTRIC_DOMAIN;
            descriptor = (dimension >= 3) ? PositionVector7Param3D.PARAMETERS
                                          : PositionVector7Param2D.PARAMETERS;
        } else switch (method) {
            case GEOCENTRIC_DOMAIN: {
                descriptor = (dimension >= 3) ? GeocentricTranslation3D.PARAMETERS
                                              : GeocentricTranslation2D.PARAMETERS;
                break;
            }
            case MOLODENSKY: {
                descriptor = Molodensky.PARAMETERS;
                break;
            }
            case ABRIDGED_MOLODENSKY: {
                descriptor = AbridgedMolodensky.PARAMETERS;
                break;
            }
            default: throw new AssertionError(method);
        }
        /*
         * Following lines will set all Bursa-Wolf parameter values (scale, translation
         * and rotation terms). In the particular case of Molodensky method, we have an
         * additional parameter for the number of source and target dimensions (2 or 3).
         * If the number of source and target dimensions are not the same, set to 3 and
         * let the caller adds or removes an ellipsoidal height as needed.
         */
        final Parameters values = createParameters(descriptor, parameters, isTranslation);
        switch (method) {
            case MOLODENSKY:
            case ABRIDGED_MOLODENSKY: {
                values.getOrCreate(Molodensky.DIMENSION).setValue(Math.min(dimension, 3));
                break;
            }
        }
        return values;
    }

    /**
     * Given a transformation chain, conditionally replaces the affine transform elements by an alternative object
     * showing the Bursa-Wolf parameters. The replacement is applied if and only if the affine transform is a scale,
     * translation or rotation in the geocentric domain.
     *
     * <p>This method is invoked only by {@code ConcatenatedTransform.getPseudoSteps()} for the need of WKT formatting.
     * The purpose of this method is very similar to the purpose of {@code AbstractMathTransform.beforeFormat(List, int,
     * boolean)} except that we need to perform the {@code forDatumShift(…)} work only after {@code beforeFormat(…)}
     * finished its work for all {@code ContextualParameters}, including the {@code EllipsoidToCentricTransform}'s one.</p>
     *
     * @param  transforms  the full chain of concatenated transforms.
     */
    public static void asDatumShift(final List<Object> transforms) {
        for (int i=transforms.size() - 2; --i >= 0;) {
            if (isOperation(GeographicToGeocentric.NAME, transforms.get(i)) &&
                isOperation(GeocentricToGeographic.NAME, transforms.get(i+2)))
            {
                final Object step = transforms.get(i+1);
                if (step instanceof LinearTransform) {
                    final BursaWolfParameters parameters = new BursaWolfParameters(null, null);
                    try {
                        parameters.setPositionVectorTransformation(((LinearTransform) step).getMatrix(), BURSAWOLF_TOLERANCE);
                    } catch (IllegalArgumentException e) {
                        /*
                         * Should not occur, except sometimes on inverse transform of relatively complex datum shifts
                         * (more than just translation terms). We can fallback on formatting the full matrix.
                         */
                        Logging.recoverableException(WKTUtilities.LOGGER, GeocentricAffine.class, "asDatumShift", e);
                        continue;
                    }
                    final boolean isTranslation = parameters.isTranslation();
                    final Parameters values = createParameters(isTranslation ? GeocentricTranslation.PARAMETERS
                                            : PositionVector7Param.PARAMETERS, parameters, isTranslation);
                    transforms.set(i+1, new FormattableObject() {
                        @Override protected String formatTo(final Formatter formatter) {
                            WKTUtilities.appendParamMT(values, formatter);
                            return WKTKeywords.Param_MT;
                        }
                    });
                }
            }
        }
    }

    /**
     * Returns {@code true} if the given object is an operation of the given name.
     */
    private static boolean isOperation(final String expected, final Object actual) {
        return (actual instanceof Parameterized) &&
               IdentifiedObjects.isHeuristicMatchForName(((Parameterized) actual).getParameterDescriptors(), expected);
    }
}
