/**
 * NurbsCurve -- Implementation and interface in common to all NURBS curves.
 *
 * Copyright (C) 2009-2025, 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.*;
import jahuwaldt.js.param.Parameter;
import jahuwaldt.tools.math.BracketRoot1D;
import jahuwaldt.tools.math.Evaluatable1D;
import jahuwaldt.tools.math.MathTools;
import jahuwaldt.tools.math.RootException;
import java.text.MessageFormat;
import java.util.ArrayList;
import java.util.Collections;
import java.util.List;
import static java.util.Objects.nonNull;
import static java.util.Objects.requireNonNull;
import java.util.logging.Level;
import java.util.logging.Logger;
import javax.measure.quantity.Area;
import javax.measure.quantity.Dimensionless;
import javax.measure.quantity.Length;
import javax.measure.unit.Unit;
import javolution.context.ObjectFactory;
import javolution.context.StackContext;
import javolution.lang.MathLib;
import javolution.text.Text;
import javolution.text.TextBuilder;
import javolution.util.FastList;
import javolution.util.FastMap;
import javolution.util.FastTable;
import javolution.xml.XMLFormat;
import javolution.xml.stream.XMLStreamException;
import org.jscience.mathematics.number.Float64;
import org.jscience.mathematics.vector.Float64Vector;

/**
 * The interface and implementation in common to all NURBS curves.
 *
 * <p> Modified by: Joseph A. Huwaldt </p>
 *
 * @author Joseph A. Huwaldt, Date: May 28, 2009
 * @version February 17, 2025
 */
@SuppressWarnings({"serial", "CloneableImplementsClone"})
public abstract class NurbsCurve extends AbstractCurve<NurbsCurve> {

    /**
     * Return a list of control points for this curve.
     *
     * @return the ordered control points
     */
    public abstract List<ControlPoint> getControlPoints();

    /**
     * Return the knot vector of this curve.
     *
     * @return The knot vector.
     */
    public abstract KnotVector getKnotVector();

    /**
     * Return the degree of the NURBS curve.
     *
     * @return degree of curve
     */
    public abstract int getDegree();

    /**
     * Return an immutable version of this NURBS curve.
     *
     * @return an immutable version of this curve.
     */
    public abstract BasicNurbsCurve immutable();
    
    /**
     * Return a new curve that is identical to this one, but with the parameterization
     * reversed.
     *
     * @return A new curve that is identical to this one but with the parameterization
     *         reversed.
     */
    @Override
    public NurbsCurve reverse() {
        List<ControlPoint> oldCPs = getControlPoints();

        //  Reverse the order of the control points.
        FastTable<ControlPoint> cps = FastTable.newInstance();
        for (int i = oldCPs.size() - 1; i >= 0; --i)
            cps.add(oldCPs.get(i));
        FastTable.recycle((FastTable)oldCPs);

        //  Reverse the order of the knot vector.
        KnotVector kv = getKnotVector().reverse();

        //  Return a new curve.
        BasicNurbsCurve curve = BasicNurbsCurve.newInstance(cps, kv);
        copyState(curve);   //  Copy over the super-class state for this curve to the new one.

        FastTable.recycle(cps);
        return curve;
    }

    /**
     * Create and return a new {@link BasicNurbsCurve NURBS curve} that is geometrically
     * identical to this one but with a new knot inserted at the specified parametric
     * location. The parameterization of the new curve is identical to this curve, but the
     * new curve will have a new knot at the specified location and a new set of control
     * points.
     *
     * @param s             The parametric location where the knot is to be inserted.
     * @param numInsertions The number of times that the knot should be inserted at the
     *                      specified position.
     * @return A new NURBS curve that is identical to this one, but with a new knot
     *         inserted at the specified location.
     */
    public BasicNurbsCurve insertKnot(double s, int numInsertions) {
        //  Algorithm A5.1, pg. 151, Piegl, L., Tiller, W., The Nurbs Book, 2nd Edition, Springer-Verlag, Berlin, 1997.

        validateParameter(s);
        if (numInsertions < 1)
            throw new IllegalArgumentException(RESOURCES.getString("numInsertionsLE0Err"));

        StackContext.enter();
        try {
            //  Construct the new knot vector by inserting the specified parametric position
            //  into the list of knots.
            KnotVector SP = getKnotVector();
            int numKnots = SP.length();
            int p = SP.getDegree();
            int order = p + 1;
            int r = numInsertions;

            //  Find where in the knot vector the new knot lies.
            int k = 0;
            for (int i = 0; i < numKnots; ++i) {
                if (SP.getValue(i) > s) {
                    k = i - 1;
                    break;
                }
            }

            //  Determine the multiplicity of this knot.
            int sm = 0;
            if (s <= SP.getValue(k))
                sm = SP.findMultiplicity(k);

            if ((r + sm) > order)
                r = order - sm;

            if (r <= 0)
                return StackContext.outerCopy((BasicNurbsCurve)this.copyToReal());

            //  Load the new knot vector.
            FastTable<Float64> SQList = FastTable.newInstance();
            for (int i = 0; i <= k; ++i) {
                SQList.add(SP.get(i));
            }

            for (int i = 0; i < r; ++i) {
                SQList.add(Float64.valueOf(s));
            }

            for (int i = k + 1; i < numKnots; ++i) {
                SQList.add(SP.get(i));
            }

            KnotVector SQ = KnotVector.newInstance(p, Float64Vector.valueOf(SQList));

            //  Save unaltered control points.
            List<ControlPoint> Pw = getControlPoints();
            int np = Pw.size();
            FastTable<ControlPoint> Qw = FastTable.newInstance();
            Qw.setSize(np + r);

            for (int i = 0; i <= k - p; ++i) {
                Qw.set(i, Pw.get(i));
            }

            for (int i = k - sm; i < np; ++i) {
                Qw.set(i + r, Pw.get(i));
            }

            //  Temporary storage for altered control points.
            FastTable<ControlPoint> Rw = FastTable.newInstance();
            Rw.setSize(order);
            for (int i = 0; i <= p - sm; ++i) {
                Rw.set(i, Pw.get(k - p + i).applyWeight());     //  Project them from NURBS to B-Spline.
            }

            //  Insert the knot r times.
            int L = 0;
            for (int j = 1; j <= r; ++j) {
                L = k - p + j;
                int pmjms = p - j - sm;
                for (int i = 0; i <= pmjms; ++i) {
                    //  Defining new control points in terms of old.
                    //  alpha = (s - SP(L+i))/(SP(i+k+1) - SP(L+i))
                    //  Rw(i) = alpha*Rw[i+1] + (1.0 - alpha)*Rw[i]
                    double alpha = (s - SP.getValue(L + i)) / (SP.getValue(i + k + 1) - SP.getValue(L + i));

                    ControlPoint Rwi = Rw.get(i);
                    Rwi = Rwi.times(1.0 - alpha);
                    ControlPoint Rwip1 = Rw.get(i + 1);
                    Rwip1 = Rwip1.times(alpha);
                    Rwip1 = Rwip1.plus(Rwi);

                    //  Store the new control point.
                    Rw.set(i, Rwip1);
                }

                //  Store the new control points (converting back from B-Spline to NURBS).
                Qw.set(L, Rw.get(0).getHomogeneous());
                if (pmjms > 0)
                    Qw.set(k + r - j - sm, Rw.get(pmjms).getHomogeneous());

            }

            //  Load the remaining control points.
            for (int i = L + 1; i < k - sm; ++i) {
                Qw.set(i, Rw.get(i - L).getHomogeneous());  //  Converting back from B-Spline to NURBS.
            }

            //  Construct and return a new curve.
            BasicNurbsCurve curve = BasicNurbsCurve.newInstance(Qw, SQ);
            return StackContext.outerCopy(curve);

        } finally {
            StackContext.exit();
        }

    }

    /**
     * Create and return a new {@link BasicNurbsCurve NURBS curve} that is geometrically
     * identical to this one but with the specified list of knots inserted into it's knot
     * vector. The parameterization of the new curve is identical to this curve, but the
     * new curve will have the new knots at the specified locations and a new set of
     * control points. This is more efficient than repeatedly calling insertKnot() for
     * multiple knot insertions.
     *
     * @param newKnots The parametric locations where the knots are to be inserted.
     * @return A new NURBS curve that is identical to this one, but with the new knots
     *         inserted at the specified locations.
     */
    public BasicNurbsCurve refineKnotVector(double[] newKnots) {
        return refineKnotVector(Float64Vector.valueOf(newKnots));
    }

