/**
 * BasicNurbsCurve -- A basic implementation of NURBS curve.
 *
 * Copyright (C) 2009-2015 by Joseph A. Huwaldt. All rights reserved.
 *
 * This library is free software; you can redistribute it and/or modify it under the terms
 * of the GNU Lesser General Public License as published by the Free Software Foundation;
 * either version 2.1 of the License, or (at your option) any later version.
 *
 * This library is distributed in the hope that it will be useful, but WITHOUT ANY
 * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A
 * PARTICULAR PURPOSE. See the GNU Library General Public License for more details.
 *
 * You should have received a copy of the GNU Lesser General Public License along with
 * this program; if not, write to the Free Software Foundation, Inc., 59 Temple Place -
 * Suite 330, Boston, MA 02111-1307, USA. Or visit: http://www.gnu.org/licenses/lgpl.html
 *
 * This source file is based on, but heavily modified from, the LGPL jgeom (Geometry
 * Library for Java) code by Samuel Gerber, Copyright (C) 2005.
 */
package geomss.geom.nurbs;

import geomss.geom.Point;
import geomss.geom.Vector;
import jahuwaldt.js.math.BinomialCoef;
import jahuwaldt.tools.math.MathTools;
import java.text.MessageFormat;
import java.util.List;
import java.util.Objects;
import static java.util.Objects.isNull;
import static java.util.Objects.nonNull;
import javax.measure.converter.ConversionException;
import javax.measure.quantity.Length;
import javax.measure.unit.Unit;
import javolution.context.ArrayFactory;
import javolution.context.ObjectFactory;
import javolution.context.StackContext;
import javolution.lang.MathLib;
import javolution.lang.ValueType;
import javolution.text.Text;
import javolution.text.TextBuilder;
import javolution.util.FastTable;
import javolution.xml.XMLFormat;
import javolution.xml.stream.XMLStreamException;

/**
 * A basic implementation of a parametric NURBS curve.
 * 
 * <p> Modified by: Joseph A. Huwaldt </p>
 *
 * @author Samuel Gerber, Date: May 14, 2009, Version 1.0.
 * @version November 28, 2015
 */
@SuppressWarnings({"serial", "CloneableImplementsClone"})
public final class BasicNurbsCurve extends NurbsCurve implements ValueType {

    private static final double TOL = 1.0e-14;          //  Toleranace for avoiding 0/0 in binomial.

    private FastTable<ControlPoint> _cpoly;
    private KnotVector _uKnots;

    /*
     *  References:
     *  1.) Michael E. Mortenson, "Geometric Modeling, Third Edition", ISBN:9780831132989, 2008.
     *  2.) Piegl, L., Tiller, W., The Nurbs Book, 2nd Edition, Springer-Verlag, Berlin, 1997.
     */
    
    /**
     * Create a NURBS curve from the given control points, knots and degree.
     *
     * @param cps    Array of control points.  May not be null.
     * @param degree Degree of the NURBS curve
     * @param uK     Knot values.  May not be null.
     * @return A new NURBS curve
     * @throws IllegalArgumentException if the knot vector is not valid.
     */
    public static BasicNurbsCurve newInstance(ControlPoint cps[], int degree, double... uK) {
        if (cps.length < 1)
            throw new IllegalArgumentException(RESOURCES.getString("noControlPointsErr"));

        if (uK.length != degree + cps.length + 1)
            throw new IllegalArgumentException(
                    MessageFormat.format(RESOURCES.getString("crvWrongNumKnotsErr"),
                            uK.length, degree + cps.length + 1));

        //  Convert the units of the input control points.
        Unit<Length> refUnit = cps[0].getUnit();
        FastTable<ControlPoint> cpList = FastTable.newInstance();
        int length = cps.length;
        for (int i = 0; i < length; ++i) {
            cpList.add(cps[i].to(refUnit));
        }
        KnotVector knots = KnotVector.newInstance(degree, uK);

        //  Build the object.
        BasicNurbsCurve o = FACTORY.object();
        o._cpoly = cpList;
        o._uKnots = knots;

        return o;
    }

