001/**
002 * CurveUtils -- Utility methods for working with NURBS curves.
003 *
004 * Copyright (C) 2009-2025, by Joseph A. Huwaldt. All rights reserved.
005 *
006 * This library is free software; you can redistribute it and/or modify it under the terms
007 * of the GNU Lesser General Public License as published by the Free Software Foundation;
008 * either version 2 of the License, or (at your option) any later version.
009 *
010 * This library is distributed in the hope that it will be useful, but WITHOUT ANY
011 * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A
012 * PARTICULAR PURPOSE. See the GNU Library General Public License for more details.
013 *
014 * You should have received a copy of the GNU Lesser General Public License along with
015 * this program; if not, write to the Free Software Foundation, Inc., 59 Temple Place -
016 * Suite 330, Boston, MA 02111-1307, USA. Or visit: http://www.gnu.org/licenses/lgpl.html
017 *
018 * This source file is based on, but heavily modified from, the LGPL jgeom (Geometry
019 * Library for Java) code by Samuel Gerber, Copyright (C) 2005.
020 */
021package geomss.geom.nurbs;
022
023import geomss.geom.*;
024import jahuwaldt.js.param.Parameter;
025import jahuwaldt.tools.math.MathTools;
026import java.text.MessageFormat;
027import java.util.Arrays;
028import java.util.Collections;
029import java.util.List;
030import static java.util.Objects.requireNonNull;
031import java.util.ResourceBundle;
032import javax.measure.quantity.Angle;
033import javax.measure.quantity.Dimensionless;
034import javax.measure.quantity.Length;
035import javax.measure.unit.SI;
036import javax.measure.unit.Unit;
037import javolution.context.ArrayFactory;
038import javolution.context.ImmortalContext;
039import javolution.context.StackContext;
040import static javolution.lang.MathLib.max;
041import static javolution.lang.MathLib.min;
042import javolution.util.FastSet;
043import javolution.util.FastTable;
044import org.jscience.mathematics.number.Float64;
045import org.jscience.mathematics.vector.Float64Vector;
046
047/**
048 * A collection of utility methods for working with NURBS curves.
049 *
050 * <p> Modified by: Joseph A. Huwaldt </p>
051 *
052 * @author Samuel Gerber, Date: May 14, 2009, Version 1.0.
053 * @version February 17, 2025
054 */
055public final class CurveUtils {
056
057    /**
058     * The resource bundle for this package.
059     */
060    private static final ResourceBundle RESOURCES = geomss.geom.AbstractGeomElement.RESOURCES;
061
062    private static final Parameter<Angle> ANGULAR_TOLERANCE;
063
064    static {
065        ImmortalContext.enter();
066        try {
067            ANGULAR_TOLERANCE = Parameter.valueOf(0.01, SI.RADIAN);
068
069        } finally {
070            ImmortalContext.exit();
071        }
072    }
073
074    /**
075     * Elevate the degree of a curve to the specified degree. If the specified degree is
076     * &le; the current degree of the curve, then no change is made and a reference to the
077     * input curve is returned.
078     *
079     * @param curve  The curve to have the degree elevated. May not be null.
080     * @param degree The desired degree for the curve to be elevated to.
081     * @return The input NurbsCurve with the degree elevated to the specified degree.
082     */
083    public static NurbsCurve elevateToDegree(NurbsCurve curve, int degree) {
084        int p = curve.getDegree();
085        if (degree <= p)
086            return curve;
087
088        int t = degree - p;
089        return degreeElevate(curve, t);
090    }
091
092    /**
093     * Elevate the degree of a curve by the specified number of times.
094     *
095     * @param curve The curve to have the degree elevated. May not be null.
096     * @param t     The number of times to elevate the degree.
097     * @return The input NurbsCurve with the degree elevated by the specified amount.
098     */
099    public static NurbsCurve degreeElevate(NurbsCurve curve, int t) {
100        //  This is A5.9 from The NURBS Book, pg. 206.
101        requireNonNull(curve);
102        if (t <= 0)
103            return curve;
104
105        StackContext.enter();
106        try {
107            Unit<Length> unit = curve.getUnit();
108            int dim = curve.getPhyDimension();
109            KnotVector U = curve.getKnotVector();
110            List<ControlPoint> Pw = curve.getControlPoints();
111
112            int n = Pw.size() - 1;
113            int p = curve.getDegree();
114            int m = n + p + 1;
115            int ph = p + t;
116            int ph2 = ph / 2;
117
118            // Coefficients for degree elevating the Bezier segment
119            double[][] bezalfs = allocate2DArray(ph + 1, p + 1);
120            // pth-degree Bezier control points of the current segment
121            ControlPoint[] bpts = CONTROLPOINTARRAY_FACTORY.array(p + 1);
122            // (p+t)th-degree Bezier control points of the  current segment
123            ControlPoint[] ebpts = CONTROLPOINTARRAY_FACTORY.array(ph + 1);
124            // leftmost control points of the next Bezier segment
125            ControlPoint[] Nextbpts = CONTROLPOINTARRAY_FACTORY.array(p - 1);
126            // Knot instertion alphas
127            double[] alphas = ArrayFactory.DOUBLES_FACTORY.array(p - 1);
128
129            // Compute the binomial coefficients
130            double[][] Bin = allocate2DArray(ph + 1, ph2 + 1);
131            binomialCoef(Bin, ph + 1, ph2 + 1);
132
133            // Compute Bezier degree elevation coefficients
134            bezalfs[0][0] = bezalfs[ph][p] = 1.0;
135            for (int i = 1; i <= ph2; i++) {
136                double inv = 1.0 / Bin[ph][i];
137                int mpi = min(p, i);
138                for (int j = max(0, i - t); j <= mpi; j++) {
139                    bezalfs[i][j] = inv * Bin[p][j] * Bin[t][i - j];
140                }
141            }
142
143            for (int i = ph2 + 1; i < ph; i++) {
144                int mpi = min(p, i);
145                for (int j = max(0, i - t); j <= mpi; j++) {
146                    bezalfs[i][j] = bezalfs[ph - i][p - j];
147                }
148            }
149
150            // Allocate more control points and knot vector space than necessary
151            ControlPoint[] Qw = CONTROLPOINTARRAY_FACTORY.array(Pw.size() + Pw.size() * t);
152            double[] Uh = ArrayFactory.DOUBLES_FACTORY.array(Pw.size() + Pw.size() * t + ph + 1);
153
154            int mh = ph;
155            int kind = ph + 1;
156            double ua = U.getValue(0);
157            double ub;
158            int r = -1;
159            int oldr;
160            int a = p;
161            int b = p + 1;
162            int cind = 1;
163            int rbz, lbz;
164
165            Qw[0] = Pw.get(0);
166            for (int i = 0; i <= ph; i++) {
167                Uh[i] = ua;
168            }
169
170            // Initialize the first Bezier segment
171            for (int i = 0; i <= p; i++) {
172                bpts[i] = Pw.get(i);
173            }
174
175            while (b < m) { // Big loop thru knot vector
176                int i = b;
177                //  while (b < m && U.getValue(b) == U.getValue(b+1))
178                while (b < m && MathTools.isApproxEqual(U.getValue(b), U.getValue(b + 1))) {
179                    ++b;
180                }
181                int mul = b - i + 1;
182                mh += mul + t;
183                ub = U.getValue(b);
184                oldr = r;
185                r = p - mul;
186                //  Insert knot u(b) r times
187                if (oldr > 0)
188                    lbz = (oldr + 2) / 2;
189                else
190                    lbz = 1;
191                if (r > 0)
192                    rbz = ph - (r + 1) / 2;
193                else
194                    rbz = ph;
195                if (r > 0) {
196                    // Insert knot to get Bezier segment
197                    double numer = ub - ua;
198                    for (int k = p; k > mul; k--) {
199                        alphas[k - mul - 1] = numer / (U.getValue(a + k) - ua);
200                    }
201                    for (int j = 1; j <= r; j++) {
202                        int save = r - j;
203                        int s = mul + j;
204                        for (int k = p; k >= s; k--) {
205                            //  bpts[k] = alphas[k-s]*bpts[k] + (1.0-alphas[k-s])*bpts[k-1];
206                            bpts[k] = applyAlpha(bpts[k], bpts[k - 1], alphas[k - s]);
207                        }
208                        Nextbpts[save] = bpts[p];
209                    }
210                }   //  End of inert knot
211
212                for (i = lbz; i <= ph; i++) { // Degree elevate Bezier,  only the points lbz,...,ph are used
213                    ebpts[i] = ControlPoint.newInstance(dim, unit); //  ebpts[i] = 0.0;
214                    int mpi = min(p, i);
215                    for (int j = max(0, i - t); j <= mpi; j++) {
216                        //  ebpts[i] += bezalfs(i,j)*bpts[j]
217                        ControlPoint q1 = bpts[j].applyWeight();
218                        q1 = q1.times(bezalfs[i][j]);
219                        ControlPoint q2 = ebpts[i].applyWeight();
220                        q2 = q2.plus(q1);
221                        ebpts[i] = q2.getHomogeneous();
222                    }
223                }   //  End of degree elevating Bezier.
224
225                if (oldr > 1) { // Must remove knot u=U[a] oldr times
226                    // if(oldr>2) // Alphas on the right do not change
227                    //  alfj = (ua-U[kind-1])/(ub-U[kind-1]);
228                    int first = kind - 2;
229                    int last = kind;
230                    double den = ub - ua;
231                    double bet = (ub - Uh[kind - 1]) / den;
232                    for (int tr = 1; tr < oldr; tr++) { // Knot removal loop
233                        i = first;
234                        int j = last;
235                        int kj = j - kind + 1;
236                        while (j - i > tr) { // Loop and compute the new control points for one removal step
237                            if (i < cind) {
238                                double alf = (ub - Uh[i]) / (ua - Uh[i]);
239                                //  Qw[i] = alf*Qw[i] + (1.0-alf)*Qw[i-1]
240                                Qw[i] = applyAlpha(Qw[i], Qw[i - 1], alf);
241                            }
242                            if (j >= lbz) {
243                                if (j - tr <= kind - ph + oldr) {
244                                    double gam = (ub - Uh[j - tr]) / den;
245                                    //  ebpts[kj] = gam*ebpts[kj] + (1.0-gam)*ebpts[kj+1];
246                                    ebpts[kj] = applyAlpha(ebpts[kj], ebpts[kj + 1], gam);
247
248                                } else {
249                                    //  ebpts[kj] = bet*ebpts[kj]+(1.0-bet)*ebpts[kj+1];
250                                    ebpts[kj] = applyAlpha(ebpts[kj], ebpts[kj + 1], bet);
251                                }
252                            }
253                            ++i;
254                            --j;
255                            --kj;
256                        }
257                        --first;
258                        ++last;
259                    }
260                }   //  End removing knot u=U[a]
261
262                if (a != p) // load the knot u=U[a]
263                    for (i = 0; i < ph - oldr; i++) {
264                        Uh[kind++] = ua;
265                    }
266                for (int j = lbz; j <= rbz; j++) { // load control points onto the curve
267                    Qw[cind++] = ebpts[j];
268                }
269
270                if (b < m) {
271                    System.arraycopy(Nextbpts, 0, bpts, 0, r);
272                    for (int j = r; j <= p; j++) {
273                        bpts[j] = Pw.get(b - p + j);
274                    }
275                    a = b;
276                    ++b;
277                    ua = ub;
278
279                } else {
280                    for (i = 0; i <= ph; i++) {
281                        Uh[kind + i] = ub;
282                    }
283                }
284            }   //  End of while-loop (b < m)
285
286            //  Create the final lists of control points and knot vector.
287            FastTable<ControlPoint> newP = FastTable.newInstance();
288            n = mh - ph;
289            for (int i = 0; i < n; ++i) {
290                newP.add(Qw[i]);
291            }
292
293            FastTable<Float64> newKVList = FastTable.newInstance();
294            m = n + ph + 1;
295            for (int i = 0; i < m; ++i) {
296                newKVList.add(Float64.valueOf(Uh[i]));
297            }
298            KnotVector newKV = KnotVector.newInstance(ph, Float64Vector.valueOf(newKVList));
299
300            //  Create the new curve.
301            BasicNurbsCurve newCrv = BasicNurbsCurve.newInstance(newP, newKV);
302            return StackContext.outerCopy(newCrv);
303
304        } finally {
305            StackContext.exit();
306        }
307    }
308
309    /**
310     * Decompose a NURBS curve into a list of Bezier segments. The returned list contains
311     * NURBS curves each of which represents a single Bezier segment.
312     *
313     * @param curve The curve to be decomposed into Bezier segments. May not be null.
314     * @return A list of Bezier curve segments (will be at least degree = 2).
315     */
316    public static GeomList<BasicNurbsCurve> decomposeToBezier(NurbsCurve curve) {
317        //  Algorithm A5.6, pg. 173, Ref. 1.
318
319        GeomList<BasicNurbsCurve> output = GeomList.newInstance();    //  Created outside of StackContext.
320        StackContext.enter();
321        try {
322            //  Make sure the curve degree is > 1.
323            if (curve.getDegree() == 1)
324                curve = degreeElevate(curve, 1);
325
326            List<ControlPoint> Pw = curve.getControlPoints();
327            KnotVector SP = curve.getKnotVector();
328            int numKnots = SP.length();
329            int p = SP.getDegree();
330
331            FastTable<FastTable<ControlPoint>> Qw = FastTable.newInstance();
332            FastTable<ControlPoint> Qwnb = FastTable.newInstance();
333            Qw.add(Qwnb);
334            for (int i = 0; i <= p; ++i) {
335                Qwnb.add(Pw.get(i));
336            }
337
338            int m = numKnots - 1;
339            int a = p;
340            int b = p + 1;
341            int nb = 0;
342            while (b < m) {
343                int i = b;
344                // while (b < m && SP[b] == SP[b+1])
345                while (b < m && MathTools.isApproxEqual(SP.getValue(b), SP.getValue(b + 1))) {
346                    ++b;
347                }
348                int mult = b - i + 1;
349                if (mult < p) {
350                    double numer = SP.getValue(b) - SP.getValue(a); //  Numerator of alpha.
351
352                    //  Compute and store alphas.
353                    double[] alphas = ArrayFactory.DOUBLES_FACTORY.array(p - mult);
354                    for (int j = p; j > mult; --j) {
355                        alphas[j - mult - 1] = numer / (SP.getValue(a + j) - SP.getValue(a));
356                    }
357
358                    //  Insert knot r times.
359                    int r = p - mult;
360                    for (int j = 1; j <= r; ++j) {
361                        int save = r - j;
362                        int si = mult + j;  //  This many new points.
363                        for (int k = p; k >= si; --k) {
364                            double alpha = alphas[k - si];
365                            //  Qwnb[k] = alpha*Qwnb[k] + (1.0 - alpha)*Qwnb[k-1];
366                            Qwnb.set(k, applyAlpha(Qwnb.get(k), Qwnb.get(k - 1), alpha));
367                        }
368                        if (b < numKnots) {
369                            //  Control point of next segment.
370                            if (Qw.size() <= nb + 1) {
371                                FastTable<ControlPoint> tmp = FastTable.newInstance();
372                                Qw.add(tmp);
373                            }
374                            FastTable<ControlPoint> Qwnbp1 = Qw.get(nb + 1);
375                            if (Qwnbp1.size() <= save)
376                                Qwnbp1.setSize(save + 1);
377                            Qwnbp1.set(save, Qwnb.get(p));
378                        }
379                    }
380                    ArrayFactory.DOUBLES_FACTORY.recycle(alphas);
381
382                    //  Bezier segment completed.
383                    ++nb;
384                    if (Qw.size() <= nb) {
385                        FastTable tmp = FastTable.newInstance();
386                        tmp.setSize(p + 1);
387                        Qw.add(tmp);
388                    }
389                    Qwnb = Qw.get(nb);
390                    if (b < numKnots) {
391                        //  Initialize for next segment.
392                        if (Qwnb.size() <= p)
393                            Qwnb.setSize(p + 1);
394                        for (i = p - mult; i <= p; ++i) {
395                            Qwnb.set(i, Pw.get(b - p + i));
396                        }
397                        a = b;
398                        ++b;
399                    }
400                }
401            }
402
403            //  Turn the list of lists of control points into a list of curves.
404            //  Build up a knot vector for the new Bezier segments.
405            KnotVector bezierKV = bezierKnotVector(p);
406
407            int size = Qw.size();
408            for (int i = 0; i < size; ++i) {
409                FastTable<ControlPoint> cpLst = Qw.get(i);
410                BasicNurbsCurve crv = BasicNurbsCurve.newInstance(cpLst, bezierKV);
411                output.add(StackContext.outerCopy(crv));        //  "output" is outside of StackContext.
412            }
413
414        } finally {
415            StackContext.exit();
416        }
417
418        return output;
419    }
420
421    /**
422     * Returns a list containing the parameters at the start (and end) of each Bezier
423     * segment of the specified NURBS curve. This first element of this list will always
424     * be zero and the last will always be 1.
425     *
426     * @param crv The NURBS curve to return the Bezier segment start parameters for. May
427     *            not be null.
428     * @return A list containing the parameters of the start and end of each Bezier
429     *         segment.
430     */
431    public static FastTable<Double> getBezierStartParameters(NurbsCurve crv) {
432
433        //  Extract the knot vector and curve degree.
434        KnotVector kv = crv.getKnotVector();
435        int degree = kv.getDegree();
436
437        //  Create a list of interior knots.
438        FastTable<Double> knots = FastTable.newInstance();
439        int size = kv.length() - degree;
440        double oldS = -1;
441        for (int i = degree; i < size; ++i) {
442            double s = kv.getValue(i);
443            if (s != oldS) {
444                knots.add(s);
445                oldS = s;
446            }
447        }
448
449        return knots;
450    }
451
452    /**
453     * Return a knot vector that can be used to create a Bezier curve segment of the
454     * specified degree using the BasicNurbsCurve class.
455     *
456     * @param degree The degree of the knot vector to return.
457     * @return A knot vector for creating a Bezier curve segment of the specified degree.
458     */
459    public static KnotVector bezierKnotVector(int degree) {
460        return KnotVector.bezierKnotVector(degree);
461    }
462
463    /**
464     * Method that applies the specified factor to the two specified breakpoints in the
465     * form: <code>P = alpha*P1 + (1-alpha)*P2</code>.
466     */
467    private static ControlPoint applyAlpha(ControlPoint P1, ControlPoint P2, double alpha) {
468        ControlPoint term1 = P1.times(alpha);
469        ControlPoint term2 = P2.times(1.0 - alpha);
470        ControlPoint P = term1.plus(term2);
471        return P;
472    }
473
474    /**
475     * Connect together or "concatenate" a list of curves end-to-start. The end of each
476     * curve (u=1) is connected to the start (u=0) of the next curve. No check is made for
477     * the curve ends being near each other. The resulting curve will have the physical
478     * dimension of the highest dimension curve in the input list. The resulting curve
479     * will have the degree of the highest degree curve in the input list.
480     *
481     * @param curves The list of curves to be connected together. May not be null.
482     * @return A single {@link BasicNurbsCurve} made up of the input curves connected
483     *         together end-to-start.
484     * @throws IllegalArgumentException if there are less than 2 curves.
485     */
486    public static BasicNurbsCurve connectCurves(List<NurbsCurve> curves) throws IllegalArgumentException {
487        //  Reference:  Yu, T., Soni, B., "Handbook of Grid Generation", Section 30.3.4
488
489        int numCurves = curves.size();
490        if (numCurves < 2) {
491            throw new IllegalArgumentException(RESOURCES.getString("reqTwoCurves"));
492        }
493
494        StackContext.enter();
495        try {
496            //  Determine the highest degree and dimension curve in the list.
497            int degree = 0;
498            int dim = 0;
499            for (int i = 0; i < numCurves; ++i) {
500                NurbsCurve crv = curves.get(i);
501                if (crv.getDegree() > degree) {
502                    degree = crv.getDegree();
503                }
504                if (crv.getPhyDimension() > dim) {
505                    dim = crv.getPhyDimension();
506                }
507            }
508
509            //  Create a knot list and put a degree+1 touple knot at the start.
510            FastTable<Float64> knots = FastTable.newInstance();
511            for (int i = 0; i <= degree; i++) {
512                knots.add(Float64.ZERO);
513            }
514
515            //  Create a list of control points.
516            FastTable<ControlPoint> cps = FastTable.newInstance();
517
518            //  Loop over all the curves in the input list.
519            boolean firstPass = true;
520            double w1 = 1;
521            for (int i = 0; i < numCurves; ++i) {
522                NurbsCurve crv = curves.get(i);
523
524                //  Make sure the curve has the highest physical dimension.
525                crv = crv.toDimension(dim);
526
527                //  Elevate the curve to the highest degree (if required).
528                crv = elevateToDegree(crv, degree);
529
530                //  Get the knot vector for this curve.
531                KnotVector u = crv.getKnotVector();
532                int uLength = u.length();
533
534                //  Add to the knot vector we are building.
535                double start = knots.getLast().doubleValue() - u.getValue(0);
536                int end = uLength - degree - 1;
537                for (int j = degree + 1; j < end; j++) {
538                    knots.addLast(Float64.valueOf(start + u.getValue(j)));
539                }
540
541                //  Add a degree+1 touple knot at the junction of the curves and at the end.
542                Float64 value = Float64.valueOf(start + u.getValue(uLength - 1));
543                for (int j = 0; j < degree; j++) {
544                    knots.addLast(value);
545                }
546
547                //  Create the list of control points we are building.
548                List<ControlPoint> pts = crv.getControlPoints();
549                double wfactor = 1;
550                if (firstPass) {
551                    //  Add the 1st control point on the 1st curve only.
552                    cps.addLast(pts.get(0));
553                    firstPass = false;
554
555                } else {
556                    //  Determine the weight scaling factor for the next curve segment.
557                    double w2 = pts.get(0).getWeight();
558                    wfactor = w1 / w2;
559                }
560
561                //  Add the remaining control points while scaling the weights to match at the joint point.
562                int numCP = pts.size();
563                for (int j = 1; j < numCP; ++j) {
564                    ControlPoint pt = pts.get(j);
565                    double w2 = pt.getWeight();
566                    pt = pt.changeWeight(w2 * wfactor);
567                    cps.addLast(pt);
568                }
569
570                //  Keep the weight of the last point for use in the next curve.
571                w1 = cps.getLast().getWeight();
572            }
573            knots.addLast(knots.getLast());
574
575            //  Convert list of knot values into a KnotVector.
576            Float64Vector knotsV = Float64Vector.valueOf(knots);
577            knotsV = knotsV.times(1. / knots.getLast().doubleValue());  //  Normalize to a parameter length of 1.0.
578            KnotVector kv = KnotVector.newInstance(degree, knotsV);
579
580            //  Create and return the concatenated curve.
581            BasicNurbsCurve crv = BasicNurbsCurve.newInstance(cps, kv);
582            return StackContext.outerCopy(crv);
583
584        } finally {
585            StackContext.exit();
586        }
587    }
588
589    /**
590     * Connect together or "concatenate" a list of curves end-to-start. The end of each
591     * curve (u=1) is connected to the start (u=0) of the next curve. No check is made for
592     * the curve ends being near each other. The resulting curve will have the physical
593     * dimension of the highest dimension curve in the input list. The resulting curve
594     * will have the degree of the highest degree curve in the input list.
595     *
596     * @param curves The list of curves to be connected together. May not be null.
597     * @return A single {@link BasicNurbsCurve} made up of the input curves connected
598     *         together end-to-start.
599     * @throws IllegalArgumentException if there are less than 2 curves.
600     */
601    public static BasicNurbsCurve connectCurves(NurbsCurve... curves)
602            throws IllegalArgumentException {
603        FastTable<NurbsCurve> list = FastTable.newInstance();
604        list.addAll(Arrays.asList(requireNonNull(curves)));
605
606        BasicNurbsCurve crv = connectCurves(list);
607
608        return crv;
609    }
610
611    /**
612     * Return a fillet curve between two curves. Parameter values indicate the points
613     * between which the fillet is to be produced. The output curve will have the physical
614     * dimension of the highest dimension of the two input curves.
615     *
616     * @param crv1   The 1st input curve. May not be null.
617     * @param crv2   The 2nd input curve. May not be null.
618     * @param eps    The geometry tolerance to use when determining if this a planar
619     *               curve. If null is passed, then the curve would have be be exactly
620     *               planar to be a circular or conic segment.
621     * @param sEnd1  Parameter value on the first curve telling that the part of the curve
622     *               lying on this side of sFil1 shall not be replaced by the fillet.
623     * @param sFil1  Parameter value of the starting point of the fillet on the first
624     *               curve.
625     * @param sEnd2  Parameter value on the second curve telling that the part of the
626     *               curve lying on this side of sFil2 shall not be replaced by the
627     *               fillet.
628     * @param sFil2  Parameter value of the starting point of the fillet on the second
629     *               curve.
630     * @param itype  Indicator of type of fillet. = 1 - Circle, interpolating point and
631     *               tangent on first curve and the point on the 2nd curve, but not the
632     *               tangent on curve 2. = 2 - Planar conic if possible else - A generic
633     *               polynomial segment.
634     * @param degree Degree of fillet curve to return.
635     * @return A {@link BasicNurbsCurve} fillet curve.
636     */
637    @SuppressWarnings("null")
638    public static BasicNurbsCurve createFillet(Curve crv1, Curve crv2, Parameter<Length> eps,
639            double sEnd1, double sFil1, double sEnd2, double sFil2,
640            int itype, int degree) {
641
642        //  The points and tangents specified by the input parameter values
643        //  are calculated, then the interpolation curve of type given by
644        //  the  fillet indicator is calculated if possible. If a conic is
645        //  specified and not possible to make, a cubic spline segment is
646        //  made in stead.
647        
648        StackContext.enter();
649        try {
650            //  Convert both curves to the same physical dimension (the highest of the two).
651            int dim = Math.max(crv1.getPhyDimension(), crv2.getPhyDimension());
652            crv1 = crv1.toDimension(dim);
653            crv2 = crv2.toDimension(dim);
654
655            //  Get position and tangent of fillet point on the 1st curve.
656            List<GeomVector<Length>> ders = crv1.getSDerivatives(sFil1, 1);
657            GeomVector<Length> p1v = ders.get(0);
658            Point p1 = p1v.getOrigin().copyToReal();
659            GeomVector<Length> pu1 = ders.get(1);
660
661            if (sEnd1 >= sFil1) {
662                //  Reverse the direction of the tangent.
663                pu1 = pu1.opposite();
664            }
665
666            //  Get position and tangent of fillet point on the 2nd curve.
667            ders = crv2.getSDerivatives(sFil2, 1);
668            GeomVector<Length> p2v = ders.get(0);
669            Point p2 = p2v.getOrigin().copyToReal();
670            GeomVector<Length> pu2 = ders.get(1);
671
672            if (sEnd2 < sFil2) {
673                //  Reverse the direction of the tangent.
674                pu2 = pu2.opposite();
675            }
676
677            //  First, make sure p1 and p2 are not colocated.
678            //  If they are colocated, return a degenerate curve.
679            Parameter<Length> dend = p1.distance(p2);
680            if (dend.isApproxZero())
681                return StackContext.outerCopy(CurveFactory.createPoint(degree, p1));
682
683            //  Test if the points and tangents describe a planar curve or do we have a 3D curve?
684            Vector sdum = null;
685            boolean plane = true;
686            if (dim > 2) {
687                //  Find deviation from a plane by projecting a vector from p1 to p2
688                //  onto the normal formed by the cross product of the tangent vectors.
689                //  If the projection is anything other than zero, the points and tangent
690                //  vectors are not in the same plane.
691                sdum = pu1.cross(pu2).toUnitVector();
692                Vector p1p2 = p1.minus(p2).toGeomVector();
693                Parameter proj = p1p2.dot(sdum);
694                if (proj.isLargerThan(eps))
695                    plane = false;
696            }
697
698            //  Test if conic is in fact a circle.
699            if (plane && itype != 1) {
700                //  Project pu1 onto pu2.
701                Parameter dum = pu1.dot(pu2);
702                Parameter sum1 = pu1.mag();
703                Parameter sum2 = pu2.mag();
704
705                Parameter<Angle> sang;
706                if (sum1.isApproxZero() || sum2.isApproxZero())
707                    sang = Parameter.ZERO_ANGLE;
708                else {
709                    Parameter tcos = dum.divide(sum1.times(sum2));
710                    tcos = tcos.min(Parameter.ONE);
711                    sang = Parameter.acos(tcos);
712                }
713
714                if (Parameter.PI_ANGLE.minus(sang).isLessThan(ANGULAR_TOLERANCE)) {
715                    //  Opposite parallel start and end tangents.
716                    Vector norder1 = p2.minus(p1).toGeomVector();
717
718                    if (norder1.angle(pu1).minus(Parameter.HALFPI_ANGLE).abs().isLessThan(ANGULAR_TOLERANCE))
719                        itype = 1;
720
721                } else if (sang.isLargerThan(ANGULAR_TOLERANCE)) {
722                    //  Not parallel start and end tangents.
723
724                    Vector norder1, norder2;
725                    if (dim == 2) {
726                        norder1 = Vector.valueOf(pu1.get(1), pu1.get(0).opposite()).toUnitVector();
727                        norder2 = Vector.valueOf(pu2.get(1), pu2.get(0).opposite()).toUnitVector();
728
729                    } else {
730                        norder1 = sdum.cross(pu1);
731                        norder2 = sdum.cross(pu2);
732                    }
733
734                    norder1 = norder1.minus(norder2);
735                    norder2 = p2.minus(p1).toGeomVector();
736
737                    if (norder1.angle(norder2).isLessThan(ANGULAR_TOLERANCE))
738                        itype = 1;
739                }
740            }
741
742            BasicNurbsCurve crv;
743            if (plane && (itype == 2 || itype == 1)) {
744                //  The curve is planar, interpolate with a conic.
745                if (itype == 1)
746                    //  Circular arc.
747                    crv = CurveFactory.createCircularArc(p1, pu1, p2);
748                else
749                    //  Parabolic arc.
750                    crv = CurveFactory.createParabolicArc(p1, pu1.toUnitVector(), p2, pu2.toUnitVector());
751
752            } else {
753                //  Polynomial fillet.
754                crv = CurveFactory.createBlend(p1, pu1.toUnitVector(), p2, pu2.toUnitVector(), true, 0.5);
755            }
756
757            //  Elevate to the desired degree.
758            BasicNurbsCurve output = (BasicNurbsCurve)elevateToDegree(crv, degree);
759            return StackContext.outerCopy(output);
760
761        } finally {
762            StackContext.exit();
763        }
764    }
765
766    /**
767     * Return a circular fillet curve between two planar non-parallel lines. The output
768     * curve will have the physical dimension of the highest dimension of the two input
769     * lines. The corner where the fillet is applied depends on the directions of the
770     * input lines. It is the corner that is to the left of each input line when standing
771     * at the start of the line and looking toward the end of the line.
772     *
773     * @param p1     The 1st end of the 1st input line. May not be null.
774     * @param p2     The 2nd end of the 1st input line. May not be null.
775     * @param p3     The 1st end of the 2nd input line. May not be null.
776     * @param p4     The 2nd end of the 2nd input line. May not be null.
777     * @param radius The radius of the circular fillet. May not be null.
778     * @param eps    The geometry tolerance to use in determining if the input line
779     *               segments are co-planar. If null is passed, the input lines must be
780     *               exactly co-planar.
781     * @return A {@link BasicNurbsCurve} circular arc fillet curve.
782     */
783    public static BasicNurbsCurve createFillet(GeomPoint p1, GeomPoint p2, GeomPoint p3, GeomPoint p4,
784            Parameter<Length> radius, Parameter<Length> eps) {
785        requireNonNull(p1);
786        requireNonNull(p2);
787        requireNonNull(p3);
788        requireNonNull(p4);
789        requireNonNull(radius);
790
791        //  Convert both lines to the same physical dimension (the highest of the two).
792        int dim = GeomUtil.maxPhyDimension(p1, p2, p3, p4);
793        if (dim < 2)
794            throw new DimensionException(
795                    MessageFormat.format(RESOURCES.getString("dimensionNotAtLeast2"),
796                            "input points", dim));
797
798        StackContext.enter();
799        try {
800            int numDims = dim;
801            if (dim < 3)
802                dim = 3;    //  Use a min. of 3 dimension in calculations.
803            p1 = p1.toDimension(dim);
804            p2 = p2.toDimension(dim);
805            p3 = p3.toDimension(dim);
806            p4 = p4.toDimension(dim);
807
808            //  Make sure the lines are not degenerate.
809            if (p1.distance(p2).isLessThan(eps) || p3.distance(p4).isLessThan(eps)) {
810                throw new IllegalArgumentException(RESOURCES.getString("degLineInput"));
811            }
812
813            //  Extract the tangents of the two lines.
814            Vector<Length> t1 = p2.toGeomVector().minus(p1.toGeomVector());
815            Vector<Length> t2 = p4.toGeomVector().minus(p3.toGeomVector());
816
817            //  Test if the points and tangents describe co-planar lines?
818            //  Find deviation from a plane by projecting a vector from p1 to p3
819            //  onto the normal formed by the cross product of the tangent vectors.
820            //  If the projection is anything other than zero, the points and tangent
821            //  vectors are not in the same plane.
822            Vector<Dimensionless> nhat = t1.cross(t2).toUnitVector();
823            Vector p1p3 = p1.minus(p3).toGeomVector();
824            if (p1p3.mag().isApproxZero())
825                p1p3 = p1.minus(p4).toGeomVector();
826            Parameter proj = p1p3.dot(nhat);
827            if (proj.isLargerThan(eps))
828                throw new IllegalArgumentException(RESOURCES.getString("nonCoPlanarLines"));
829
830            //  Determine if the two lines are parallel.
831            Parameter a = t1.dot(t1);
832            Parameter b = t1.dot(t2);
833            Parameter c = t2.dot(t2);
834            Parameter denom = a.times(c).minus(b.times(b));     //  denom = a*c - b^2
835            if (denom.isApproxZero())
836                throw new IllegalArgumentException(RESOURCES.getString("parallelLinesInput"));
837
838            //  Offset the 1st line by "radius" amount in the direction of the 2nd line.
839            Vector<Dimensionless> yhat = nhat.cross(t1).toUnitVector();
840            Point offset = Point.valueOf(yhat.times(radius));
841            Point p1p = p1.plus(offset);
842
843            //  Offset the 2nd line by "radius" amount in the direction of the 1st line.
844            yhat = nhat.cross(t2).toUnitVector();
845            offset = Point.valueOf(yhat.times(radius));
846            Point p3p = p3.plus(offset);
847
848            //  The intersection of the two lines is the center of the fillet arc.
849            MutablePoint pc = MutablePoint.newInstance(dim, p1.getUnit());
850            GeomUtil.lineLineIntersect(p1p, t1, p3p, t2, eps, null, pc, null);
851
852            //  Find the point on the 1st line that is closest to the fillet arc center.
853            Point pa = GeomUtil.pointLineClosest(pc, p1, t1.toUnitVector());
854
855            //  Find the point on the 2nd line that is closest to the fillet arc center.
856            Point pb = GeomUtil.pointLineClosest(pc, p3, t2.toUnitVector());
857
858            //  Create a circular arc from the center and end points.
859            BasicNurbsCurve crv = CurveFactory.createCircularArcO(pc.immutable(), pa, pb);
860
861            crv = crv.toDimension(numDims); //  May need to convert to 2D at the end.
862            return StackContext.outerCopy(crv);
863
864        } finally {
865            StackContext.exit();
866        }
867    }
868
869    /**
870     * Return a blending curve between two curve ends. Parameter values indicate which
871     * ends of the input curves the blend curve is to touch. The output curve will be a
872     * conic section if possible otherwise it will be a polynomial curve and will have the
873     * physical dimension of the highest dimension of the two input curves.
874     *
875     * @param crv1   The 1st input curve. May not be null.
876     * @param crv2   The 2nd input curve. May not be null.
877     * @param eps    The geometry tolerance to use when determining if the output is a
878     *               planar curve. If null is passed, then the curve would have be be
879     *               exactly planar to be a circular or conic segment.
880     * @param sEnd1  Parameter value on crv1 indicating which end of curve 1 should start
881     *               the blend.
882     * @param sEnd2  Parameter value on crv2 indicating which end of curve 2 should end
883     *               the blend.
884     * @param degree Degree of blend curve to return.
885     * @return A {@link BasicNurbsCurve} blending curve between the input curves.
886     */
887    public static BasicNurbsCurve createBlend(Curve crv1, Curve crv2, Parameter<Length> eps,
888            double sEnd1, double sEnd2, int degree) {
889        requireNonNull(crv1);
890        requireNonNull(crv2);
891
892        //  Sort out end points.
893        double fil1 = 1, end1 = 0;
894        if (sEnd1 < 0.5) {
895            fil1 = 0;
896            end1 = 1;
897        }
898        double fil2 = 1, end2 = 0;
899        if (sEnd2 < 0.5) {
900            fil2 = 0;
901            end2 = 1;
902        }
903
904        return createFillet(crv1, crv2, eps, end1, fil1, end2, fil2, 2, degree);
905    }
906
907    /**
908     * Return a NurbsCurve, with the same degree as the input curve, which is an
909     * approximate arc-length re-parameterization of the input NurbsCurve.
910     *
911     * @param crv The NurbsCurve to be approximately re-parameterized. May not be null.
912     * @param tol The geometric tolerance to use when re-parameterizing the input curve.
913     *  May not be null or zero.
914     * @return The arc-length parameterized curve that is approximately the same as the
915     *         input curve.
916     */
917    public static BasicNurbsCurve arcLengthParameterize(NurbsCurve crv, Parameter<Length> tol) {
918        requireNonNull(tol);
919        
920        StackContext.enter();
921        try {
922            //  Grid points onto the curve.
923            PointString<SubrangePoint> str = crv.gridToTolerance(tol);
924            
925            //  Get the parametric positions of each of the points.
926            FastTable<Double> pPos = FastTable.newInstance();
927            for (SubrangePoint p : str) {
928                pPos.add(p.getParPosition().getValue(Point.X));
929            }
930            
931            //  Create a list of Bezier segment starting parametric positions (and the last end position, 1).
932            FastTable<Double> bParams = getBezierStartParameters(crv);
933            int numSegs = bParams.size() - 1;
934            
935            //  Create some points in between the Bezier start/end points and add to the list of parametric positions.
936            int p = crv.getDegree();
937            FastTable<Double> fPos = (FastTable<Double>)GridSpacing.linear(p + 3);
938            int numSubPnts = fPos.size() - 1;
939            for (int i=0; i < numSegs; ++i) {
940                double s0 = bParams.get(i);
941                double s1 = bParams.get(i+1);
942                for (int j=0; j < numSubPnts; ++j) {
943                    double s = MathTools.lineInterp(0, s0, 1, s1, fPos.get(j));
944                    pPos.add(s);
945                }
946            }
947            
948            //  Remove any exact duplicates.
949            FastSet<Double> fset = FastSet.newInstance();
950            fset.addAll(pPos);
951            pPos.clear();
952            pPos.addAll(fset);
953            
954            //  Sort the parametric position points from 0 to 1.
955            Collections.sort(pPos);
956            
957            //  Create a new set of subrange points on the curve.
958            str.clear();
959            for (double s : pPos) {
960                str.add(crv.getPoint(s));
961            }
962            
963            //  Convert the points back to a NURBS curve.
964            BasicNurbsCurve ncrv = CurveFactory.fitPoints(p, str);
965            
966            //  Remove any unnecessary knots on the new curve.
967            ncrv = thinKnotsToTolerance(ncrv, tol.divide(10));
968            
969            return StackContext.outerCopy(ncrv);
970
971        } finally {
972            StackContext.exit();
973        }
974    }
975
976    /**
977     * Attempt to minimize the number of knots in the specified NURBS curve while
978     * maintaining the same shape (to within tolerance) of the original curve.
979     *
980     * @param original The original curve to have the knots thinned to tolerance. May not be null.
981     * @param tol      The allowable tolerance between the original curve and the curve
982     *                 with the knots thinned. May not be null.
983     * @return A curve that represents the original curve to within the given tolerance
984     *         but with fewer knots. If no knots could be removed, the original curve is
985     *         returned.
986     */
987    public static BasicNurbsCurve thinKnotsToTolerance(BasicNurbsCurve original, Parameter<Length> tol) {
988
989        StackContext.enter();
990        try {
991            int degree = original.getDegree();
992            int ki = degree + 1;
993
994            //  Adjust the tolerance to conservatively make up for the multiple passes through
995            //  refinement below (where errors get stacked up).
996            Parameter<Length> tol2 = tol.divide(original.getKnotVector().length() * degree);
997
998            //  Loop over each knot and try to remove it.
999            NurbsCurve crv = original;
1000            while (ki < crv.getKnotVector().length() - degree - 1) {
1001                NurbsCurve ncrv = crv.removeKnot(ki, 3, tol2);
1002                if (crv == ncrv)
1003                    ++ki;
1004                crv = ncrv;
1005            }
1006
1007            return StackContext.outerCopy(crv.immutable());
1008            
1009        } finally {
1010            StackContext.exit();
1011        }
1012    }
1013
1014    /**
1015     * Set up a matrix containing all the binomial coefficients up to the specified row
1016     * (n) and column (k) values.
1017     *
1018     * <pre>
1019     *    bin(i,j) = (i  j)
1020     *    The binomical coefficients are defined as follow:
1021     *       (n)         n!
1022     *      ( )  =  ------------
1023     *       (k)       k!(n-k)!       0&le;k&le;n
1024     *    and the following relationship applies:
1025     *    (n+1)     (n)   ( n )
1026     *    ( k ) =   (k) + (k-1)
1027     * </pre>
1028     *
1029     * @param bin  An existing matrix that will be filled in with the binomial
1030     *             coefficients.
1031     * @param rows The number of rows (1st index) of the matrix to compute (n).
1032     * @param cols The number of columns (2nd index) of the matrix to compute (k).
1033     */
1034    private static void binomialCoef(double[][] bin, int rows, int cols) {
1035
1036        //  Make sure the array is initially zeroed.
1037        double[] row0 = bin[0];
1038        for (int j = 0; j < cols; ++j)
1039            row0[j] = 0;
1040        for (int i = 1; i < rows; ++i) {
1041            System.arraycopy(row0, 0, bin[i], 0, cols);
1042        }
1043
1044        // Setup the first line
1045        row0[0] = 1.0;
1046        //for(int k = cols-1; k > 0; --k)
1047        //  bin[0][k] = 0.0;
1048
1049        // Setup the other lines
1050        for (int n = 0; n < rows - 1; n++) {
1051            int np1 = n + 1;
1052            @SuppressWarnings("MismatchedReadAndWriteOfArray")
1053            double[] binnp1 = bin[np1];
1054            binnp1[0] = 1.0;
1055
1056            double[] binn = bin[n];
1057            for (int k = 1; k < cols; k++) {
1058                if (np1 >= k)
1059                    binnp1[k] = binn[k] + binn[k - 1];
1060                //else
1061                //  bin[n][k] = 0.0;
1062            }
1063        }
1064
1065    }
1066
1067    /**
1068     * Return true if the supplied control point polygon is simple and convex.
1069     *
1070     * @param P The list of control points to evaluate. May not be null.
1071     * @return true if the control point polygon is simple and convex.
1072     */
1073    public static boolean isSimpleConvexPolygon(List<ControlPoint> P) {
1074        //  Algorithm 1 from Ma, Hewitt, "Point inversion and projection for NURBS curve and surface", CAGD, 2003.
1075
1076        StackContext.enter();
1077        try {
1078            int n = P.size() - 1;
1079            Point P0 = P.get(0).getPoint();
1080            Point Pn = P.get(n).getPoint();
1081
1082            Point Pim1 = P.get(0).getPoint();
1083            Point Pi = P.get(1).getPoint();
1084            for (int i = 1; i < n; ++i) {
1085                Point Pip1 = P.get(i + 1).getPoint();
1086
1087                //  Calculate projection of Pi onto a line from P(i-1) to P(i+1).
1088                Point V1 = GeomUtil.pointLineClosest(Pi, Pim1, Pip1.minus(Pim1).toGeomVector().toUnitVector());
1089
1090                //  Generate V1Pi vector.
1091                Vector<Length> V1Pi = Pi.minus(V1).toGeomVector();
1092
1093                Parameter R;
1094                int no2 = n / 2;
1095                if (i < no2) {
1096                    Vector<Length> V1Pn = Pn.minus(V1).toGeomVector();
1097                    R = V1Pi.dot(V1Pn);
1098
1099                } else {
1100                    Vector<Length> V1P0 = P0.minus(V1).toGeomVector();
1101                    R = V1Pi.dot(V1P0);
1102                }
1103
1104                if (R.getValue() < -1e-14)      //  (R.getValue() < 0) accounting for typical roundoff errors.
1105                    //  The polygon is not a simple and convex one.
1106                    return false;
1107
1108                //  Prepare for next loop.
1109                Pim1 = Pi;
1110                Pi = Pip1;
1111            }
1112
1113            //  It is a simple and convex polygon
1114            return true;
1115
1116        } finally {
1117            StackContext.exit();
1118        }
1119    }
1120
1121    /**
1122     * Returns <code>true</code> if an end point is the closest point in the specified
1123     * Bezier curve segment to the target point. Returns <code>false</code> if the closest
1124     * point may be interior to the curve segment.
1125     *
1126     * @param point The point to find the closest/farthest point on this curve to/from.
1127     *              May not be null.
1128     * @param P     The control point polygon for the Bezier curve. May not be null.
1129     * @return <code>true</code> if an end point of the specified segment is closest to
1130     *         the target point.
1131     */
1132    public static boolean isBezierEndPointNearest(GeomPoint point, List<ControlPoint> P) {
1133        //  Algorithm 3 from Ma, Hewitt, "Point inversion and projection for NURBS curve and surface", CAGD, 2003
1134        //  (Logic reversed).
1135
1136        StackContext.enter();
1137        try {
1138            //  Extract points we need from the control polygon.
1139            Point P0 = P.get(0).getPoint();
1140            Point P1 = P.get(1).getPoint();
1141            int n = P.size() - 1;
1142            Point Pnm1 = P.get(n - 1).getPoint();
1143            Point Pn = P.get(n).getPoint();
1144
1145            //  Form vectors.
1146            Vector<Length> P0P = point.minus(P0).toGeomVector();
1147            Vector<Length> P0P1 = P1.minus(P0).toGeomVector();
1148            Vector<Length> PPn = Pn.minus(point).toGeomVector();
1149            Vector<Length> Pnm1Pn = Pn.minus(Pnm1).toGeomVector();
1150            Vector<Length> PnP0 = P0.minus(Pn).toGeomVector();
1151            Vector<Length> PnP = PPn.opposite();
1152
1153            //  Calculate the dot products.
1154            Parameter R1 = P0P.dot(P0P1);
1155            Parameter R2 = PPn.dot(Pnm1Pn);
1156            Parameter R3 = PnP0.dot(PnP);
1157            Parameter R4 = PnP0.dot(P0P);
1158
1159            return R1.getValue() < 0 || R2.getValue() < 0 && R3.times(R4).getValue() > 0;
1160
1161        } finally {
1162            StackContext.exit();
1163        }
1164    }
1165
1166    /**
1167     * Return the union between two knot vectors. This merges two knot vectors together to
1168     * obtain on combined knot vector. Warning: The results of this method are useless if
1169     * the knot vectors being merged are not from NURBS curves of the same degree.
1170     *
1171     * @param Ua The 1st knot vector to merge. May not be null.
1172     * @param Ub The 2nd knot vector to merge. May not be null.
1173     * @return The union of Ua and Ub.
1174     */
1175    public static Float64Vector knotVectorUnion(Float64Vector Ua, Float64Vector Ub) {
1176
1177        StackContext.enter();
1178        try {
1179            int UaSize = Ua.getDimension();
1180            int UbSize = Ub.getDimension();
1181            FastTable<Float64> U = FastTable.newInstance();
1182            U.setSize(UaSize + UbSize);
1183
1184            boolean done = false;
1185            int i = 0, ia = 0, ib = 0;
1186            while (!done) {
1187                double Uav = Ua.getValue(ia);
1188                double Ubv = Ub.getValue(ib);
1189                double t;
1190                if (MathTools.isApproxEqual(Uav,Ubv)) {
1191                    t = Uav;
1192                    ++ia;
1193                    ++ib;
1194                } else if (Uav < Ubv) {
1195                    t = Uav;
1196                    ++ia;
1197                } else {
1198                    t = Ubv;
1199                    ++ib;
1200                }
1201                U.set(i++, Float64.valueOf(t));
1202                done = (ia >= UaSize || ib >= UbSize);
1203            }
1204            U.setSize(i);
1205
1206            //  Create the merged knot vector instance.
1207            Float64Vector Uv = Float64Vector.valueOf(U);
1208            return StackContext.outerCopy(Uv);
1209
1210        } finally {
1211            StackContext.exit();
1212        }
1213    }
1214
1215    /**
1216     * Allocate a recyclable 2D array of double values using factory methods.
1217     * <p>
1218     * WARNING: The array returned may not be zeroed and may be larger than the requested
1219     * number of rows &amp; columns! The array returned by this method can be recycled by
1220     * recycle2DArray().
1221     * </p>
1222     *
1223     * @param rows The minimum number of rows (1st index) for the returned array.
1224     * @param cols The minimum number of columns (2nd index) for the returned array.
1225     * @return A 2D array of doubles allocated using factory methods.
1226     * @see #recycle2DArray
1227     */
1228    public static double[][] allocate2DArray(int rows, int cols) {
1229        double[][] arr = DOUBLEARRAY_FACTORY.array(rows);
1230        for (int i = 0; i < rows; ++i) {
1231            arr[i] = ArrayFactory.DOUBLES_FACTORY.array(cols);
1232        }
1233        return arr;
1234    }
1235
1236    /**
1237     * Recycle any 2D array of doubles that was created by allocate2DArray().
1238     *
1239     * @param arr The array to be recycled. The array must have been created by this the
1240     *            allocate2DArray() method! May not be null.
1241     * @see #allocate2DArray
1242     */
1243    public static void recycle2DArray(double[][] arr) {
1244        int length = arr.length;
1245        for (int i = 0; i < length; ++i) {
1246            if (arr[i] != null) {
1247                ArrayFactory.DOUBLES_FACTORY.recycle(arr[i]);
1248                arr[i] = null;
1249            }
1250        }
1251        DOUBLEARRAY_FACTORY.recycle(arr);
1252    }
1253
1254    /**
1255     * Factory that creates recyclable arrays that can contain ControlPoint objects.
1256     */
1257    private static final ArrayFactory<ControlPoint[]> CONTROLPOINTARRAY_FACTORY = new ArrayFactory<ControlPoint[]>() {
1258        @Override
1259        protected ControlPoint[] create(int size) {
1260            return new ControlPoint[size];
1261        }
1262    };
1263
1264    /**
1265     * Factory that creates recyclable arrays of doubles.
1266     */
1267    private static final ArrayFactory<double[][]> DOUBLEARRAY_FACTORY = new ArrayFactory<double[][]>() {
1268        @Override
1269        protected double[][] create(int size) {
1270            return new double[size][];
1271        }
1272    };
1273}