    /**
     * Create and return a new {@link BasicNurbsCurve NURBS curve} that is geometrically
     * identical to this one but with the specified list of knots inserted into it's knot
     * vector. The parameterization of the new curve is identical to this curve, but the
     * new curve will have the new knots at the specified locations and a new set of
     * control points. This is more efficient than repeatedly calling insertKnot() for
     * multiple knot insertions.
     *
     * @param newKnots The parametric locations where the knots are to be inserted. May
     *                 not be null.
     * @return A new NURBS curve that is identical to this one, but with the new knots
     *         inserted at the specified locations.
     */
    public BasicNurbsCurve refineKnotVector(Float64Vector newKnots) {
        //  Algorithm A5.4, pg. 164, Piegl, L., Tiller, W., The Nurbs Book, 2nd Edition, Springer-Verlag, Berlin, 1997.

        //  Check on the validity of the new knots.
        int r = newKnots.getDimension() - 1;
        double sOld = 0;
        for (int i = 0; i <= r; ++i) {
            double s = newKnots.getValue(i);
            validateParameter(s);
            if (s < sOld)
                throw new IllegalArgumentException(RESOURCES.getString("knotsNotIncreasingErr"));
            sOld = s;
        }

        StackContext.enter();
        try {
            //  Get the existing knot vector.
            KnotVector U = getKnotVector();
            int p = U.getDegree();

            //  Get the existing control points.
            List<ControlPoint> Pw = getControlPoints();
            int n = Pw.size() - 1;

            int m = n + p + 1;

            //  Allocate storage for the new control points.
            FastTable<ControlPoint> Qw = FastTable.newInstance();
            Qw.setSize(n + r + 2);

            //  Allocate storage for the new knot vector.
            FastTable<Float64> UbarList = FastTable.newInstance();
            UbarList.setSize(m + r + 2);

            int a = U.findSpan(newKnots.getValue(0));
            int b = U.findSpan(newKnots.getValue(r));
            ++b;

            //  Copy over the unaltered control points.
            for (int j = 0; j <= a - p; ++j) {
                Qw.set(j, Pw.get(j));
            }
            for (int j = b - 1; j <= n; ++j) {
                Qw.set(j + r + 1, Pw.get(j));
            }

            //  Copy over the unaltered knots.
            for (int j = 0; j <= a; ++j) {
                UbarList.set(j, U.get(j));
            }
            for (int j = b + p; j <= m; ++j) {
                UbarList.set(j + r + 1, U.get(j));
            }

            int i = b + p - 1;
            int k = b + p + r;
            for (int j = r; j >= 0; --j) {
                while (newKnots.getValue(j) <= U.getValue(i) && i > a) {
                    Qw.set(k - p - 1, Pw.get(i - p - 1));
                    UbarList.set(k, U.get(i));
                    --k;
                    --i;
                }

                Qw.set(k - p - 1, Qw.get(k - p));

                for (int l = 1; l <= p; ++l) {
                    int ind = k - p + l;
                    int indm1 = ind - 1;
                    int kpl = k + l;
                    double Ubarkpl = UbarList.get(kpl).doubleValue();
                    double alpha = Ubarkpl - newKnots.getValue(j);
                    if (MathTools.isApproxZero(alpha))
                        Qw.set(indm1, Qw.get(ind));
                    else {
                        alpha /= Ubarkpl - U.getValue(i - p + l);

                        //  Qw[ind-1] = alpha*Qw[ind-1] + (1.0-alpha)*Qw[ind];
                        ControlPoint Qwind = Qw.get(ind).applyWeight();
                        ControlPoint Qwindm1 = Qw.get(indm1).applyWeight();
                        Qwindm1 = Qwindm1.times(alpha);
                        Qwind = Qwind.times(1.0 - alpha);
                        Qwind = Qwind.plus(Qwindm1);
                        Qw.set(indm1, Qwind.getHomogeneous());
                    }
                }

                UbarList.set(k, newKnots.get(j));
                --k;
            }

            //  Construct and return a new curve.
            KnotVector Ubar = KnotVector.newInstance(p, Float64Vector.valueOf(UbarList));
            BasicNurbsCurve curve = BasicNurbsCurve.newInstance(Qw, Ubar);
            return StackContext.outerCopy(curve);

        } finally {
            StackContext.exit();
        }

    }

    /**
     * Return a {@link BasicNurbsCurve NURBS curve} that is geometrically identical to
     * this one but with the specified list of knots merged into it's knot vector. The
     * parameterization of the new curve is identical to this curve, but the new curve
     * will have any new knots at the specified locations added and a new set of control
     * points. If the input knot list matches the knot list of this curve, then this curve
     * is returned.
     *
     * @param Um The knot vector to merge with this curve's knot vector. May not be null.
     * @return A NURBS curve that is identical to this one, but with the merger of the
     *         input knot vector and this curve's knot vector.
     */
    public NurbsCurve mergeKnotVector(Float64Vector Um) {

        StackContext.enter();
        try {
            //  Get this curve's knot vector.
            KnotVector U = getKnotVector();
            int Usize = U.length();
            int Umsize = Um.getDimension();

            //  Find the knots that need inserting.
            FastTable<Float64> Iknots = FastTable.newInstance();

            boolean done = false;
            int ia = 0, ib = 0;
            while (!done) {
                double Umv = Um.getValue(ib);
                double Uv = U.getValue(ia);
                if (MathLib.abs(Umv - Uv) <= TOL_S) {   //  Umv == Uv
                    ++ib;
                    ++ia;
                } else {
                    Iknots.add(Float64.valueOf(Umv));
                    ++ib;
                }
                done = (ia >= Usize || ib >= Umsize);
            }

            //  Check to see if any knots need inserting.
            if (Iknots.isEmpty()) {
                return this;
            }

            BasicNurbsCurve newCrv = this.refineKnotVector(Float64Vector.valueOf(Iknots));
            return StackContext.outerCopy(newCrv);

        } finally {
            StackContext.exit();
        }

    }