    /**
     * Generate a NURBS curve from the given control points and the given knot vector.
     *
     * @param cpList List of control points.  May not be null.
     * @param uKnots Knot vector for the curve.  May not be null.
     * @return A new NURBS curve.
     * @throws IllegalArgumentException if the knot vector is not consistent with the
     * number of control points.
     */
    public static BasicNurbsCurve newInstance(List<ControlPoint> cpList, KnotVector uKnots) {
        if (cpList.size() < 1)
            throw new IllegalArgumentException(RESOURCES.getString("noControlPointsErr"));

        if (uKnots.length() != uKnots.getDegree() + cpList.size() + 1)
            throw new IllegalArgumentException(
                    MessageFormat.format(RESOURCES.getString("crvWrongNumKnotsErr"),
                            uKnots.length(), uKnots.getDegree() + cpList.size() + 1));

        //  Convert the units of the input control points.
        Unit<Length> refUnit = cpList.get(0).getUnit();
        FastTable<ControlPoint> nCpList = FastTable.newInstance();
        int size = cpList.size();
        for (int i = 0; i < size; ++i) {
            ControlPoint cp = cpList.get(i);
            nCpList.add(cp.to(refUnit));
        }

        //  Build the object.
        BasicNurbsCurve o = FACTORY.object();
        o._cpoly = nCpList;
        o._uKnots = uKnots;

        return o;
    }

    /**
     * Returns the number of child-elements that make up this geometry element. This
     * implementation returns the number of control points in this NURBS curve.
     *
     * @return The number of control points in this curve.
     */
    @Override
    public int size() {
        return _cpoly.size();
    }

    /**
     * Returns the number of physical dimensions of the geometry element.
     *
     * @return The number of physical dimensions of this geometry element.
     */
    @Override
    public int getPhyDimension() {
        return _cpoly.get(0).getPhyDimension();
    }

    /**
     * Return a list of control points for this curve.
     *
     * @return the ordered control points
     */
    @Override
    public List<ControlPoint> getControlPoints() {
        List<ControlPoint> list = FastTable.newInstance();
        list.addAll(_cpoly);
        return list;
    }

    /**
     * Return the knot vector of this curve.
     *
     * @return The knot vector.
     */
    @Override
    public KnotVector getKnotVector() {
        return _uKnots;
    }

    /**
     * Return the degree of the NURBS curve.
     *
     * @return degree of curve
     */
    @Override
    public int getDegree() {
        return _uKnots.getDegree();
    }

    /**
     * Calculate a point on the curve for the given parametric distance along the curve,
     * <code>p(s)</code>.
     *
     * @param s parametric distance to calculate a point for (0.0 to 1.0 inclusive).
     * @return the calculated point
     */
    @Override
    public Point getRealPoint(double s) {
        validateParameter(s);

        //  End-point interpolation: C(0) = P0, and C(1) = Pn.
        if (parNearStart(s, TOL_S))         //  if (s == 0)
            return _cpoly.get(0).getPoint();
        else if (parNearEnd(s, TOL_S))      //  if (s == 1)
            return _cpoly.get(_cpoly.size() - 1).getPoint();

        StackContext.enter();
        try {
            //  Algorithm: A4.1, Ref. 2.
            /*
             This routine calculates the following:
             P(s) = sum_{i=1}^{N} ( Pi*Wi*Nik(s) ) / sum_{i=1}^{N} ( Wi*Nik(s) )
             */
            Unit<Length> unit = getUnit();
            int span = _uKnots.findSpan(s);
            int degree = _uKnots.getDegree();
            int sind = span - degree;

            //  Get the basis function values (Nik).
            double[] Nk = _uKnots.basisFunctions(span, s);

            //  Apply the basis functions to each degree of the curve.
            //  Calculate "cw" point as sum_{i=1}^N (Pi*Wi*Nik(s)) and "cw" weight as sum_{i=1}^N (Wi*Ni,k(s)).
            ControlPoint cw = ControlPoint.newInstance(getPhyDimension(), unit);
            for (int i = 0; i <= degree; i++) {
                ControlPoint tmp = _cpoly.get(sind + i).applyWeight();
                tmp = tmp.times(Nk[i]);
                cw = cw.plus(tmp);
            }

            //  Convert the control point to a geometric point.
            //  The "cw" weight is the sum of Wi*Ni,k calculated above.
            Point outPoint = cw.getPoint();
            if (outPoint.normValue() > MathTools.SQRT_EPS)
                outPoint = outPoint.divide(cw.getWeight());
            outPoint = outPoint.to(unit);

            return StackContext.outerCopy(outPoint);
        } finally {
            StackContext.exit();
        }
    }