    /**
     * Attempts to remove the knot with the specified index from this NURBS curve the
     * specified number of times. If the knot can not be removed without changing the
     * shape of the curve, then a reference to this curve is returned:
     * <code>if (crv == crv.removeKnot(...)) { no knots removed }</code>.
     *
     * @param index       The index of the knot to be removed (degree + 1 &le; index &lt;
     *                    knot vector length - degree).
     * @param numRemovals The number of times that the knot at "index" should be removed
     *                    (1 &le; num &le; multiplicity of the knot). If numRemovals is &gt;
     *                    multiplicity of the knot, then numRemovals is set to the
     *                    multiplicity of the knot.
     * @param tolerance   The knot(s) specified will be removed as long as the deviation
     *                    of the curve is everywhere less than this tolerance. May not be
     *                    null.
     * @return A new NURBS curve that is identical to this one (to within the specified
     *         tolerance), but with the specified knot removed or a reference to this
     *         curve if the knot could not be removed.
     */
    public NurbsCurve removeKnot(int index, int numRemovals, Parameter<Length> tolerance) {
        //  Algorithm: A5.8, pg. 185, Piegl, L., Tiller, W., The Nurbs Book, 2nd Edition, Springer-Verlag, Berlin, 1997.
        requireNonNull(tolerance);
        
        KnotVector Ukv = getKnotVector();
        int degree = Ukv.getDegree();
        int order = degree + 1;
        int m = Ukv.length();

        //  Check to make sure the knot to be removed is a legal one.
        if (index < order || index >= m - degree)
            throw new IllegalArgumentException(
                    MessageFormat.format(RESOURCES.getString("illegalKnotIndexForRemoval"),
                            index, order, m - degree));
        if (numRemovals < 0)
            throw new IllegalArgumentException(RESOURCES.getString("negNumKnotsRemoveErr"));
        if (numRemovals == 0)
            return this;

        StackContext.enter();
        try {
            //  Check if this is a multiple knot.
            //  If it is, reset "index" to be the largest index for the multiple knot.
            //  Want: Ukv(index) != Ukv(index+1)
            while (index < m - degree - 1 && Ukv.getValue(index + 1) == Ukv.getValue(index)) {
                ++index;
            }
            double u = Ukv.getValue(index);

            //  Determine multiplicity.
            int s = Ukv.findMultiplicity(index);
            if (numRemovals > s)
                numRemovals = s;

            //  Get the original control points.
            List<ControlPoint> Pw = getControlPoints();
            int n = Pw.size();

            //  Determine an appropriate and meaningful tolerance for knot removal.
            double tol = tolerance.getValue(getUnit());
            boolean isRational = false;
            double wmin = Double.MAX_VALUE;
            double Pmax = -Double.MAX_VALUE;
            for (int i = 0; i < n; ++i) {
                ControlPoint Pwi = Pw.get(i);
                double w = Pwi.getWeight();
                wmin = MathLib.min(wmin, w);
                if (MathTools.isApproxEqual(w, 1))
                    isRational = true;
                Pmax = MathLib.max(Pmax, Pwi.normValue());
            }
            if (isRational)
                tol = tol * (1 + Pmax) / wmin;

            //  Allocate storage for calculated control points.
            FastTable<ControlPoint> temp = FastTable.newInstance();
            temp.setSize(2 * degree + 1);

            //  First control point out.
            int fout = (2 * index - s - degree) / 2;
            int last = index - s;
            int first = index - degree;
            int i, j, t;

            //  This loop is Eqn. 5.28
            for (t = 0; t < numRemovals; ++t) {
                //  Diff. in index between temp and Pw.
                int off = first - 1;
                temp.set(0, Pw.get(off).applyWeight());
                ControlPoint cp = Pw.get(last + 1).applyWeight();
                temp.set(last + 1 - off, cp);
                i = first;
                j = last;
                int ii = 1;
                int jj = last - off;
                boolean remflag = false;

                while (j - i > t) {
                    //  Compute new control points for one removal step.
                    double alfi = (u - Ukv.getValue(i)) / (Ukv.getValue(i + order + t) - Ukv.getValue(i));
                    double alfj = (u - Ukv.getValue(j - t)) / (Ukv.getValue(j + order) - Ukv.getValue(j - t));
                    ControlPoint Pwi = Pw.get(i).applyWeight();
                    ControlPoint Pwj = Pw.get(j).applyWeight();

                    //  temp[ii] = (Pw[i] - (1.0-alfi)*temp[ii-1])/alfi;
                    ControlPoint tmpCP = temp.get(ii - 1);
                    tmpCP = tmpCP.times(1.0 - alfi);
                    Pwi = Pwi.minus(tmpCP);
                    Pwi = Pwi.divide(alfi);
                    temp.set(ii, Pwi);

                    //  temp[jj] = (Pw[j] - alfj*temp[jj+1])/(1 - alfi);
                    tmpCP = temp.get(jj + 1);
                    tmpCP = tmpCP.times(alfj);
                    Pwj = Pwj.minus(tmpCP);
                    Pwj = Pwj.divide(1.0 - alfj);
                    temp.set(jj, Pwj);

                    ++i;
                    ++ii;
                    --j;
                    --jj;
                }

                //  Check if knot removable.
                if (j - i < t) {
                    if (temp.get(ii - 1).getHomogeneous().distance(temp.get(jj + 1).getHomogeneous()) <= tol)
                        remflag = true;
                } else {
                    double alfi = (u - Ukv.getValue(i)) / (Ukv.getValue(i + order + t) - Ukv.getValue(i));
                    //  tmp = alfi*temp[ii+t+1] + (1.0-alfi)*temp[ii-1];
                    ControlPoint tmp2 = temp.get(ii - 1);
                    tmp2 = tmp2.times(1.0 - alfi);
                    ControlPoint tmp = temp.get(ii + t + 1);
                    tmp = tmp.times(alfi);
                    tmp = tmp.plus(tmp2);
                    if (Pw.get(i).distance(tmp.getHomogeneous()) <= tol)
                        remflag = true;
                }
                if (!remflag)
                    //  Can not remove any more knots.
                    break;  //  Get out of the for-loop.
                else {
                    //  Successfull removal.  Save new control points.
                    i = first;
                    j = last;

                    while (j - i > t) {
                        Pw.set(i, temp.get(i - off).getHomogeneous());
                        Pw.set(j, temp.get(j - off).getHomogeneous());
                        ++i;
                        --j;
                    }
                }
                --first;
                ++last;
            }   //  end of for-loop.

            if (t == 0) {
                //  No knots removed.
                return this;
            }

            //  Construct a new knot vector by shifting the old knots over.
            FastTable<Float64> newKV = FastTable.newInstance();
            int end = index + 1 - t;
            for (int k = 0; k < end; ++k) {
                newKV.add(Ukv.get(k));
            }
            for (int k = index + 1; k < m; ++k) {
                newKV.add(Ukv.get(k));
            }
            KnotVector kv = KnotVector.newInstance(degree, Float64Vector.valueOf(newKV));

            //  Pj thru Pi will be overwritten
            j = fout;
            i = j;
            for (int k = 1; k < t; ++k) {
                if ((k % 2) == 1)   //  k modulo 2
                    ++i;
                else
                    --j;
            }

            //  Construct a new control point list.
            FastTable<ControlPoint> newCPs = FastTable.newInstance();
            for (int k = 0; k < j; ++k) {
                newCPs.add(Pw.get(k));
            }
            for (int k = i + 1; k < n; ++k) {
                newCPs.add(Pw.get(k));
            }

            //  Construct a new curve.
            BasicNurbsCurve crv = BasicNurbsCurve.newInstance(newCPs, kv);
            return StackContext.outerCopy(crv);

        } finally {
            StackContext.exit();
        }

    }

    /**
     * Return the total arc length of this curve.
     *
     * @param eps The desired fractional accuracy on the arc-length.
     * @return The total arc length of the curve.
     */
    @Override
    public Parameter<Length> getArcLength(double eps) {

        if (getDegree() == 1)
            //  A polygon is a special case that can be dealt with easily.
            return getCPPolyLength();

        //  Otherwise, just use the standard approach.
        return getArcLength(0, 1, eps);
    }

    /**
     * Return the length of the control point polygon for this curve.
     */
    private Parameter<Length> getCPPolyLength() {

        StackContext.enter();
        try {
            Parameter<Length> length = Parameter.ZERO_LENGTH.to(getUnit());
            List<ControlPoint> Pw = getControlPoints();

            int size = Pw.size();
            Point Pim1 = Pw.get(0).getPoint();
            for (int i = 1; i < size; ++i) {
                Point Pi = Pw.get(i).getPoint();
                length = length.plus(Pim1.distance(Pi));
                Pim1 = Pi;
            }

            return StackContext.outerCopy(length);

        } finally {
            StackContext.exit();
        }
    }

    /**
     * Return the arc length of a segment of this curve.
     *
     * @param s1  The starting parametric distance along the curve to begin the arc-length
     *            calculation from.
     * @param s2  The ending parametric distance along the curve to end the arc-length
     *            calculation at.
     * @param eps The desired fractional accuracy on the arc-length.
     * @return The arc length of the specified segment of the curve.
     */
    @Override
    public Parameter<Length> getArcLength(double s1, double s2, double eps) {
        validateParameter(s1);
        validateParameter(s2);

        //  Method:  Find arc-lengths of each Bezier segment and add them up.
        //           This avoids problems with discontinuities and improves accuracy.
        //           Adjusts eps value for each segment assuming "s" is proportional to arc-length
        //           which is not exactly true (also applied a factor of 2 conservatism).
        
        //  Create a list of Bezier segment starting parametric positions ( and the last end position, 1).
        FastTable<Double> bParams = getSegmentParameters();

        //  Find the arc-length of each segment but the last in turn and sum them up.
        Parameter<Length> L = Parameter.ZERO_LENGTH.to(getUnit());
        int nSegs = bParams.size() - 1;
        double epsf = 2 * (s2 - s1);
        double sprev = s1;
        for (int i = 1; i < nSegs; ++i) {
            double si = bParams.get(i);
            if (si > s1 && si < s2) {
                L = L.plus(super.getArcLength(sprev, si, eps / (si - sprev) * 0.5));
                sprev = si;
                epsf = 1;
            }
        }

        //  Add in the last segment.
        L = L.plus(super.getArcLength(sprev, s2, eps * epsf / (s2 - sprev) * 0.5));

        return L;
    }

    /**
     * Return a string of points that are gridded onto the curve with the number of points
     * and spacing chosen to result in straight line segments between the points that have
     * mid-points that are all within the specified tolerance of this curve.
     *
     * @param tol The maximum distance that a straight line between gridded points may
     *            deviate from this curve. May not be null.
     * @return A string of subrange points gridded onto the curve.
     */
    @Override
    public PointString<SubrangePoint> gridToTolerance(Parameter<Length> tol) {
        if (getDegree() > 1)
            return super.gridToTolerance(tol);
        
        //  A more efficient algorithm is possible for degree=1 (linear) curves.
        //  Since a linear curve consists of straight line segments between knots,
        //  we can start with the knots and then "thin to tolerance" from there.
        
        if (isDegenerate(requireNonNull(tol))) {
            //  If this curve is a single point, just return that point twice.
            PointString<SubrangePoint> str = PointString.newInstance();
            str.add(getPoint(0));
            str.add(getPoint(1));
            return str;
        }
        
        //  Create a list of Bezier segment starting parametric positions ( and the last end position, 1).
        FastTable<Double> bParams = getSegmentParameters();
        
        //  Convert the parametric positions into SubrangePoints.
        PointString<SubrangePoint> pntList = PointString.newInstance();
        int numPnts = bParams.size();
        for (int i = 0; i < numPnts; ++i) {
            Double s = bParams.get(i);
            pntList.add(getPoint(s));
        }
        
        //  Loop over the points and remove any that are closer than tol to line segments
        //  between adjacent points.
        Parameter<Area> tol2 = tol.times(tol);
        for (int i = numPnts - 2; i >= 1; --i) {
            SubrangePoint pi = pntList.get(i);
            double s = pi.getParPosition().getValue(0);

            //  Calculate a point, pa, along the line between p(i+1) and p(i-1)
            //  at the position of "s" (or pi).
            SubrangePoint pip1 = pntList.get(i + 1);
            SubrangePoint pim1 = pntList.get(i - 1);
            Point pa = interpSPoint(pim1, pip1, s);

            //  Compute the distance between the point at "s" and the line seg point
            //  and compare with the tolerance.
            boolean distLTtol = pa.distanceSq(pi).isLessThan(tol2);
            if (distLTtol) {
                //  Remove this point from the list.
                SubrangePoint pnt = pntList.remove(i);
                SubrangePoint.recycle(pnt);
            }
        }
        
        return pntList;
    }
    
    /**
     * Guess an initial set of points to use when gridding to tolerance. The initial
     * points must include the start and end of the curve as well as any discontinuities
     * in the curve. Other than that, they should be distributed as best as possible to
     * indicate areas of higher curvature, but should be generated as efficiently as
     * possible. The guessed points MUST be monotonically increasing from 0 to 1.
     *
     * @param tol2 The square of the maximum distance that a straight line
     *             between gridded points may deviate from this curve.
     * @return A FastList of subrange points to use as a starting point for the grid to
     *         tolerance algorithms.
     */
    @Override
    protected FastList<SubrangePoint> guessGrid2TolPoints(Parameter<Area> tol2) {
        if (getPhyDimension() < 2)
            return super.guessGrid2TolPoints(tol2);

        //  Extract the boundaries of Bezier segments and place grid points at those.
        //  Find the control point that is furthest from a straight line from the start to
        //  the end of each Bezier segment and place a grid point that that approximate location.
        FastTable<Double> parS = FastTable.newInstance();
        StackContext.enter();
        try {
            //  Create a list of Bezier segment starting parametric positions ( and the last end position, 1).
            FastTable<Double> bParams = getSegmentParameters();

            //  Create a list of Bezier segments for the curve.
            GeomList<BasicNurbsCurve> segs = CurveUtils.decomposeToBezier(this);

            //  Loop over the segment starting positions.
            int nSegs = segs.size();
            for (int i = 0; i < nSegs; ++i) {
                //  Always add the start point for each segment to the output points.
                Double s0 = bParams.get(i);
                parS.add(s0);

                //  Get this Bezier segment.
                BasicNurbsCurve seg = segs.get(i);

                //  Get this segment's control points.
                FastTable<ControlPoint> cps = (FastTable<ControlPoint>)seg.getControlPoints();

                //  Find the control point that is furthest away from a line segment between the 1st and last CP.
                Point p0 = cps.get(0).getPoint();
                Point p1 = cps.getLast().getPoint();
                Unit<Length> unit = getUnit();
                double Lv = p1.distance(p0).getValue(unit);
                if (!MathTools.isApproxZero(Lv)) {
                    Vector<Dimensionless> u = p1.minus(p0).toGeomVector().toUnitVector();
                    double dmax2 = 0;
                    double sdmax = 0.5;
                    int numCP = cps.size() - 1;
                    for (int cpIdx = 1; cpIdx < numCP; ++cpIdx) {
                        Point cpi = cps.get(cpIdx).getPoint();
                        Point pi = GeomUtil.pointLineClosest(cpi, p0, u);
                        double dist2 = cpi.distanceSqValue(pi);
                        if (dist2 > dmax2) {
                            //  Store maximum distance (squared).
                            dmax2 = dist2;
                            //  Approx. parametric location of maximum point.
                            sdmax = pi.distance(p0).getValue(unit) / Lv;
                        }
                    }

                    //  Add the location of the maximum lateral offset to the starting points.
                    if (dmax2 > tol2.getValue(unit.pow(2).asType(Area.class))) {
                        //  Convert from segment to full curve parametric position.
                        Double s1 = bParams.get(i + 1);
                        sdmax = segmentPos2Parent(sdmax, s0, s1);
                        parS.add(sdmax);
                    }
                }
            }

            //  Add the last point to the list.
            parS.add(bParams.getLast());

        } finally {
            StackContext.exit();
        }

        //  Convert the parametric positions into SubrangePoints.
        FastList<SubrangePoint> pntList = FastList.newInstance();
        int nSegs = parS.size();
        for (int i = 0; i < nSegs; ++i) {
            Double s = parS.get(i);
            pntList.add(getPoint(s));
        }
        FastTable.recycle(parS);

        return pntList;
    }

    /**
     * Return a list of brackets that surround possible closest/farthest points. Each of
     * the brackets will be refined using a root solver in order to find the actual
     * closest or farthest points.
     *
     * @param point   The point to find the closest/farthest point on this curve to/from.
     * @param closest Set to <code>true</code> to return the closest point or
     *                <code>false</code> to return the farthest point.
     * @param eval    A function that calculates the dot product of (p(s)-q) and ps(s) on
     *                the curve.
     * @return A list of brackets that each surround a potential closest/farthest point.
     */
    @Override
    protected List<BracketRoot1D> guessClosestFarthest(GeomPoint point, boolean closest,
            Evaluatable1D eval) {
        requireNonNull(point);
        requireNonNull(eval);
        
        //  Inspired by: Ma, Hewitt, "Point inversion and projection for NURBS curve and surface", CAGD, 2003.
        //  But not quite as complex (or efficient or robust).
        //  Extract the curve degree.
        int degree = getDegree();
        int order = degree + 1;

        //  Create a list of Bezier segment starting parametric positions ( and the last end position, 1).
        FastTable<Double> bParams = getSegmentParameters();

        //  Decompose the curve into Bezier segments.
        GeomList<BasicNurbsCurve> segs = CurveUtils.decomposeToBezier(this);

        //  Create a list of Bezier end parameters that may be closest/farthest points.
        FastTable<Double> closestEnds = FastTable.newInstance();

        //  Create a list of segments to sub-divide into more segments.
        FastTable<Integer> segIndexes = FastTable.newInstance();

        //  Loop over all the Bezier segments.
        int size = segs.size();
        for (int i = 0; i < size; ++i) {
            BasicNurbsCurve seg = segs.get(i);
            List<ControlPoint> segCP = seg.getControlPoints();

            //  Determine if this segment is a "valid" segment (simple and convex).
            if (CurveUtils.isSimpleConvexPolygon(segCP)) {
                //  Determine if the end-points are closest to the target point.
                if (CurveUtils.isBezierEndPointNearest(point, segCP)) {
                    //  One of the ends is closest.
                    double s = bParams.get(i);  //  Start of segment "i".
                    closestEnds.add(s);
                    s = bParams.get(i + 1);     //  End of segment "i".
                    closestEnds.add(s);

                } else {
                    //  An interior point is likely closest.
                    segIndexes.add(i);
                }

            } else {
                //  Handle invalid segments by just searching for interior brackets.
                segIndexes.add(i);
            }

            FastTable.recycle((FastTable)segCP);
        }
        GeomList.recycle(segs);

        //  Look for brackets in each of the Bezier segments identified.
        FastTable<BracketRoot1D> brackets = findInteriorBrackets(order + 1, segIndexes, bParams, eval);

        //  Finally, look for brackets near the end points of other segments.
        brackets.addAll(findEndBrackets(closestEnds, point, closest, eval));

        //  Look for and remove overlapping brackets that are bracketing the same root.
        removeOverlappingBrackets(brackets, eval);

        //  Cleanup before leaving.
        FastTable.recycle(closestEnds);
        FastTable.recycle(segIndexes);
        FastTable.recycle(bParams);

        return brackets;
    }