    /**
     * Calculate all the derivatives from <code>0</code> to <code>grade</code> with
     * respect to parametric distance on the curve for the given parametric distance along
     * the curve, <code>d^{grade}p(s)/d^{grade}s</code>.
     * <p>
     * Example:<br>
     * 1st derivative (grade = 1), this returns <code>[p(s), dp(s)/ds]</code>;<br>
     * 2nd derivative (grade = 2), this returns <code>[p(s), dp(s)/ds, d^2p(s)/d^2s]</code>; etc.
     * </p>
     *
     * @param s     Parametric distance to calculate derivatives for (0.0 to 1.0
     *              inclusive).
     * @param grade The maximum grade to calculate the derivatives for (1=1st derivative,
     *              2=2nd derivative, etc)
     * @return A list of derivatives up to the specified grade of the curve at the
     *         specified parametric position.
     * @throws IllegalArgumentException if the grade is &lt; 0.
     */
    @Override
    public List<Vector<Length>> getSDerivatives(double s, int grade) {
        validateParameter(s);
        s = roundParNearEnds(s);

        if (grade < 0)
            throw new IllegalArgumentException(RESOURCES.getString("gradeLTZeroErr"));

        //  Uses Algorithm: A4.2, Ref. 2, pg. 127
        //  First find unweighted derivatives of a B-Spline
        //  Algorithm: A3.2, Ref. 2, pg. 93.
        int dim = getPhyDimension();
        int degree = _uKnots.getDegree();
        int span = _uKnots.findSpan(s);
        int sind = span - degree;
        int ds = MathLib.min(grade, degree);
        int gradeP1 = grade + 1;

        Unit<Length> unit = getUnit();
        Vector[] cOut = Vector.allocateArray(gradeP1);   //  Created outside of StackContext.
        StackContext.enter();
        try {
            //  Get the basis function derivative values for a B-Spline (unweighted derivatives).
            double[][] Nik = _uKnots.basisFunctionDerivatives(span, s, ds);

            //  Calculate z = sum_i^N ( Wi*Nik ) for each basis function derivative and
            //  p = sum_i^N (Pi*Wi*Nik) for each control point and basis function derivative.
            int order = degree + 1;
            double[] z = ArrayFactory.DOUBLES_FACTORY.array(gradeP1);
            Point[] p = Point.allocateArray(gradeP1);
            for (int k = order; k <= grade; ++k) {
                z[k] = 0.;
                p[k] = Point.newInstance(dim, unit);
            }
            for (int k = 0; k <= ds; ++k) {
                z[k] = 0.;
                p[k] = Point.newInstance(dim, unit);
                for (int i = 0; i < order; i++) {
                    ControlPoint cpi = _cpoly.get(sind + i);
                    double Wi = cpi.getWeight();
                    double WiNik = Wi * Nik[k][i];
                    z[k] += WiNik;
                    p[k] = p[k].plus(cpi.getPoint().times(WiNik));
                }
            }

            //  Compute the weighted derivatives by stepping through a binomial expansion for each derivative.
            //    Algorithm: A4.2, Ref. 2, pg. 127      
            Point[] c = binomial(p, z, gradeP1);

            //  Copy the final points out of the StackContext.
            for (int i = gradeP1-1; i >= 0; --i) {
                cOut[i] = StackContext.outerCopy(c[i].to(unit).toGeomVector());
            }

        } finally {
            StackContext.exit();
        }

        //  Convert to a list of Vector objets.
        FastTable<Vector<Length>> output = FastTable.newInstance();
        Point origin = Point.valueOf(cOut[0]);
        for (int i = 0; i <= grade; ++i) {
            Vector<Length> v = cOut[i];
            v.setOrigin(origin);
            output.add(v);
        }

        Vector.recycleArray(cOut);

        return output;
    }