    /**
     * Return a list of brackets that surround possible plane intersection points. Each of
     * the brackets will be refined using a root solver in order to find the actual
     * intersection points.
     *
     * @param plane The plane to find the intersection points on this curve to/from.
     * @param eval  A function that calculates the <code>n dot (p(s) - pp)</code> on the
     *              curve.
     * @return A list of brackets that each surround a potential intersection point.
     *         Returns an empty list if there are no intersections found.
     */
    @Override
    protected List<BracketRoot1D> guessIntersect(GeomPlane plane, Evaluatable1D eval) {
        requireNonNull(plane);
        requireNonNull(eval);
        
        //  Get the Bezier segment start parameters.
        FastTable<Double> bParams = getSegmentParameters();

        //  Extract the curve degree.
        int degree = getDegree();
        int numPoints = degree + 3;

        //  Find multiple intersection points (roots) by sub-dividing each Bezier segment into pieces and
        //  bracketing any roots that occur in any of those segments.
        FastTable<BracketRoot1D> outputs = FastTable.newInstance();
        try {
            double s0 = bParams.get(0);
            double fp = eval.function(s0);
            int numSegs = bParams.size() - 1;
            for (int j = 0; j < numSegs; ++j) {
                //  Get the next segment end.
                double s1 = bParams.get(j + 1);

                //  Search for brackets between s0 and s1.
                double dx = (s1 - s0) / numPoints;
                if (dx != 0) {
                    double x = s0;
                    double xp = x;
                    for (int i = 0; i < numPoints; ++i) {
                        x += dx;
                        if (x > s1)
                            x = s1;
                        double fc = eval.function(x);
                        if (fc * fp <= 0.0) {
                            //  Create a bracket (while dealing with roundoff)
                            double xpt = xp - MathTools.epsilon(xp);
                            if (xpt < 0.0)
                                xpt = 0.0;
                            double xt = x + MathTools.epsilon(x);
                            if (xt > 1.0)
                                xt = 1.0;
                            outputs.add(new BracketRoot1D(xpt, xt));
                        }
                        xp = x;
                        fp = fc;
                    }
                }

                //  Set the next segment start to this segment end.
                s0 = s1;
            }
        } catch (RootException err) {
            Logger.getLogger(NurbsCurve.class.getName()).log(Level.WARNING,
                    "Failed to bracket plane intersection points.", err);
        }

        return outputs;
    }

    /**
     * Return a list of brackets that surround possible surface intersection points. Each
     * of the brackets will be refined using a root solver in order to find the actual
     * intersection points.
     *
     * @param surface The surface to find the intersection points on this curve to/from.
     * @param eval    A function that calculates the <code>n dot (p(s) - pp)</code> on the
     *                curve.
     * @return A list of brackets that each surround a potential intersection point.
     *         Returns an empty list if there are no intersections found.
     */
    @Override
    protected List<BracketRoot1D> guessIntersect(Surface surface, Evaluatable1D eval) {

        //  Get the Bezier segment start parameters.
        FastTable<Double> bParams = getSegmentParameters();
        int numSegs = bParams.size() - 1;

        //  Determine the number of points per segment.
        int numPoints = getNumPointsPerSegment();

        //  Don't allow an excessive number of bracket tests.
        if (numPoints * numSegs > 50) {
            numPoints = 50 / numSegs;
            if (numPoints == 0)
                numPoints = 1;
        }

        //  Find multiple intersection points (roots) by sub-dividing curve into pieces and
        //  bracketing any roots that occur in any of those segments.
        FastTable<BracketRoot1D> brackets = FastTable.newInstance();
        try {
            double fp = eval.function(bParams.get(0));
            for (int i = 0; i < numSegs; ++i) {
                double s0 = bParams.get(i);
                double s1 = bParams.get(i + 1);
                fp = findNextBrackets(eval, s0, s1, numPoints, fp, brackets);
            }

        } catch (RootException err) {
            Logger.getLogger(NurbsCurve.class.getName()).log(Level.WARNING,
                    "Failed to bracket surface intersection points.", err);
        }

        return brackets;
    }

    /**
     * Given a function, <code>func</code>, defined on the interval <code>x1</code> to
     * <code>x2</code>, this routine subdivides the interval into <code>n</code> equally
     * spaced segments, and searches for zero crossings of the function. Brackets around
     * any zero crossings found are placed in the input "brackets" list. The function
     * value at x1 is also input to eliminate unnecessary function evaluations when
     * calling this function in a loop.
     *
     * @param func     The function that is being search for zero crossings.
     * @param x1       The start of the interval to be searched.
     * @param x2       The end of the interval to be searched.
     * @param n        The number segments to divide the interval into.
     * @param fx1      The function value at x1. This is used to eliminate a redundant
     *                 function evaluation when this function is called repeatedly in a
     *                 loop.
     * @param brackets An existing list that has any new brackets added to it.
     * @return The value of the function evaluated at x2.
     * @throws RootException if the Evaluatable1D throws a exception.
     */
    private static double findNextBrackets(Evaluatable1D func, double x1, double x2, int n,
            double fx1, List<BracketRoot1D> brackets) throws RootException {

        double dx = (x2 - x1) / n;
        if (MathTools.isApproxZero(dx))
            return fx1;

        double x = x1;
        double fp = fx1;
        for (int i = 0; i < n; ++i) {
            x += dx;
            if (x > x2)
                x = x2;
            double fc = func.function(x);
            if (fc * fp <= 0.0) {
                brackets.add(new BracketRoot1D(x - dx, x));
            }
            fp = fc;
        }

        return fp;
    }

    /**
     * Bracket any closest/farthest points interior to the listed segments by
     * sub-dividing.
     *
     * @param numSubSegs The number of segments to divide the Bezier segment into for root
     *                   bracketing.
     * @param segIndexes The indexes for the Bezier segments to be searched for
     *                   closest/farthest points.
     * @param knotList   A list containing the parametric positions of the start of each
     *                   Bezier segment.
     * @param eval       A function that calculates the dot product of (p(s)-q) and ps(s)
     *                   on the curve.
     * @return A list of bracketed closest/farthest points in the segments indicated by
     *         segIndexes.
     */
    private static FastTable<BracketRoot1D> findInteriorBrackets(int numSubSegs, List<Integer> segIndexes, List<Double> knotList, Evaluatable1D eval) {

        FastTable<BracketRoot1D> brackets = FastTable.newInstance();
        for (int i : segIndexes) {
            double s1 = knotList.get(i);        //  S for start of segment "i".
            double s2 = knotList.get(i + 1);    //  S for end of segment "i".
            try {
                brackets.addAll(BracketRoot1D.findBrackets(eval, s1, s2, numSubSegs));
            } catch (RootException err) {
                Logger.getLogger(NurbsCurve.class.getName()).log(Level.WARNING,
                        "Failed to find interior bracket points.", err);
            }
        }

        return brackets;
    }

    /**
     * Bracket any closest/farthest points at the end of the closest/farthest Bezier
     * segment end.
     *
     * @param segEnds A list of ends of Bezier segments to search for the closest/farthest
     *                point.
     * @param point   The point to find the closest/farthest point on this curve to/from.
     * @param closest Set to <code>true</code> to return the closest point or
     *                <code>false</code> to return the farthest point.
     * @param eval    A function that calculates the dot product of (p(s)-q) and ps(s) on
     *                the curve.
     * @return A list of brackets surrounding the closest/farthest end point of the
     *         supplied Bezier segment ends.
     */
    private List<BracketRoot1D> findEndBrackets(List<Double> segEnds, GeomPoint point, boolean closest,
            Evaluatable1D eval) {

        int size = segEnds.size();
        if (size > 0) {
            //  Find the closest/farthest segment end-point.  The other's don't matter.
            double sOpt = segEnds.get(0);
            StackContext.enter();
            try {
                double dOpt = getRealPoint(sOpt).distanceValue(point);
                double oldS = -1;
                for (int i = 1; i < size; ++i) {
                    double s = segEnds.get(i);
                    if (s != oldS) {
                        double dist = getRealPoint(s).distanceValue(point);
                        if (closest) {
                            if (dist < dOpt) {
                                dOpt = dist;
                                sOpt = s;
                            }
                        } else {
                            if (dist > dOpt) {
                                dOpt = dist;
                                sOpt = s;
                            }
                        }
                    }
                    oldS = s;
                }
            } finally {
                StackContext.exit();
            }

            //  Bracket around the closest end-point.
            try {
                return bracketNear(eval, sOpt);
            } catch (RootException err) {
                Logger.getLogger(NurbsCurve.class.getName()).log(Level.WARNING,
                        "Failed to find end bracket points.", err);
            }
        }

        //  Return an empty list (must be ArrayList as that is what "Roots.bracket" returns).
        return new ArrayList();
    }

    /**
     * Remove any overlapping brackets that are covering the same closest/farthest point.
     *
     * @param brackets A list of root brackets. The list will be modified to delete any
     *                 brackets that are both bracketing the same root.
     * @param eval     A function that calculates the dot product of (p(s)-q) and ps(s) on
     *                 the curve.
     */
    private void removeOverlappingBrackets(List<BracketRoot1D> brackets, Evaluatable1D eval) {

        //  Look for and remove overlapping brackets that are bracketing the same root.
        int size = brackets.size();
        if (size > 1) {
            try {
                //  First sort the brackets into ascending order of the start of each bracket.
                Collections.sort(brackets);

                --size;
                for (int i = 0; i < size; ++i) {
                    BracketRoot1D b = brackets.get(i);
                    BracketRoot1D bp1 = brackets.get(i + 1);
                    double x1 = b.x2;
                    double x2 = bp1.x1;
                    if (x1 > x2) {
                        //  We have an overlapping bracket pair.
                        //  Do they bound the same root?
                        double f1 = eval.function(x1);
                        double f2 = eval.function(x2);
                        if (f1 * f2 <= 0) {
                            //  Both bracket the same root.

                            //  Modify the 1st bracket to more closely cover the root.
                            b.x1 = x2;
                            b.x2 = x1;

                            //  Delete the 2nd bracket.
                            brackets.remove(i + 1);
                            --size;
                        }
                    }
                }
            } catch (RootException err) {
                Logger.getLogger(NurbsCurve.class.getName()).log(Level.WARNING,
                        "Evaluation function failed while removing overlapping brackets.", err);
            }
        }
    }

    /**
     * Returns transformed version of this element. The returned object implements
     * {@link geomss.geom.GeomTransform} and contains this element as a child.
     *
     * @param transform The transformation to apply to this geometry. May not be null.
     * @return A new triangle that is identical to this one with the specified
     *         transformation applied.
     * @throws DimensionException if this point is not 3D.
     */
    @Override
    public NurbsCurveTrans getTransformed(GTransform transform) {
        return NurbsCurveTrans.newInstance(this, requireNonNull(transform));
    }

    /**
     * Split this {@link NurbsCurve} at the specified parametric position returning a list
     * containing two curves (a lower curve with smaller parametric positions than "s" and
     * an upper curve with larger parametric positions).
     *
     * @param s The parametric position where this curve should be split (must not be 0 or
     *          1!).
     * @return A list containing two curves: 0 == the lower curve, 1 == the upper curve.
     */
    @Override
    public GeomList<NurbsCurve> splitAt(double s) {
        validateParameter(s);
        if (parNearEnds(s, TOL_S))
            throw new IllegalArgumentException(
                    MessageFormat.format(RESOURCES.getString("canNotSplitAtEnds"), "curve"));

        GeomList<NurbsCurve> output = GeomList.newInstance();   //  Outside of StackContext.
        StackContext.enter();
        try {
            KnotVector U = getKnotVector();
            int degree = U.getDegree();
            int span = U.findSpan(s);

            //  Find the multiplicity of the split parameter.
            int sIndex;
            if (MathLib.abs(s - U.getValue(span)) <= TOL_S)
                sIndex = U.findMultiplicity(span);
            else
                sIndex = 0;

            //  Number of times to insert the knot required to split the curve.
            int t = degree + 1 - sIndex;

            //  Insert the knot where we are splitting the curve "t" times.
            BasicNurbsCurve cTmp = insertKnot(s, t);
            KnotVector kvTmp = cTmp.getKnotVector();
            List<ControlPoint> cpTmp = cTmp.getControlPoints();

            span = kvTmp.findSpan(s) - degree;      //  span is now the beginning of the upper curve.

            //  Create a knot vector for the lower curve.
            FastTable<Float64> kvList = FastTable.newInstance();
            int size = span + degree + 1;
            for (int i = 0; i < size; ++i) {
                double snew = kvTmp.getValue(i) / s;            //  Scale to range 0-1.
                //  Watch for roundoff
                if (snew <= TOL_S)
                    snew = 0;
                else if (snew >= 1 - TOL_S)
                    snew = 1;
                kvList.add(Float64.valueOf(snew));
            }
            KnotVector kv = KnotVector.newInstance(degree, Float64Vector.valueOf(kvList));

            //  Create a ControlPoint list for the lower curve.
            FastTable<ControlPoint> cpList = FastTable.newInstance();
            size = span;
            for (int i = 0; i < size; ++i) {
                cpList.add(cpTmp.get(i));
            }

            //  Create the lower curve.
            BasicNurbsCurve cl = BasicNurbsCurve.newInstance(cpList, kv);

            //  Create a knot vector for the upper curve.
            double oms = 1 - s;
            kvList.clear();
            size = kvTmp.length();
            for (int i = span; i < size; ++i) {
                double snew = (kvTmp.getValue(i) - s) / oms;    //  Scale to range 0-1.
                //  Watch for roundoff.
                if (snew <= TOL_S)
                    snew = 0;
                else if (snew >= 1 - TOL_S)
                    snew = 1;
                kvList.add(Float64.valueOf(snew));
            }
            kv = KnotVector.newInstance(degree, Float64Vector.valueOf(kvList));

            //  Create a ControlPoint list for the upper curve.
            cpList.clear();
            size = cpTmp.size();
            for (int i = span; i < size; ++i) {
                cpList.add(cpTmp.get(i));
            }

            //  Create the upper curve.
            BasicNurbsCurve cu = BasicNurbsCurve.newInstance(cpList, kv);

            //  Put in the output list (which is outside the StackContext).
            output.add(StackContext.outerCopy(cl));
            output.add(StackContext.outerCopy(cu));

        } finally {
            StackContext.exit();
        }

        return output;
    }

    /**
     * Return the intersection between a line segment and this curve.
     *
     * @param line The line segment to intersect. May not be null.
     * @param tol  Tolerance (in physical space) to refine the point positions to. May not be null.
     * @return A list containing two PointString instances each containing zero or more
     *         subrange points, on this curve and the input line segment respectively,
     *         made by the intersection of this curve with the specified line segment. If
     *         no intersections are found a list of two empty PointStrings are returned.
     */
    @Override
    public GeomList<PointString<SubrangePoint>> intersect(LineSegment line, Parameter<Length> tol) {
        requireNonNull(line);
        requireNonNull(tol);
        
        //  First check to see if the line segment even comes near the curve.
        GeomList<PointString<SubrangePoint>> output = GeomList.newInstance();
        output.add(PointString.newInstance());
        output.add(PointString.newInstance());
        if (!GeomUtil.lineSegBoundsIntersect(line.getStart(), line.getEnd(), this))
            return output;

        //  Check for a degenerate curve case.
        if (isDegenerate(tol)) {
            //  Just see if the start point on this curve intersects that curve (is on that curve).
            SubrangePoint p1 = getPoint(0);
            Point p1r = p1.copyToReal();

            //  Find the closest point on the input curve.
            SubrangePoint p2 = line.getClosest(p1r, GTOL);

            if (p1r.distance(p2).isLessThan(tol)) {
                //  Add the points to the lists of intersections being collected.
                output.get(0).add(p1);
                output.get(1).add(p2);
            }

            return output;
        }

        //  Create a line segment-curve intersector object.
        LineSegNurbsCurveIntersect intersector = LineSegNurbsCurveIntersect.newInstance(this, tol, line, output);

        //  Intersect the line segment and curve (intersections are added to "output").
        output = intersector.intersect();

        //  Clean up.
        LineSegNurbsCurveIntersect.recycle(intersector);

        return output;
    }