    /**
     * Compute the weighted derivatives by stepping through a binomial expansion for each
     * derivative. e.g.:
     * <pre>
     *   surface point, c(0) = p(0)/z(0)
     *      1st derivative, c(1) = (p(1) - c(0)*z(1))/z(0),
     *      2nd derivative, c(2) = (p(2) - c(0)*z(2) - 2*c(1)*z(1))/z(0),
     *      3rd derivative, c(3) = (p(3) - c(0)*z(3) - 3*c(1)*z(2) - 3*c(2)*z(1))/z(0),
     *      etc
     * </pre>
     *
     * @param pnts    The list of weighted control points: p[0..grade+1].
     * @param z       The derivative weights (Wi*Nik): z[0..grade+1].
     * @param gradeP1 The maximum grade to calculate the v-derivatives for (1=1st
     *                derivative, 2=2nd derivative, etc) + 1.
     */
    private static Point[] binomial(Point[] pnts, double[] z, int gradeP1) {

        BinomialCoef bin = BinomialCoef.newInstance(gradeP1);   //  Get the binomial coefficients.

        Point[] c = Point.allocateArray(gradeP1);
        Point ZERO = Point.newInstance(pnts[0].getPhyDimension(), pnts[0].getUnit());
        double z0inv = 1.0 / z[0];
        for (int k = 0; k < gradeP1; ++k) {
            Point v = pnts[k];
            for (int i = k; i > 0; --i) {
                v = v.minus(c[k - i].times(bin.get(k, i) * z[i]));
            }
            if (v.normValue() > TOL)        //  Handles 0/0 == 0.
                c[k] = v.times(z0inv);
            else
                c[k] = ZERO;  //  v.times(0);
        }

        return c;
    }

    /**
     * Return the coordinate point representing the minimum bounding box corner of this
     * geometry element (e.g.: min X, min Y, min Z).
     *
     * @return The minimum bounding box coordinate for this geometry element.
     * @throws IndexOutOfBoundsException if this list contains no geometry.
     */
    @Override
    public Point getBoundsMin() {

        Point minPoint = _cpoly.get(0).getBoundsMin();

        int size = _cpoly.size();
        for (int i = 0; i < size; ++i) {
            ControlPoint cp = _cpoly.get(i);
            minPoint = minPoint.min(cp.getBoundsMin());
        }

        return minPoint;
    }

    /**
     * Return the coordinate point representing the maximum bounding box corner (e.g.: max
     * X, max Y, max Z).
     *
     * @return The maximum bounding box coordinate for this geometry element.
     * @throws IndexOutOfBoundsException if this list contains no elements.
     */
    @Override
    public Point getBoundsMax() {

        Point maxPoint = _cpoly.get(0).getBoundsMax();

        int size = _cpoly.size();
        for (int i = 0; i < size; ++i) {
            ControlPoint cp = _cpoly.get(i);
            maxPoint = maxPoint.max(cp.getBoundsMax());
        }

        return maxPoint;
    }

    /**
     * Returns the unit in which the control points in this curve are stated.
     *
     * @return The unit that this curves points are stated in.
     */
    @Override
    public Unit<Length> getUnit() {
        return _cpoly.get(0).getPoint().getUnit();
    }

    /**
     * Returns the equivalent to this curve but stated in the specified unit.
     *
     * @param unit The length unit of the curve to be returned. May not be null.
     * @return An equivalent to this curve but stated in the specified unit.
     * @throws ConversionException if the the input unit is not a length unit.
     */
    @Override
    public BasicNurbsCurve to(Unit<Length> unit) throws ConversionException {
        if (unit.equals(getUnit()))
            return this;

        //  Convert the control points.
        FastTable<ControlPoint> nCpoly = FastTable.newInstance();
        int size = _cpoly.size();
        for (int i = 0; i < size; ++i) {
            ControlPoint cp = _cpoly.get(i);
            nCpoly.add(cp.to(unit));
        }

        BasicNurbsCurve curve = BasicNurbsCurve.newInstance(nCpoly, getKnotVector());
        return copyState(curve);    //  Copy over the super-class state for this object to the new one.
    }