    /**
     * Return a NURBS curve representation of this curve to within the specified
     * tolerance. This implementation returns a reference to this curve itself since it is
     * already a NURBS curve and the tolerance parameter is ignored.
     *
     * @param tol Ignored. May pass <code>null</code>.
     * @return A reference to this NURBS curve is returned.
     */
    @Override
    public NurbsCurve toNurbs(Parameter<Length> tol) {
        return this;
    }

    /**
     * Return a NURBS curve representation of this curve. This implementation returns a
     * reference to this curve itself since it is already a NURBS curve.
     *
     * @return A reference to this NURBS curve is returned.
     */
    public NurbsCurve toNurbs() {
        return this;
    }

    /**
     * Return <code>true</code> if this NurbsCurve contains valid and finite numerical
     * components. A value of <code>false</code> will be returned if any of the control
     * point values or knots are NaN or Inf.
     *
     * @return true if this curve contains valid and finite values.
     */
    @Override
    public boolean isValid() {

        //  Check the knot vector.
        if (!getKnotVector().isValid())
            return false;

        //  Check the control points.
        List<ControlPoint> controlPoints = getControlPoints();
        int size = controlPoints.size();
        for (int i = 0; i < size; ++i) {
            ControlPoint cp = controlPoints.get(i);
            if (!cp.isValid())
                return false;
        }

        return true;
    }

    /**
     * Return <code>true</code> if this curve is degenerate (i.e.: has length less than
     * the specified tolerance).
     *
     * @param tol The tolerance for determining if this curve is degenerate. May not be null.
     * @return true if this curve is degenerate.
     * @see #getArcLength
     */
    @Override
    public boolean isDegenerate(Parameter<Length> tol) {
        requireNonNull(tol);
        
        //  For a NURBS curve this is much faster than computing
        //  the length of the curve using getArcLength() and then comparing
        //  with tol.
        Parameter<Length> cLength = getCPPolyLength();
        return cLength.compareTo(tol) <= 0;
    }

    /**
     * Returns <code>true</code> if this curve is a line to within the specified
     * tolerance.
     *
     * @param tol The tolerance for determining if this curve is a line. If null is
     *            passed, then exact linearity is required.
     * @return true if this curve is a line.
     */
    @Override
    public boolean isLine(Parameter<Length> tol) {

        //  Make some simple checks on linearity first.
        if (getDegree() == 1) {
            StackContext.enter();
            try {

                //  Get the control points.
                List<ControlPoint> controlPoints = getControlPoints();

                //  First check for an obvious case.
                int numCp = controlPoints.size();
                if (numCp <= 2) {
                    return true;
                }

                //  Are all the weights equal to 1?
                boolean isLinear = true;
                int size = controlPoints.size();
                for (int i = 0; i < size; ++i) {
                    ControlPoint cp = controlPoints.get(i);
                    if (MathTools.isApproxEqual(cp.getWeight(), 1)) {
                        isLinear = false;
                        break;
                    }
                }

                //  If degree == 1, weights == 1, and all the control points are co-linear,
                //  then the curve is a line.
                if (isLinear) {
                    Point cp0 = controlPoints.get(0).getPoint();
                    Point cp1 = controlPoints.get(1).getPoint();
                    Vector<Dimensionless> u = cp1.minus(cp0).toGeomVector().toUnitVector();
                    for (int i = 2; i < numCp; ++i) {
                        Point cp = controlPoints.get(i).getPoint();
                        Parameter<Length> d = GeomUtil.pointLineDistance(cp, cp0, u);
                        if (d.isGreaterThan(tol)) {
                            isLinear = false;
                            break;
                        }
                    }
                    if (isLinear)
                        return true;
                }

            } finally {
                StackContext.exit();
            }
        }

        //  If simple checks fail, try more complex tests for linearity.
        return super.isLine(tol);
    }

    /**
     * Return <code>true</code> if this curve is planar or <code>false</code> if it is
     * not.
     *
     * @param tol The geometric tolerance to use in determining if the curve is planar.
     *            May not be null.
     * @return true if this curve is planar.
     */
    @Override
    public boolean isPlanar(Parameter<Length> tol) {
        //  If the curve is less than 3D, then it must be planar.
        int numDims = getPhyDimension();
        if (numDims <= 2)
            return true;

        //  Is this curve degenerate?
        if (isDegenerate(requireNonNull(tol)))
            return true;

        //  Is the curve linear?
        if (isLine(tol))
            return true;

        //  The curve is planar if the control points are co-planar.
        StackContext.enter();
        try {

            //  Get the control points.
            List<ControlPoint> cps = getControlPoints();

            //  If there are only two or three control points, then the curve is planar by definition
            //  (it is a line or triangle).
            int numCP = cps.size();
            if (numCP <= 3)
                return true;

            //  Form a vector with the 1st two control points.
            Point p0 = cps.get(0).getPoint();
            Point p1 = cps.get(1).getPoint();
            Vector<Length> v1 = p1.toGeomVector().minus(p0.toGeomVector());

            //  Form a vector with the 2nd and 3rd control points.
            Point p2 = cps.get(2).getPoint();
            Vector<Length> v2 = p2.toGeomVector().minus(p1.toGeomVector());

            //  Cross product of these two vectors defines the plane.
            Point ZERO_POINT = Point.newInstance(numDims);
            Vector<Dimensionless> nhat = v1.cross(v2).toUnitVector();
            nhat.setOrigin(ZERO_POINT);

            //  Find deviation from a plane for each point by projecting a vector from p1 to pi
            //  onto the normal formed by the cross product of the 1st two vectors.
            //  If the projection is anything other than zero, the points and tangent
            //  vectors are not in the same plane.
            for (int i = 3; i < numCP; ++i) {
                Point pi = cps.get(i).getPoint();
                Vector<Length> p1pi = pi.minus(p1).toGeomVector();
                if (!p1pi.mag().isApproxZero()) {
                    Parameter proj = p1pi.dot(nhat);
                    if (proj.isLargerThan(tol))
                        return false;
                }
            }

        } finally {
            StackContext.exit();
        }

        return true;
    }

    /**
     * Return a list containing parameters at the start of logical segments of this curve.
     * The first element in this list must always be 0.0 and the last element 1.0. These
     * should be good places to break the curve up into pieces for analysis, but otherwise
     * are arbitrary.
     * <p>
     * This implementation returns the parameters at the start of each Bezier segment of
     * this NURBS curve.
     * </p>
     *
     * @return A list of parameters at the start and end of each Bezier segment.
     */
    @Override
    protected FastTable<Double> getSegmentParameters() {
        return CurveUtils.getBezierStartParameters(this);
    }

    /**
     * Return the number of points per segment that should be used when analyzing curve
     * segments by certain methods.
     * <p>
     * This implementation always returns 2*(p + 1) where p is the degree of the curve.
     * </p>
     *
     * @return The number of points per segment that should be used when analyzing this
     *         curve.
     */
    @Override
    protected int getNumPointsPerSegment() {
        return 2 * (getDegree() + 1);
    }

    /**
     * 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(" = {");
        }
        List<ControlPoint> cpList = getControlPoints();
        int size = cpList.size();
        for (int i = 0; i < size; i++) {
            tmp.append(cpList.get(i).toText());
            if (i != size - 1) {
                tmp.append(", ");
            }
        }
        tmp.append("},");

        tmp.append(getKnotVector().toText());

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

    /**
     * A class used to intersect an line segment and this NURBS curve.
     */
    private static class LineSegNurbsCurveIntersect {

        private NurbsCurve _thisCurve;
        private Parameter<Length> _tol;
        private LineSegment _line;
        private GeomPoint _p0;
        private GeomPoint _p1;
        private Vector<Length> _Lu;
        private GeomList<PointString<SubrangePoint>> _output;