    /**
     * Return the equivalent of this curve converted to the specified number of physical
     * dimensions. If the number of dimensions is greater than this element, then zeros
     * are added to the additional dimensions. If the number of dimensions is less than
     * this element, then the extra dimensions are simply dropped (truncated). If the new
     * dimensions are the same as the dimension of this element, then this element is
     * simply returned.
     *
     * @param newDim The dimension of the curve to return.
     * @return The equivalent of this curve converted to the new dimensions.
     */
    @Override
    public BasicNurbsCurve toDimension(int newDim) {
        if (getPhyDimension() == newDim)
            return this;

        FastTable<ControlPoint> newCpLst = FastTable.newInstance();
        int size = _cpoly.size();
        for (int i = 0; i < size; ++i) {
            ControlPoint cp = _cpoly.get(i);
            Point pnt = cp.getPoint().toDimension(newDim);
            ControlPoint newCp = ControlPoint.valueOf(pnt, cp.getWeight());
            newCpLst.add(newCp);
        }

        BasicNurbsCurve newCrv = BasicNurbsCurve.newInstance(newCpLst, getKnotVector());

        return copyState(newCrv);   //  Copy over the super-class state for this object to the new one.
    }

    /**
     * Returns a copy of this {@link BasicNurbsCurve} instance
     * {@link javolution.context.AllocatorContext allocated} by the calling thread
     * (possibly on the stack).
     *
     * @return an identical and independent copy of this object.
     */
    @Override
    public BasicNurbsCurve copy() {
        return copyOf(this);
    }

    /**
     * Return a copy of this object with any transformations or subranges removed
     * (applied).
     *
     * @return A copy of this object with any transformations or subranges removed
     *         (applied).
     */
    @Override
    public BasicNurbsCurve copyToReal() {
        return copy();
    }

    /**
     * Return an immutable version of this NURBS curve. This implementation simply returns
     * this BasicNurbsCurve instance.
     *
     * @return an immutable version of this curve.
     */
    @Override
    public BasicNurbsCurve immutable() {
        return this;
    }
    
    /**
     * Resets the internal state of this object to its default values.
     */
    @Override
    public void reset() {
        _cpoly.reset();
        super.reset();
    }

    /**
     * Compares this BasicNurbsCurve against the specified object for strict equality
     * (same values and same units).
     *
     * @param obj the object to compare with.
     * @return <code>true</code> if this object is identical to that object;
     *         <code>false</code> otherwise.
     */
    @Override
    public boolean equals(Object obj) {
        if (this == obj)
            return true;
        if ((obj == null) || (obj.getClass() != this.getClass()))
            return false;

        BasicNurbsCurve that = (BasicNurbsCurve)obj;
        return this._cpoly.equals(that._cpoly)
                && this._uKnots.equals(that._uKnots)
                && super.equals(obj);
    }

    /**
     * Returns the hash code for this parameter.
     *
     * @return the hash code value.
     */
    @Override
    public int hashCode() {
        return 31*super.hashCode() + Objects.hash(_cpoly, _uKnots);
    }

    /**
     * Returns the text representation of this geometry element that consists of the name
     * followed by the control point values, followed by the knot vector. For example:
     * <pre>
     *   {aCurve = {{{1 ft, 0 ft}, 1.0}, {{0 ft, 1 ft}, 0.25}, {{-1 ft, 0 ft}, 1.0}},{degree=2,{0.0, 0.0, 0.0, 1.0, 1.0, 1.0}}}
     * </pre>
     * If there is no name, then the output looks like this:
     * <pre>
     *   {{{{1 ft, 0 ft}, 1.0}, {{0 ft, 1 ft}, 0.25}, {{-1 ft, 0 ft}, 1.0}},{degree=2,{0.0, 0.0, 0.0, 1.0, 1.0, 1.0}}}
     * </pre>
     *
     * @return the text representation of this geometry element.
     */
    @Override
    public Text toText() {
        TextBuilder tmp = TextBuilder.newInstance();
        tmp.append('{');
        String nameStr = getName();
        boolean hasName = nonNull(nameStr);
        if (hasName) {
            tmp.append(nameStr);
            tmp.append(" = {");
        }
        int size = _cpoly.size();
        for (int i = 0; i < size; i++) {
            tmp.append(_cpoly.get(i).toText());
            if (i != size - 1) {
                tmp.append(", ");
            }
        }
        tmp.append("},");

        tmp.append(_uKnots.toText());

        if (hasName)
            tmp.append('}');
        tmp.append('}');
        Text txt = tmp.toText();
        TextBuilder.recycle(tmp);
        return txt;
    }

    /**
     * Recycles a <code>BasicNurbsCurve</code> instance immediately (on the stack when
     * executing in a StackContext).
     *
     * @param instance The instance to be recycled.
     */
    public static void recycle(BasicNurbsCurve instance) {
        FACTORY.recycle(instance);
    }

    /**
     * Holds the default XML representation for this object.
     */
    @SuppressWarnings("FieldNameHidesFieldInSuperclass")
    protected static final XMLFormat<BasicNurbsCurve> XML = new XMLFormat<BasicNurbsCurve>(BasicNurbsCurve.class) {
        @Override
        public BasicNurbsCurve newInstance(Class<BasicNurbsCurve> cls, InputElement xml) throws XMLStreamException {
            return FACTORY.object();
        }

        @SuppressWarnings("unchecked")
        @Override
        public void read(InputElement xml, BasicNurbsCurve obj) throws XMLStreamException {
            NurbsCurve.XML.read(xml, obj);     // Call parent read.

            FastTable<ControlPoint> cpList = xml.get("CPoly", FastTable.class);
            KnotVector knots = xml.get("Knots", KnotVector.class);

            if (isNull(cpList) || cpList.size() < 1)
                throw new XMLStreamException(RESOURCES.getString("noControlPointsErr"));

            if (knots.length() != knots.getDegree() + cpList.size() + 1)
                throw new XMLStreamException(
                        MessageFormat.format(RESOURCES.getString("crvWrongNumKnotsErr"),
                                knots.length(), knots.getDegree() + cpList.size() + 1));

            //  Convert the units of the input control points.
            Unit<Length> refUnit = cpList.get(0).getUnit();
            FastTable<ControlPoint> nCpList = FastTable.newInstance();
            int size = cpList.size();
            for (int i = 0; i < size; ++i) {
                ControlPoint cp = cpList.get(i);
                nCpList.add(cp.to(refUnit));
            }
            obj._cpoly = nCpList;
            obj._uKnots = knots;

        }

        @Override
        public void write(BasicNurbsCurve obj, OutputElement xml) throws XMLStreamException {
            NurbsCurve.XML.write(obj, xml);    // Call parent write.

            xml.add(obj._cpoly, "CPoly", FastTable.class);
            xml.add(obj._uKnots, "Knots", KnotVector.class);

        }
    };

    ///////////////////////
    // Factory creation. //
    ///////////////////////
    private BasicNurbsCurve() { }

    @SuppressWarnings("unchecked")
    private static final ObjectFactory<BasicNurbsCurve> FACTORY = new ObjectFactory<BasicNurbsCurve>() {
        @Override
        protected BasicNurbsCurve create() {
            return new BasicNurbsCurve();
        }

        @Override
        protected void cleanup(BasicNurbsCurve obj) {
            obj.reset();
            obj._cpoly = null;
            obj._uKnots = null;
        }
    };

    @SuppressWarnings("unchecked")
    private static BasicNurbsCurve copyOf(BasicNurbsCurve original) {
        BasicNurbsCurve o = FACTORY.object();
        FastTable<ControlPoint> newCPList = FastTable.newInstance();
        int size = original._cpoly.size();
        for (int i = 0; i < size; ++i) {
            ControlPoint cp = original._cpoly.get(i);
            newCPList.add(cp.copy());
        }
        o._cpoly = newCPList;
        o._uKnots = original._uKnots.copy();
        original.copyState(o);
        return o;
    }

}