        public static LineSegNurbsCurveIntersect newInstance(NurbsCurve thisCurve, Parameter<Length> tol,
                LineSegment line, GeomList<PointString<SubrangePoint>> output) {

            //  Make sure the line and curve have the same dimensions.
            int numDims = thisCurve.getPhyDimension();
            if (line.getPhyDimension() > numDims)
                throw new IllegalArgumentException(MessageFormat.format(
                        RESOURCES.getString("sameOrFewerDimErr"), "line", "curve"));
            else if (line.getPhyDimension() < numDims)
                line = line.toDimension(numDims);

            if (thisCurve instanceof GeomTransform)
                thisCurve = thisCurve.copyToReal();
            if (line instanceof GeomTransform)
                line = line.copyToReal();

            //  Create an instance of this class and fill in the class variables.
            LineSegNurbsCurveIntersect o = FACTORY.object();
            o._thisCurve = thisCurve;
            o._tol = tol;
            o._line = line;
            o._p0 = line.getStart();
            o._p1 = line.getEnd();
            o._Lu = line.getDerivativeVector().immutable();
            o._output = output;

            return o;
        }

        /**
         * Method that actually carries out the intersection adding the intersection
         * points to the list of lists provided in the constructor.
         *
         * @return A list containing two PointString instances each containing zero or
         *         more subrange points, on this curve and the input line segment
         *         respectively, made by the intersection of this curve with the specified
         *         line segment. If no intersections are found a list of two empty
         *         PointStrings are returned.
         */
        public GeomList<PointString<SubrangePoint>> intersect() {

            //  Decompose the NurbsCurve into a set of Bezier curve segments.
            GeomList<BasicNurbsCurve> segs = CurveUtils.decomposeToBezier(_thisCurve);

            //  Create a list of segment starting par. positions (and the last end position, 1).
            FastTable<Double> bParams = _thisCurve.getSegmentParameters();

            //  Work with each segment separately.
            int numSegs = segs.size();
            for (int i = 0; i < numSegs; ++i) {
                NurbsCurve curveSeg = segs.get(i);
                if (GeomUtil.lineSegBoundsIntersect(_p0, _p1, curveSeg)) {
                    double s0 = bParams.get(i);
                    double s1 = bParams.get(i + 1);
                    divideAndConquerLineSeg(curveSeg, s0, s1);
                }
            }

            //  Clean up.
            FastTable.recycle(bParams);
            GeomList.recycle(segs);

            return _output;
        }

        /**
         * Uses a recursive "Divide and Conquer" approach to intersecting a line segment
         * with a curve. On each call, the following happens:
         * <pre>
         *      1) The curve is checked to see if it is approx. a straight line.
         *          1b) If it is, then a line-line intersection is performed, the
         *              intersection point added to the output list and the method exited.
         *      2) The curve is divided into two halves.
         *          2b) Each half is tested to see if it's bounding box is intersected by the line.
         *              If it is, then this method is called with that segment recursively.
         *              Otherwise, the method is exited.
         * </pre>
         * A number of class variables are used to pass information to this
         * recursion:
         * <pre>
         *      _tol  The tolerance to use when determining if the curve is locally linear.
         *      _line The line segment instance we are intersecting the curve with.
         *      _p0   A point at the start of the line being intersected with this curve.
         *      _p1   A point at the end of the line being intersected with this curve.
         *      _Lu   A vector indicating the direction of the line.
         *      _output A list used to collect the output subrange intersection points.
         * </pre>
         *
         * @param crv The curve segment being checked for intersection with the line.
         * @param s0  The parametric position of the start of this curve segment on the
         *            master curve.
         * @param s1  The parametric position of the end of this curve segment on the
         *            master curve.
         */
        private void divideAndConquerLineSeg(NurbsCurve crv, double s0, double s1) {

            //  Check to see if the input curve is within tolerance of a straight line.
            Point p0 = crv.getRealPoint(0);
            Point p1 = crv.getRealPoint(1);
            Point pAvg = p0.plus(p1).divide(2);
            Point pm = crv.getRealPoint(0.5);
            if (pAvg.distance(pm).isLessThan(_tol)) {
                //  The input curve segment is nearly a straight line.

                //  Intersect the input line and the line approximating the curve.
                Vector<Length> crvDir = p1.minus(p0).toGeomVector();
                MutablePoint sLn = MutablePoint.newInstance(2);
                int numDims = p0.getPhyDimension();
                Unit<Length> unit = p0.getUnit();
                MutablePoint pi1 = MutablePoint.newInstance(numDims, unit);
                MutablePoint pi2 = MutablePoint.newInstance(numDims, unit);
                IntersectType type = GeomUtil.lineLineIntersect(_p0, _Lu, p0, crvDir, _tol, sLn, pi1, pi2);
                Vector.recycle(crvDir);

                //  Make sure the intersection points fall on the line segments.
                double sl_1 = sLn.getValue(0);
                if (sl_1 > 1) {
                    pi1.set(_line.getEnd());
                    sl_1 = 1;
                } else if (sl_1 < 0) {
                    pi1.set(_line.getStart());
                    sl_1 = 0;
                }
                double sl_2 = sLn.getValue(1);
                if (sl_2 > 1) {
                    pi2.set(crv.getRealPoint(1));
                    sl_2 = 1;
                } else if (sl_2 < 0) {
                    pi2.set(crv.getRealPoint(0));
                    sl_2 = 0;
                }
                if (pi1.distance(pi2).isLessThan(_tol)) {
                    //  Add an intersection point on each curve to the output lists.
                    _output.get(1).add(_line.getPoint(sl_1));
                    double s = s0 + (s1 - s0) * sl_2;   //  Convert from segment par. pos. to full curve pos.
                    s = roundParNearEnds(s);            //  Deal with roundoff issues.
                    _output.get(0).add(_thisCurve.getPoint(s));
                }

                //  Clean up.
                MutablePoint.recycle(sLn);
                MutablePoint.recycle(pi1);
                MutablePoint.recycle(pi2);

            } else {
                //  Sub-divide the curve into two segments.
                GeomList<NurbsCurve> segments = crv.splitAt(0.5);

                //  Check for possible intersections on the left segment.
                NurbsCurve seg = segments.get(0);
                if (GeomUtil.lineSegBoundsIntersect(_p0, _p1, seg)) {
                    //  May be an intersection.
                    double sl = s0;
                    double sh = (s0 + s1) / 2;

                    //  Recurse down to a finer level of detail.
                    divideAndConquerLineSeg(seg, sl, sh);
                }
                BasicNurbsCurve.recycle((BasicNurbsCurve)seg);

                //  Check for possible intersections on the right segment.
                seg = segments.get(1);
                if (GeomUtil.lineSegBoundsIntersect(_p0, _p1, seg)) {
                    //  May be an intersection.
                    double sl = (s0 + s1) / 2;
                    double sh = s1;

                    //  Recurse down to a finer level of detail.
                    divideAndConquerLineSeg(seg, sl, sh);
                }
                BasicNurbsCurve.recycle((BasicNurbsCurve)seg);

                //  Clean up.
                GeomList.recycle(segments);
            }

        }

        public static void recycle(LineSegNurbsCurveIntersect instance) {
            FACTORY.recycle(instance);
        }

        private LineSegNurbsCurveIntersect() { }

        private static final ObjectFactory<LineSegNurbsCurveIntersect> FACTORY = new ObjectFactory<LineSegNurbsCurveIntersect>() {
            @Override
            protected LineSegNurbsCurveIntersect create() {
                return new LineSegNurbsCurveIntersect();
            }

            @Override
            protected void cleanup(LineSegNurbsCurveIntersect obj) {
                obj._thisCurve = null;
                obj._tol = null;
                obj._line = null;
                obj._p0 = null;
                obj._p1 = null;
                obj._Lu = null;
                obj._output = null;
            }
        };
    }

    /**
     * Holds the default XML representation for this object.
     */
    @SuppressWarnings("FieldNameHidesFieldInSuperclass")
    protected static final XMLFormat<NurbsCurve> XML = new XMLFormat<NurbsCurve>(NurbsCurve.class) {

        @Override
        public void read(XMLFormat.InputElement xml, NurbsCurve obj) throws XMLStreamException {
            String name = xml.getAttribute("name", (String)null);
            if (nonNull(name))
                obj.setName(name);
            FastMap userData = xml.get("UserData", FastMap.class);
            if (nonNull(userData)) {
                obj.putAllUserData(userData);
            }
        }

        @Override
        public void write(NurbsCurve obj, XMLFormat.OutputElement xml) throws XMLStreamException {
            String name = obj.getName();
            if (nonNull(name))
                xml.setAttribute("name", name);
            FastMap userData = (FastMap)obj.getAllUserData();
            if (!userData.isEmpty())
                xml.add(userData, "UserData", FastMap.class);
        }
    };

}