Source code

001/**
002 * AbstractCurve -- Represents the implementation in common to all curves in nD space.
003 *
004 * Copyright (C) 2009-2025, 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.1 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 Lesser 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 */
018package geomss.geom;
019
020import jahuwaldt.js.param.Parameter;
021import jahuwaldt.tools.math.*;
022import static jahuwaldt.tools.math.MathTools.isApproxEqual;
023import java.text.MessageFormat;
024import java.util.Collections;
025import java.util.Comparator;
026import java.util.List;
027import static java.util.Objects.isNull;
028import static java.util.Objects.requireNonNull;
029import java.util.logging.Level;
030import java.util.logging.Logger;
031import javax.measure.quantity.Area;
032import javax.measure.quantity.Dimensionless;
033import javax.measure.quantity.Length;
034import javax.measure.unit.SI;
035import javax.measure.unit.Unit;
036import javolution.context.ObjectFactory;
037import javolution.context.StackContext;
038import javolution.lang.MathLib;
039import javolution.util.FastList;
040import javolution.util.FastTable;
041
042/**
043 * The interface and implementation in common to all curves.
044 *
045 * <p> Modified by: Joseph A. Huwaldt </p>
046 *
047 * @author Joseph A. Huwaldt, Date: May 15, 2009
048 * @version February 18, 2025
049 *
050 * @param <T> The sub-type of this AbstractCurve.
051 */
052@SuppressWarnings({"serial", "CloneableImplementsClone"})
053public abstract class AbstractCurve<T extends AbstractCurve> extends AbstractGeomElement<T> implements Curve<T> {
054
055    /*
056     *  Reference 1:  Michael E. Mortenson, "Geometric Modeling, Third Edition", ISBN:9780831132989, 2008.
057     */
058    
059    /**
060     * Generic/default tolerance for use in root finders, etc.
061     */
062    public static final double GTOL = 1e-6;
063
064    /**
065     * Tolerance allowed on parametric spacing being == 0 or == 1.
066     */
067    protected static final double TOL_S = Parameter.EPS;
068
069    /**
070     * Returns the number of parametric dimensions of the geometry element. This
071     * implementation always returns 1.
072     */
073    @Override
074    public int getParDimension() {
075        return 1;
076    }
077
078    /**
079     * Checks the input parameter for validity (in the range 0 to 1 for example) and
080     * throws an exception if it is not valid.
081     *
082     * @param s parametric distance (0.0 to 1.0 inclusive).
083     * @throws IllegalArgumentException if there is any problem with the parameter value.
084     */
085    protected void validateParameter(double s) {
086        if (s < -TOL_S || s > 1 + TOL_S)
087            throw new IllegalArgumentException(
088                    MessageFormat.format(RESOURCES.getString("crvInvalidSValue"), s));
089    }
090
091    /**
092     * Checks the input parameter for validity (in the range 0 to 1 for example) and
093     * throws an exception if it is not valid.
094     *
095     * @param s parametric distance. Must be a 1-dimensional point with a value in the
096     *          range 0 to 1.0. Units are ignored.
097     * @throws IllegalArgumentException if there is any problem with the parameter value.
098     */
099    protected void validateParameter(GeomPoint s) {
100        if (isNull(s) || s.getPhyDimension() != 1)
101            throw new IllegalArgumentException(MessageFormat.format(
102                    RESOURCES.getString("invalidParamDim"), "Curve", 1));
103        validateParameter(s.getValue(0));
104    }
105
106    /**
107     * If the input parameter is within tol of 0 or 1, or outside the bounds of 0 or 1,
108     * then the output value will be set to exactly 0 or 1. Otherwise the input value will
109     * be returned.
110     *
111     * @param s The parameter to check for approximate end values.
112     * @return The input value or 0 and 1 if the input value is within tolerance of the
113     *         ends.
114     */
115    protected static final double roundParNearEnds(double s) {
116        return (s < TOL_S ? 0. : (s > 1. - TOL_S ? 1. : s));
117    }
118    
119    /**
120     * Return true if the input parameter is within tolerance of or less than the start of the
121     * parameter range (s &le; tol).
122     *
123     * @param s   The parameter to check.
124     * @param tol The tolerance to use for the check.
125     * @return true if the parameter is within tolerance of zero.
126     */
127    protected static final boolean parNearStart(double s, double tol) {
128        return s <= tol;
129    }
130
131    /**
132     * Return true if the input parameter is within tolerance of or greater than the end
133     * of the parameter range (s &ge; 1.0 - tol).
134     *
135     * @param s   The parameter to check.
136     * @param tol The tolerance to use for the check.
137     * @return true if the parameter is within tolerance of 1.0.
138     */
139    protected static final boolean parNearEnd(double s, double tol) {
140        return s >= 1.0 - tol;
141    }
142
143    /**
144     * Return true if the input parameter is within tolerance of or outside of the bounds
145     * of either the end of the parameter range (s &le; tol || s &ge; 1.0 - tol).
146     *
147     * @param s   The parameter to check.
148     * @param tol The tolerance to use for the check.
149     * @return true if the parameter is within tolerance of 0.0 or 1.0.
150     */
151    protected static final boolean parNearEnds(double s, double tol) {
152        return s <= tol || s >= 1.0 - tol;
153    }
154
155    /**
156     * Split this curve at the specified parametric position returning a list containing
157     * two curves (a lower curve with smaller parametric positions than "s" and an upper
158     * curve with larger parametric positions).
159     *
160     * @param s The parametric position where this curve should be split (must not be 0 or
161     *          1!). Must be a 1-dimensional point with a value in the range 0 to 1.0,
162     *          non-inclusive. Units are ignored.
163     * @return A list containing two curves: 0 == the lower curve, 1 == the upper curve.
164     */
165    @Override
166    public GeomList<T> splitAt(GeomPoint s) {
167        return splitAt(s.getValue(0));
168    }
169
170    /**
171     * Return a subrange point on the curve for the given parametric distance along the
172     * curve.
173     *
174     * @param s parametric distance to calculate a point for (0.0 to 1.0 inclusive).
175     * @return The subrange point or <code>null</code> if the parameter is not in a valid
176     *         range.
177     */
178    @Override
179    public SubrangePoint getPoint(double s) {
180        validateParameter(s);
181        s = roundParNearEnds(s);
182        return SubrangePoint.newInstance(this, Point.valueOf(s));
183    }
184
185    /**
186     * Return a subrange point on the curve for the given parametric distance along the
187     * curve.
188     *
189     * @param s parametric distance to calculate a point for. Must be a 1-dimensional
190     *          point with a value in the range 0 to 1.0. Units are ignored.
191     * @return The subrange point
192     */
193    @Override
194    public SubrangePoint getPoint(GeomPoint s) {
195        validateParameter(s);
196        return SubrangePoint.newInstance(this, s);
197    }
198
199    /**
200     * Calculate a point on the curve for the given parametric distance along the curve.
201     *
202     * @param s parametric distance to calculate a point for. Must be a 1-dimensional
203     *          point with a value in the range 0 to 1.0. Units are ignored.
204     * @return calculated point
205     */
206    @Override
207    public Point getRealPoint(GeomPoint s) {
208        validateParameter(s);
209        return getRealPoint(s.getValue(0));
210    }
211
212    /**
213     * Calculate all the derivatives from <code>0</code> to <code>grade</code> with
214     * respect to parametric position(s) on a parametric object for the given parametric
215     * position on the object, <code>d^{grade}p(s)/d^{grade}s</code>.
216     * <p>
217     * Example:<br>
218     * 1st derivative (grade = 1), this returns <code>[p(s), dp(s)/ds]</code>;<br>
219     * 2nd derivative (grade = 2), this returns <code>[p(s), dp(s)/ds, d^2p(s)/d^2s]</code>; etc.
220     * </p>
221     *
222     * @param s     parametric position to calculate the derivatives for. Must match the
223     *              parametric dimension of this parametric surface and have each value in
224     *              the range 0 to 1.0. Units are ignored.
225     * @param grade The maximum grade to calculate the derivatives for (1=1st derivative,
226     *              2=2nd derivative, etc)
227     * @return A list of lists of derivatives (one list for each parametric dimension).
228     *         Each list contains derivatives up to the specified grade at the specified
229     *         parametric position. Example: for a surface list element 0 is the array of
230     *         derivatives in the 1st parametric dimension (s), then 2nd list element is
231     *         the array of derivatives in the 2nd parametric dimension (t).
232     * @throws IllegalArgumentException if the grade is &lt; 0.
233     */
234    @Override
235    public List<List<Vector<Length>>> getDerivatives(GeomPoint s, int grade) {
236        validateParameter(s);
237
238        List<List<Vector<Length>>> list = FastTable.newInstance();
239        list.add(getSDerivatives(s.getValue(0), grade));
240
241        return list;
242    }
243
244    /**
245     * Calculate all the derivatives from <code>0</code> to <code>grade</code> with
246     * respect to parametric distance on the curve for the given parametric distance along
247     * the curve, <code>d^{grade}p(s)/d^{grade}s</code>.
248     * <p>
249     * Example:<br>
250     * 1st derivative (grade = 1), this returns <code>[p(s), dp(s)/ds]</code>;<br>
251     * 2nd derivative (grade = 2), this returns <code>[p(s), dp(s)/ds, d^2p(s)/d^2s]</code>; etc.
252     * </p>
253     *
254     * @param s     parametric distance to calculate the derivatives for. Must be a
255     *              1-dimensional point with a value in the range 0 to 1.0. Units are
256     *              ignored.
257     * @param grade The maximum grade to calculate the derivatives for (1=1st derivative,
258     *              2=2nd derivative, etc)
259     * @return A list of derivatives up to the specified grade of the curve at the
260     *         specified parametric position.
261     * @throws IllegalArgumentException if the grade is &lt; 0.
262     */
263    @Override
264    public List<Vector<Length>> getSDerivatives(GeomPoint s, int grade) {
265        validateParameter(s);
266        return getSDerivatives(s.getValue(0), grade);
267    }
268
269    /**
270     * Calculate a derivative with respect to parametric distance of the given grade on
271     * the curve for the given parametric distance along the curve,
272     * <code>d^{grade}p(s)/d^{grade}s</code>.
273     * <p>
274     * Example:<br>
275     * 1st derivative (grade = 1), this returns <code>dp(s)/ds</code>; <br>
276     * 2nd derivative (grade = 2), this returns <code>d^2p(s)/d^2s</code>; etc.
277     * </p>
278     * <p>
279     * Note: Cartesian space derivatives can be calculated as follows:
280     * <pre>
281     *      derivative = curve.getSDerivative(s, 1);
282     *      dy/dx = (dy/ds)/(dx/ds) = derivative.getValue(1)/derivative.getValue(0),
283     *      dy/dz = (dy/ds)/(dz/ds) = derivative.getValue(1)/derivative.getValue(2),
284     *      etc
285     * </pre>
286     *
287     * @param s     Parametric distance to calculate a derivative for (0.0 to 1.0
288     *              inclusive).
289     * @param grade The grade to calculate the derivative for (1=1st derivative, 2=2nd
290     *              derivative, etc)
291     * @return The specified derivative of the curve at the specified parametric position.
292     * @throws IllegalArgumentException if the grade is &lt; 0.
293     */
294    @Override
295    public Vector<Length> getSDerivative(double s, int grade) {
296        validateParameter(s);
297        StackContext.enter();
298        try {
299
300            List<Vector<Length>> ders = getSDerivatives(s, grade);
301            Vector<Length> derivative = ders.get(grade);
302            return StackContext.outerCopy(derivative);
303
304        } finally {
305            StackContext.exit();
306        }
307    }
308
309    /**
310     * Return the tangent vector for the given parametric distance along the curve. This
311     * vector contains the normalized 1st derivative of the curve with respect to
312     * parametric distance in each component direction: dx/ds/|dp(s)/ds|,
313     * dy/ds/|dp(s)/ds|, dz/ds/|dp(s)/ds|, etc.
314     * <p>
315     * Note: Cartesian space derivatives can be calculated as follows:
316     * <pre>
317     *      tangent = curve.getTangent(s);
318     *      dy/dx = (dy/ds)/(dx/ds) = tangent.getValue(1)/tangent.getValue(0),
319     *      dy/dz = (dy/ds)/(dz/ds) = tangent.getValue(1)/tangent.getValue(2),
320     *      etc
321     * </pre>
322     *
323     * @param s Parametric distance to calculate a tangent vector for (0.0 to 1.0
324     *          inclusive).
325     * @return The tangent vector of the curve at the specified parametric position.
326     */
327    @Override
328    public Vector<Dimensionless> getTangent(double s) {
329        validateParameter(s);
330        StackContext.enter();
331        try {
332
333            //  Ref. 1, Eqn. 12.1.
334            List<Vector<Length>> ders = getSDerivatives(s, 1);
335            Vector<Dimensionless> derivative = ders.get(1).toUnitVector();
336
337            return StackContext.outerCopy(derivative);
338
339        } finally {
340            StackContext.exit();
341        }
342    }
343
344    /**
345     * Return the tangent vector for the given parametric distance along the curve. This
346     * vector contains the normalized 1st derivative of the curve with respect to
347     * parametric distance in each component direction: dx/ds/|dp(s)/ds|,
348     * dy/ds/|dp(s)/ds|, dz/ds/|dp(s)/ds|, etc.
349     * <p>
350     * Note: Cartesian space derivatives can be calculated as follows:
351     * <pre>
352     *      tangent = curve.getTangent(s);
353     *      dy/dx = (dy/ds)/(dx/ds) = tangent.getValue(1)/tangent.getValue(0),
354     *      dy/dz = (dy/ds)/(dz/ds) = tangent.getValue(1)/tangent.getValue(2),
355     *      etc
356     * </pre>
357     *
358     * @param s Parametric distance to calculate a tangent vector for. Must be a
359     *          1-dimensional point with a value in the range 0 to 1.0. Units are ignored.
360     * @return The tangent vector of the curve at the specified parametric position.
361     */
362    @Override
363    public Vector<Dimensionless> getTangent(GeomPoint s) {
364        validateParameter(s);
365        return getTangent(s.getValue(0));
366    }
367
368    /**
369     * Return the principal normal vector for the given parametric distance,
370     * <code>s</code>, along the curve. This vector is normal to the curve, lies in the
371     * normal plane (is perpendicular to the tangent vector), and points toward the center
372     * of curvature at the specified point.
373     *
374     * @param s Parametric distance to calculate the principal normal vector for (0.0 to
375     *          1.0 inclusive).
376     * @return The principal normal vector of the curve at the specified parametric
377     *         position.
378     * @throws IllegalArgumentException if the curve does not have a 2nd derivative with
379     * respect to <code>s</code> (if the curve is locally a straight line segment -- this
380     * means there are an infinite number of possible principal normal vectors).
381     */
382    @Override
383    public Vector<Dimensionless> getPrincipalNormal(double s) {
384        validateParameter(s);
385
386        StackContext.enter();
387        try {
388            //  Ref. 1, Eqn. 12.15.
389            List<Vector<Length>> ders = getSDerivatives(s, 2);
390            Vector<Length> piu = ders.get(1);
391            Vector<Length> piuu = ders.get(2);
392
393            if (piuu.mag().isApproxZero())
394                throw new IllegalArgumentException(RESOURCES.getString("undefinedPrincipleNormal"));
395
396            //  ni = k/|k|;  k = piuu - piuu DOT piu / |piu|^2 * piu;
397            Parameter piuMag = piu.mag();
398            Parameter piuuDotpiu = piuu.dot(piu);
399            Vector k = piu.times(piuuDotpiu.divide(piuMag).divide(piuMag));
400            k = piuu.minus(k);
401
402            return StackContext.outerCopy(k.toUnitVector());
403
404        } finally {
405            StackContext.exit();
406        }
407    }
408
409    /**
410     * Return the principal normal vector for the given parametric distance,
411     * <code>s</code>, along the curve. This vector is normal to the curve, lies in the
412     * normal plane (is perpendicular to the tangent vector), and points toward the center
413     * of curvature at the specified point.
414     *
415     * @param s Parametric distance to calculate the principal normal vector for. Must be
416     *          a 1-dimensional point with a value in the range 0 to 1.0. Units are
417     *          ignored.
418     * @return The principal normal vector of the curve at the specified parametric
419     *         position.
420     * @throws IllegalArgumentException if the curve does not have a 2nd derivative with
421     * respect to <code>s</code> (if the curve is locally a straight line segment -- this
422     * means there are an infinite number of possible principal normal vectors).
423     */
424    @Override
425    public Vector<Dimensionless> getPrincipalNormal(GeomPoint s) {
426        validateParameter(s);
427        return getPrincipalNormal(s.getValue(0));
428    }
429
430    /**
431     * Return the binormal vector for the given parametric distance along the curve. This
432     * vector is normal to the curve at <code>s</code>, lies in the normal plane (is
433     * perpendicular to the tangent vector), and is perpendicular to the principal normal
434     * vector. Together, the tangent vector, principal normal vector, and binormal vector
435     * form a useful orthogonal frame.
436     *
437     * @param s Parametric distance to calculate the binormal vector for (0.0 to 1.0
438     *          inclusive).
439     * @return The binormal vector of the curve at the specified parametric position.
440     * @throws IllegalArgumentException if the curve does not have a 2nd derivative with
441     * respect to <code>s</code> (if the curve is locally a straight line segment).
442     * @throws DimensionException if the curve does not have at least 3 physical
443     * dimensions.
444     */
445    @Override
446    public Vector<Dimensionless> getBinormal(double s) throws DimensionException {
447        validateParameter(s);
448        if (getPhyDimension() < 3)
449            throw new DimensionException(RESOURCES.getString("undefinedBinormal"));
450
451        StackContext.enter();
452        try {
453            
454            //  Ref. 1, Eqn. 12.18.
455            //  bi = ti X ni;
456            Vector ni = getPrincipalNormal(s);
457            Vector ti = getTangent(s);
458            Vector bi = ti.cross(ni);
459            return StackContext.outerCopy(bi);
460            
461        } finally {
462            StackContext.exit();
463        }
464    }
465
466    /**
467     * Return the binormal vector for the given parametric distance along the curve. This
468     * vector is normal to the curve at <code>s</code>, lies in the normal plane (is
469     * perpendicular to the tangent vector), and is perpendicular to the principal normal
470     * vector. Together, the tangent vector, principal normal vector, and binormal vector
471     * form a useful orthogonal frame.
472     *
473     * @param s Parametric distance to calculate the binormal vector for. Must be a
474     *          1-dimensional point with a value in the range 0 to 1.0. Units are ignored.
475     * @return The binormal vector of the curve at the specified parametric position.
476     * @throws IllegalArgumentException if the curve does not have a 2nd derivative with
477     * respect to <code>s</code> (if the curve is locally a straight line segment).
478     * @throws DimensionException if the curve does not have at least 3 physical
479     * dimensions.
480     */
481    @Override
482    public Vector<Dimensionless> getBinormal(GeomPoint s) throws DimensionException {
483        validateParameter(s);
484        return getBinormal(s.getValue(0));
485    }
486
487    /**
488     * Return the curvature (kappa = 1/rho; where rho = the radius of curvature) of the
489     * curve at the parametric position <code>s</code>. The curvature vector (vector from
490     * p(s) to the center of curvature) can be constructed from:  <code>k = rhoi*ni</code>
491     * or <code>k = curve.getPrincipalNormal(s).divide(curve.getCurvature(s));</code>
492     *
493     * @param s Parametric distance to calculate the curvature for (0.0 to 1.0 inclusive).
494     * @return The curvature of the curve at the specified parametric position in units of
495     *         1/Length.
496     */
497    @Override
498    public Parameter getCurvature(double s) {
499        validateParameter(s);
500        int dim = getPhyDimension();
501        if (dim < 2)
502            return Parameter.valueOf(0, getUnit().inverse());
503
504        StackContext.enter();
505        try {
506            //  Ref. 1, Eqn. 12.24.
507            List<Vector<Length>> ders = getSDerivatives(s, 2);
508            Vector<Length> piu = ders.get(1);
509            Vector<Length> piuu = ders.get(2);
510
511            Parameter piuMag = piu.mag();
512            if (piuu.mag().isApproxZero() || piuMag.isApproxZero())
513                return StackContext.outerCopy(Parameter.valueOf(0, getUnit().inverse()));
514
515            Parameter kappa;
516
517            if (dim > 2) {
518                //  kappa = |piu X piuu| / |piu|^3
519                Vector piuXpiuu = piu.cross(piuu);
520                kappa = piuXpiuu.mag().divide(piuMag.pow(3));
521
522            } else {
523                Parameter xu = piu.get(0);
524                Parameter yu = piu.get(1);
525                Parameter xuu = piuu.get(0);
526                Parameter yuu = piuu.get(1);
527
528                //  kappa = |xu*yuu - xuu*yu|/(xu^2 + yu^2)^(3/2) = |xu*yuu - xuu*yu|/|piu|^3
529                Parameter num = xu.times(yuu).minus(xuu.times(yu)).abs();
530                kappa = num.divide(piuMag.pow(3));
531            }
532
533            kappa = kappa.to(getUnit().inverse());
534            return StackContext.outerCopy(kappa);
535
536        } finally {
537            StackContext.exit();
538        }
539    }
540
541    /**
542     * Return the curvature (kappa = 1/rho; where rho = the radius of curvature) of the
543     * curve at the parametric position <code>s</code>. The curvature vector (vector from
544     * p(s) to the center of curvature) can be constructed from:  <code>k = rhoi*ni</code>
545     * or <code>k = curve.getPrincipalNormal(s).divide(curve.getCurvature(s));</code>
546     *
547     * @param s Parametric distance to calculate the curvature for. Must be a
548     *          1-dimensional point with a value in the range 0 to 1.0. Units are ignored.
549     * @return The curvature of the curve at the specified parametric position in units of
550     *         1/Length.
551     */
552    @Override
553    public Parameter getCurvature(GeomPoint s) {
554        validateParameter(s);
555        return getCurvature(s.getValue(0));
556    }
557
558    /**
559     * Return the variation of curvature or rate of change of curvature (VOC or
560     * dKappa(s)/ds) at the parametric position <code>s</code>.
561     *
562     * @param s Parametric distance to calculate the variation of curvature for (0.0 to
563     *          1.0 inclusive).
564     * @return The variation of curvature of the curve at the specified parametric
565     *         position in units of 1/Length.
566     */
567    @Override
568    public Parameter getVariationOfCurvature(double s) {
569        validateParameter(s);
570        int dim = getPhyDimension();
571        if (dim < 2)
572            return Parameter.valueOf(0, getUnit().inverse());
573
574        //  Ref: Moreton, H.P., "Minimum Curvature Variation Curves, Networks,
575        //  and Surfaces for Fair Free-Form Shape Design", Dissertation,
576        //  University of California Berkely, 1992, pg. 79.
577        
578        StackContext.enter();
579        try {
580            //  Get the needed derivatives.
581            List<Vector<Length>> ders = getSDerivatives(s, 3);
582            Vector<Length> piu = ders.get(1);
583            Vector<Length> piuu = ders.get(2);
584            Vector<Length> piuuu = ders.get(3);
585
586            Parameter piuMag = piu.mag();
587            if (piuMag.isApproxZero())
588                return StackContext.outerCopy(Parameter.valueOf(0, getUnit().inverse()));
589
590            Parameter dKds;
591            if (dim > 2) {
592                //  u = piu x piuu
593                Vector u = piu.cross(piuu);
594
595                //  kappa = (piu X piuu) / |piu|^3
596                Vector kappa = u.divide(piuMag.pow(3));
597
598                //  du = piu x piuuu
599                Vector du = piu.cross(piuuu);
600
601                //  v = |piu|^3 = (piu dot piu)^(3/2)
602                Parameter v = piuMag.pow(3);
603
604                //  dv = 3*|piu|*(piu dot piuu)
605                Parameter dv = piu.dot(piuu).times(piuMag).times(3);
606
607                //  dKv/ds = (v*du - u*dv)/v^2
608                Vector num = du.times(v).minus(u.times(dv));
609                Vector dKvds = num.divide(v.pow(2));
610
611                //  dK/ds = (kappa/|kappa|) dot dKv/ds
612                dKds = dKvds.dot(kappa.toUnitVector());
613
614            } else {
615                Parameter xu = piu.get(0);
616                Parameter yu = piu.get(1);
617                Parameter xuu = piuu.get(0);
618                Parameter yuu = piuu.get(1);
619                Parameter xuuu = piuuu.get(0);
620                Parameter yuuu = piuuu.get(1);
621
622                //  u = xu*yuu - xuu*yu
623                Parameter u = xu.times(yuu).minus(xuu.times(yu));
624
625                //  du = xu*yuuu - xuuu*yu
626                Parameter du = xu.times(yuuu).minus(xuuu.times(yu));
627
628                //  v = |piu|^3 = (piu dot piu)^(3/2)
629                Parameter v = piuMag.pow(3);
630
631                //  dv = 3*|piu|*(piu dot piuu)
632                Parameter dv = xu.times(xuu).plus(yu.times(yuu)).times(piuMag).times(3);
633
634                //  dK/ds = (v*du - u*dv)/v^2
635                dKds = v.times(du).minus(u.times(dv)).divide(v.pow(2));
636            }
637
638            Parameter voc = dKds.to(getUnit().inverse());
639            return StackContext.outerCopy(voc);
640
641        } finally {
642            StackContext.exit();
643        }
644    }
645
646    /**
647     * Return the variation of curvature or rate of change of curvature (VOC or
648     * dKappa(s)/ds) at the parametric position <code>s</code>.
649     *
650     * @param s Parametric distance to calculate the variation of curvature for. Must be a
651     *          1-dimensional point with a value in the range 0 to 1.0. Units are ignored.
652     * @return The variation of curvature of the curve at the specified parametric
653     *         position in units of 1/Length.
654     */
655    @Override
656    public Parameter getVariationOfCurvature(GeomPoint s) {
657        validateParameter(s);
658        return getVariationOfCurvature(s.getValue(0));
659    }
660
661    /**
662     * Return the torsion of the curve at the parametric position <code>s</code>. The
663     * torsion is a measure of the rotation or twist about the tangent vector.
664     *
665     * @param s Parametric distance to calculate the torsion for (0.0 to 1.0 inclusive).
666     * @return The torsion of the curve at the specified parametric position in units of
667     *         1/Length.
668     */
669    @Override
670    public Parameter getTorsion(double s) {
671        validateParameter(s);
672        int dim = getPhyDimension();
673        if (dim < 3)
674            return Parameter.valueOf(0, getUnit().inverse());
675
676        StackContext.enter();
677        try {
678            //  Ref. 1, Eqn. 12.30.
679            List<Vector<Length>> ders = getSDerivatives(s, 3);
680            Vector<Length> piu = ders.get(1);
681            Vector<Length> piuu = ders.get(2);
682            Vector<Length> piuuu = ders.get(3);
683
684            if (piuuu.mag().isApproxZero() || piuu.mag().isApproxZero()
685                    || piu.mag().isApproxZero())
686                return StackContext.outerCopy(Parameter.valueOf(0, getUnit().inverse()));
687
688            //  tau = piu DOT (piuu X piuuu) / |piu X piuu|^2
689            Parameter piuXpiuu = piu.cross(piuu).mag();
690            Parameter tau = piu.dot(piuu.cross(piuuu)).divide(piuXpiuu.times(piuXpiuu));
691
692            return StackContext.outerCopy(tau);
693
694        } finally {
695            StackContext.exit();
696        }
697    }
698
699    /**
700     * Return the torsion of the curve at the parametric position <code>s</code>. The
701     * torsion is a measure of the rotation or twist about the tangent vector.
702     *
703     * @param s Parametric distance to calculate the torsion for. Must be a 1-dimensional
704     *          point with a value in the range 0 to 1.0. Units are ignored.
705     * @return The torsion of the curve at the specified parametric position in units of
706     *         1/Length.
707     */
708    @Override
709    public Parameter getTorsion(GeomPoint s) {
710        validateParameter(s);
711        return getTorsion(s.getValue(0));
712    }
713
714    /**
715     * Return the total arc length of this curve.
716     *
717     * @param eps The desired fractional accuracy on the arc-length.
718     * @return The total arc length of the curve.
719     */
720    @Override
721    public Parameter<Length> getArcLength(double eps) {
722        return getArcLength(0, 1, eps);
723    }
724
725    /**
726     * Return the arc length of a segment of this curve.
727     *
728     * @param s1  The starting parametric distance along the curve to begin the arc-length
729     *            calculation from.
730     * @param s2  The ending parametric distance along the curve to end the arc-length
731     *            calculation at.
732     * @param eps The desired fractional accuracy on the arc-length.
733     * @return The arc length of the specified segment of the curve.
734     */
735    @Override
736    public Parameter<Length> getArcLength(double s1, double s2, double eps) {
737        validateParameter(s1);
738        validateParameter(s2);
739
740        //  If s1 and s2 are equal, then we know the answer already.
741        if (isApproxEqual(s1, s2, TOL_S))
742            return Parameter.ZERO_LENGTH.to(getUnit());
743
744        //  Make sure that s2 is > s1.
745        if (s1 > s2) {
746            double temp = s1;
747            s1 = s2;
748            s2 = temp;
749        }
750
751        ArcLengthEvaluatable arcLengthEvaluator = ArcLengthEvaluatable.newInstance(this);
752        double length = 0;
753        try {
754
755            //  Integrate to find the arc length:   L = int_s1^s2{ sqrt(ps(s) DOT ps(s)) ds}
756            length = Quadrature.adaptLobatto(arcLengthEvaluator, s1, s2, eps);
757
758        } catch (IntegratorException | RootException e) {
759            //  Fall back on Gaussian Quadrature.
760            Logger.getLogger(AbstractCurve.class.getName()).log(Level.WARNING,
761                    "Trying Gaussian Quadrature in getArcLength()...");
762            try {
763                length = Quadrature.gaussLegendre_Wx1N20(arcLengthEvaluator, s1, s2);
764
765            } catch (RootException err) {
766                Logger.getLogger(AbstractCurve.class.getName()).log(Level.SEVERE,
767                        "Could not find arc-length.", err);
768            }
769
770        } finally {
771            ArcLengthEvaluatable.recycle(arcLengthEvaluator);
772        }
773
774        Parameter<Length> L = Parameter.valueOf(length, getUnit());
775
776        return L;
777    }
778
779    /**
780     * Return the arc length of a segment of this curve.
781     *
782     * @param s1  The starting parametric distance along the curve to begin the arc-length
783     *            calculation from. Must be a 1-dimensional point with a value in the
784     *            range 0 to 1.0. Units are ignored.
785     * @param s2  The ending parametric distance along the curve to end the arc-length
786     *            calculation at. Must be a 1-dimensional point with a value in the range
787     *            0 to 1.0. Units are ignored.
788     * @param eps The desired fractional accuracy on the arc-length.
789     * @return The arc length of the specified segment of the curve.
790     */
791    @Override
792    public Parameter<Length> getArcLength(GeomPoint s1, GeomPoint s2, double eps) {
793        return getArcLength(s1.getValue(0), s2.getValue(0), eps);
794    }
795
796    /**
797     * Return a subrange point at the position on the curve with the specified arc-length.
798     * If the requested arc-length is &le; 0, the start of the curve is returned. If the
799     * requested arc length is &gt; the total arc length of the curve, then the end point
800     * is returned.
801     *
802     * @param targetArcLength The target arc length to find in meters.
803     * @param tol             Fractional tolerance (in parameter space) to refine the
804     *                        point position to.
805     * @return A subrange point on the curve at the specified arc-length position.
806     */
807    @Override
808    public SubrangePoint getPointAtArcLength(Parameter<Length> targetArcLength, double tol) {
809        //  Deal with an obvious case first.
810        if (targetArcLength.isApproxZero())
811            return SubrangePoint.newInstance(this, Point.valueOf(0));
812
813        double arcLengthValue = targetArcLength.getValue(SI.METER);
814
815        //  Calculate the full arc-length of the curve (to GTOL tolerance).
816        double fullArcLength = getArcLength(GTOL).getValue(SI.METER);
817
818        //  Is the requested arc length > than the actual arc length?
819        if (arcLengthValue > fullArcLength * (1 + GTOL))
820            return SubrangePoint.newInstance(this, Point.valueOf(1));
821
822        //  Calculate the arc-length tolerance to use based on the requested parameter
823        //  space tolerance and the full arc-length.
824        double eps = tol * fullArcLength / 100.;
825
826        return pointAtArcLength(0, arcLengthValue, tol, eps);
827    }
828
829    /**
830     * A point is found along this curve that when connected by a straight line to the
831     * given point in space, the resulting line is tangent to this curve. This method
832     * should only be used if this curve is a planar curve that is C1 (slope) continuous,
833     * and convex with respect to the point. The given point should be in the plane of the
834     * curve. If it is not, it will be projected into the plane of the curve.
835     *
836     * @param point The point that the tangent on the curve is to be found relative to.
837     * @param near  The parametric distance on the curve (0-1) that serves as an initial
838     *              guess at the location of the tangent point.
839     * @param tol   Fractional tolerance (in parameter space) to refine the point position
840     *              to.
841     * @return The point on this curve that is tangent to the curve and the supplied
842     *         point.
843     */
844    @Override
845    public SubrangePoint getTangencyPoint(GeomPoint point, double near, double tol) {
846        validateParameter(near);
847        requireNonNull(point, MessageFormat.format(RESOURCES.getString("paramNullErr"), "point"));
848
849        double minP = 0;
850        StackContext.enter();
851        try {
852            int cDim = this.getPhyDimension();
853            if (point.getPhyDimension() < cDim) {
854                point = point.toDimension(cDim);
855            }
856            if (point.getPhyDimension() != cDim) {
857                throw new DimensionException(RESOURCES.getString("crvPointDimensionErr"));
858            }
859
860            if (cDim > 2) {
861                //  Project the input point onto the plane of the curve.
862                Point p0 = this.getRealPoint(0);        //  Any point in plane (on curve).
863                Vector nhat = this.getBinormal(.5);     //  Plane normal vector.
864
865                //  Find the closest point on the plane.
866                Point projPoint = GeomUtil.pointPlaneClosest(point, p0, nhat);
867                point = projPoint;
868            }
869
870            //  Create a function that evaluates: 1 - |t(s) DOT ((q - p(s))/|q - p(s)|)|.
871            //  When this function is minimized, the s value may be a tangent point.
872            TangentPointEvaluatable tangentEvaluator = TangentPointEvaluatable.newInstance(this, point);
873            try {
874                //  Find the minimum of the evaluator function.
875                minP = Minimization.find(tangentEvaluator, 0, near, 1, tol, null);
876
877            } catch (RootException err) {
878                Logger.getLogger(AbstractCurve.class.getName()).log(Level.SEVERE,
879                        "Could not find tangency point.", err);
880            }
881            
882        } finally {
883            StackContext.exit();
884        }
885
886        SubrangePoint p = SubrangePoint.newInstance(this, Point.valueOf(minP));
887        return p;
888    }
889
890    /**
891     * Return a string of points that are gridded onto the curve using the specified
892     * spacing and gridding rule.
893     *
894     * @param gridRule Specifies whether the subdivision spacing is applied with respect
895     *                 to arc-length in real space or parameter space.
896     * @param spacing  A list of values used to define the subdivision spacing.
897     *                 <code>gridRule</code> specifies whether the subdivision spacing is
898     *                 applied with respect to arc-length in real space or parameter
899     *                 space. If the spacing is arc-length, then the values are
900     *                 interpreted as fractions of the arc length of the curve from 0 to
901     *                 1.
902     * @return A string of subrange points gridded onto the curve at the specified
903     *         parametric positions.
904     */
905    @Override
906    public PointString<SubrangePoint> extractGrid(GridRule gridRule, List<Double> spacing) {
907        requireNonNull(gridRule);
908        requireNonNull(spacing);
909        PointString<SubrangePoint> str = PointString.newInstance();
910
911        switch (gridRule) {
912            case PAR:
913                //  Parametric spacing.
914                for (double s : spacing) {
915                    if (s > 1)
916                        s = 1;
917                    else if (s < 0)
918                        s = 0;
919                    str.add(this.getPoint(s));
920                }
921                break;
922
923            case ARC:
924                //  Arc-length spacing.
925                double arcLength = this.getArcLength(GTOL).getValue(SI.METER);
926                double tol = GTOL / arcLength * 100.;   //  Tolerance for root finder in parameter space.
927
928                double sOld = 0;
929                for (double s : spacing) {
930                    double targetArcLength = s * arcLength;
931
932                    SubrangePoint p;
933                    if (targetArcLength >= arcLength)
934                        p = getPoint(1);
935                    else if (targetArcLength <= 0)
936                        p = getPoint(0);
937
938                    else {
939                        if (s < sOld)
940                            sOld = 0;
941                        p = pointAtArcLength(sOld, targetArcLength, tol, GTOL);
942                    }
943                    str.add(p);
944
945                    sOld = p.getParPosition().getValue(0);
946                }
947                break;
948
949            default:
950                throw new IllegalArgumentException(
951                        MessageFormat.format(RESOURCES.getString("crvUnsupportedGridRule"),
952                                gridRule.toString()));
953
954        }
955
956        return str;
957    }
958
959    /**
960     * Return a subrange point at the position on the curve with the specified arc length.
961     *
962     * @param x1              The minimum parametric position that could have the
963     *                        specified arc length.
964     * @param targetArcLength The target arc length to find in meters.
965     * @param tol             Fractional tolerance (in parameter space) to refine the
966     *                        point position to.
967     * @param eps             The desired fractional accuracy on the arc-length.
968     */
969    private SubrangePoint pointAtArcLength(double x1, double targetArcLength, double tol, double eps) {
970
971        double pos = 1;
972        StackContext.enter();
973        try {
974            //  Create a function that is used to iteratively find the target arc length.
975            ArcPosEvaluatable arcPosEvaluator = ArcPosEvaluatable.newInstance(this, targetArcLength, eps);
976
977            try {
978                //  Find the arc length position using the evaluator function.
979                arcPosEvaluator.firstPass = true;
980                pos = Roots.findRoot1D(arcPosEvaluator, x1, 1, tol);
981
982            } catch (RootException err) {
983                Logger.getLogger(AbstractCurve.class.getName()).log(Level.WARNING,
984                        "Could not find point at arc-length.", err);
985                pos = 1;
986            }
987        } finally {
988            StackContext.exit();
989        }
990        
991        SubrangePoint p = SubrangePoint.newInstance(this, Point.valueOf(pos));
992        return p;
993    }
994
995    private static final int MAX_POINTS = 10000;
996
997    /**
998     * Return a string of points that are gridded onto the curve with the number of points
999     * and spacing chosen to result in straight line segments between the points that have
1000     * mid-points that are all within the specified tolerance of this curve.
1001     *
1002     * @param tol The maximum distance that a straight line between gridded points may
1003     *            deviate from this curve. May not be null.
1004     * @return A string of subrange points gridded onto the curve.
1005     */
1006    @Override
1007    public PointString<SubrangePoint> gridToTolerance(Parameter<Length> tol) {
1008        if (isDegenerate(requireNonNull(tol))) {
1009            //  If this curve is a single point, just return that point twice.
1010            PointString<SubrangePoint> str = PointString.newInstance();
1011            str.add(getPoint(0));
1012            str.add(getPoint(1));
1013            return str;
1014        }
1015
1016        //  Start with the end-points.
1017        FastList<SubrangePoint> pntList = FastList.newInstance();
1018        pntList.add(getPoint(0));
1019        pntList.add(getPoint(1));
1020
1021        //  Refine the initial points using subdivision.
1022        Parameter<Area> tol2 = tol.times(tol);
1023        refineCurvePoints(pntList, tol2);
1024        
1025        //  Check guessed grid-to-tolerance points and add any that fall out of tolerance.
1026        FastList<SubrangePoint> guessPnts = guessGrid2TolPoints(tol2);
1027        checkGuessedGrid2TolPoints(pntList, guessPnts, tol2);
1028        FastList.recycle(guessPnts);
1029
1030        //  Refine the points after any guessed points were added using subdivision.
1031        refineCurvePoints(pntList, tol2);
1032        
1033        if (pntList.size() >= MAX_POINTS)
1034            Logger.getLogger(AbstractCurve.class.getName()).log(Level.WARNING,
1035                    MessageFormat.format(RESOURCES.getString("maxPointsWarningMsg"),
1036                    "curve", this.getID()));
1037
1038        //  Convert from a FastList to a PointString.
1039        PointString<SubrangePoint> str = PointString.valueOf(null, pntList);
1040        FastList.recycle(pntList);
1041        
1042        return str;
1043    }
1044
1045    /**
1046     * Use subdivision to refine the points gridded onto this curve until the
1047     * mid-points of all the segments meet the specified tolerance from the
1048     * curve. Points are added to the list wherever the segment mid-point
1049     * distances are greater than the tolerance.
1050     *
1051     * @param pntList The starting list of points to refine.
1052     * @param tol2 The square of the maximum distance that a straight line between
1053     *            gridded points may deviate from this curve. May not be null.
1054     */
1055    private void refineCurvePoints(FastList<SubrangePoint> pntList, Parameter<Area> tol2) {
1056        //  Start at the head of the list.
1057        FastList.Node<SubrangePoint> current = pntList.head().getNext();
1058        
1059        //  Loop until we reach the end of the list or exceed the max. number of points.
1060        FastList.Node tail = pntList.tail().getPrevious();
1061        while (!current.equals(tail) && pntList.size() < MAX_POINTS) {
1062            
1063            //  Subdivide the current segment until the mid-point is less than tol away from the curve.
1064            boolean distGTtol = true;
1065            while (distGTtol) {
1066                
1067                //  Get the current point.
1068                SubrangePoint pc = current.getValue();
1069                double sc = pc.getParPosition().getValue(0);
1070                
1071                //  Get the next point.
1072                FastList.Node<SubrangePoint> next = current.getNext();
1073                SubrangePoint pn = next.getValue();
1074                double sn = pn.getParPosition().getValue(0);
1075                
1076                //  Compute the average point half-way on a line from pc to pn.
1077                Point pa = GeomUtil.averagePoints(pc, pn);
1078                
1079                //  Find the parametric middle point on the curve.
1080                double sm = (sc + sn) / 2;
1081                SubrangePoint pm = getPoint(sm);
1082                
1083                //  Compute the distance between the middle point and the average point
1084                //  and compare with the tolerance.
1085                distGTtol = pa.distanceSq(pm).isGreaterThan(tol2);
1086                if (distGTtol)
1087                    //  Insert the new point into the list.
1088                    pntList.addBefore(next, pm);
1089                else
1090                    SubrangePoint.recycle(pm);
1091                
1092            }   //  end while()
1093            
1094            //  Move on to the next point.
1095            current = current.getNext();
1096
1097        }   //  end while()
1098    }
1099
1100    /**
1101     * A Comparator that compares the parametric S-position of a pair of
1102     * subrange points.
1103     */
1104    private static class SComparator implements Comparator<SubrangePoint> {
1105        
1106            @Override
1107            public int compare(SubrangePoint o1, SubrangePoint o2) {
1108                double s1 = o1.getParPosition().getValue(0);
1109                double s2 = o2.getParPosition().getValue(0);
1110                return Double.compare(s1, s2);
1111            }
1112    }
1113    
1114    private static final SComparator S_CMPRTR = new SComparator();
1115    
1116    /**
1117     * Check guessed grid-to-tolerance points and add them to pntList if they
1118     * are greater than the tolerance distance away from straight line segments
1119     * between the points already in pntList.
1120     * 
1121     * @param pntList The existing list of gridded points on this curve. Any
1122     * points from guessPnts that are more than tol away from straight line
1123     * segments between these points are added to this list.
1124     * @param guessPnts The existing list of segment parameter derived guessed points.
1125     * @param tol2 The square of the maximum distance that a straight line between
1126     *            gridded points may deviate from this curve. May not be null.
1127     */
1128    private void checkGuessedGrid2TolPoints(FastList<SubrangePoint> pntList, 
1129            FastList<SubrangePoint> guessPnts, Parameter<Area> tol2) {
1130        
1131        int nSegs = guessPnts.size() - 1;
1132        for (int i=1; i < nSegs; ++i) {
1133            SubrangePoint ps = guessPnts.get(i);
1134            double s = ps.getParPosition().getValue(0);
1135            
1136            //  Find where this point should fall in the gridded list of points.
1137            int idx = Collections.binarySearch(pntList, ps, S_CMPRTR);
1138            if (idx < 0) {
1139                //  ps is not already in the pntList, so see if it is in tolerance or not.
1140                idx = -idx - 1;
1141                
1142                //  Calculate a point, pa, along the line between p(i) and p(i-1)
1143                //  at the position of "s".
1144                SubrangePoint pi = pntList.get(idx);
1145                SubrangePoint pim1 = pntList.get(idx-1);
1146                Point pa = interpSPoint(pim1,pi,s);
1147                
1148                //  Compute the distance between the point at "s" and the line seg point
1149                //  and compare with the tolerance.
1150                boolean distGTtol = pa.distanceSq(ps).isGreaterThan(tol2);
1151                if (distGTtol) {
1152                    //  Insert the new point into the list.
1153                    pntList.add(idx, ps);
1154                    if (pntList.size() == MAX_POINTS)
1155                        break;
1156                } else
1157                    SubrangePoint.recycle(ps);
1158            } else
1159                SubrangePoint.recycle(ps);
1160        }
1161    }
1162
1163    /**
1164     * Interpolate a point at "s" between the points p1 and p2 (that must have
1165     * parametric positions surrounding "s").
1166     *
1167     * @param p1 Point with a parametric position less than "s".
1168     * @param p2 Point with a parametric position greater than "s".
1169     * @param s The parametric position at which a point is to be interpolated
1170     * between p1 and p2.
1171     * @return The interpolated point.
1172     */
1173    protected static Point interpSPoint(SubrangePoint p1, SubrangePoint p2, double s) {
1174        StackContext.enter();
1175        try {
1176            double s1 = p1.getParPosition().getValue(0);
1177            double s2 = p2.getParPosition().getValue(0);
1178            Point Ds = p2.minus(p1);
1179            Point pout = p1.plus(Ds.times((s - s1) / (s2 - s1)));
1180            return StackContext.outerCopy(pout);
1181        } finally {
1182            StackContext.exit();
1183        }
1184    }
1185    
1186    /**
1187     * Guess an initial set of points to use when gridding to tolerance. The initial
1188     * points must include the start and end of the curve as well as any discontinuities
1189     * in the curve. Other than that, they should be distributed as best as possible to
1190     * indicate areas of higher curvature, but should be generated as efficiently as
1191     * possible. The guessed points MUST be monotonically increasing from 0 to 1.
1192     * <p>
1193     * The default implementation places points at the positions returned by
1194     * "getSegmentParameters()". Sub-classes should override this method if they can
1195     * provide a more efficient/robust method to locate points where there may be
1196     * curvature discontinuities or areas of high curvature.
1197     * </p>
1198     *
1199     * @param tol2 The square of the maximum distance that a straight line
1200     *             between gridded points may deviate from this curve.
1201     * @return A FastList of subrange points to use as a starting point for the grid to
1202     *         tolerance algorithms.
1203     */
1204    protected FastList<SubrangePoint> guessGrid2TolPoints(Parameter<Area> tol2) {
1205
1206        //  Extract the boundaries of curve segments and place grid points at those.
1207        FastList<SubrangePoint> pntList = FastList.newInstance();
1208
1209        //  Create a list of segment starting parametric positions ( and the last end position, 1).
1210        FastTable<Double> bParams = getSegmentParameters();
1211        if (bParams.size() < 3) {
1212            //  We only have 0.0 & 1.0.  Insert 0.5 to ensure there is always at least one intermediate point.
1213            bParams.add(1, 0.5);
1214        }
1215
1216        //  Loop over the segment starting positions.
1217        int nSegs = bParams.size();
1218        for (int i = 0; i < nSegs; ++i) {
1219            double s0 = bParams.get(i);
1220            pntList.add(getPoint(s0));          //  Add the start-point of each segment.
1221        }
1222
1223        //  Cleanup before leaving.
1224        FastTable.recycle(bParams);
1225
1226        return pntList;
1227    }
1228
1229    /**
1230     * Return the intersection between a plane and this curve.
1231     *
1232     * @param plane The plane to intersect with this curve.
1233     * @param tol   Fractional tolerance (in parameter space) to refine the point
1234     *              positions to.
1235     * @return A PointString containing zero or more subrange points made by the
1236     *         intersection of this curve with the specified plane. If no intersection is
1237     *         found, an empty PointString is returned.
1238     */
1239    @Override
1240    public PointString<SubrangePoint> intersect(GeomPlane plane, double tol) {
1241
1242        //  Algorithm:  Finds places on the curve where the signed distance from a curve
1243        //  point to the plane is zero:  n dot (p(s) - pp) = 0.
1244        //  See:  http://mathworld.wolfram.com/Point-PlaneDistance.html
1245        
1246        //  Create a function that is used to iteratively find the intersections.
1247        int numDims = getPhyDimension();
1248        if (numDims < plane.getPhyDimension())
1249            numDims = plane.getPhyDimension();
1250        PlaneIntersectEvaluatable intEvaluator
1251                = PlaneIntersectEvaluatable.newInstance(this.toDimension(numDims), plane.toDimension(numDims));
1252
1253        PointString<SubrangePoint> output = PointString.newInstance();
1254        try {
1255
1256            //  Guess brackets surrounding the intersection points (if any).
1257            List<BracketRoot1D> brackets = guessIntersect(plane, intEvaluator);
1258
1259            //  Loop over all the brackets found.
1260            for (BracketRoot1D bracket : brackets) {
1261                //  Find the intersections using the evaluator function.
1262                double pos = Roots.findRoot1D(intEvaluator, bracket, tol);
1263                SubrangePoint p = SubrangePoint.newInstance(this, Point.valueOf(pos));
1264
1265                //  Add the point to the list of intersections being collected.
1266                output.add(p);
1267            }
1268
1269            //  Remove any near duplicates.
1270            int numPnts = output.size();
1271            for (int i = 0; i < numPnts; ++i) {
1272                SubrangePoint pi = output.get(i);
1273                for (int j = i + 1; j < numPnts; ++j) {
1274                    SubrangePoint pj = output.get(j);
1275                    double d = pj.getParPosition().distanceValue(pi.getParPosition());
1276                    if (d <= tol) {
1277                        output.remove(j);
1278                        --numPnts;
1279                        --j;
1280                    }
1281                }
1282            }
1283
1284        } catch (RootException err) {
1285            Logger.getLogger(AbstractCurve.class.getName()).log(Level.SEVERE,
1286                    "Could not intersect this curve witih a plane.", err);
1287
1288        } finally {
1289            PlaneIntersectEvaluatable.recycle(intEvaluator);
1290        }
1291
1292        return output;
1293    }
1294
1295    /**
1296     * Return a list of brackets that surround possible plane intersection points. Each of
1297     * the brackets will be refined using a root solver in order to find the actual
1298     * intersection points.
1299     * <p>
1300     * The default implementation sub-divides the curve into 17 equal-sized segments and
1301     * looks for a single bracket inside each segment. Sub-classes should override this
1302     * method if they can provide a more efficient/robust method to locate plane
1303     * intersection points.
1304     * </p>
1305     *
1306     * @param plane The plane to find the intersection points on this curve to/from.
1307     * @param eval  A function that calculates the <code>n dot (p(s) - pp)</code> on the
1308     *              curve.
1309     * @return A list of brackets that each surround a potential intersection point.
1310     *         Returns an empty list if there are no intersections found.
1311     */
1312    protected List<BracketRoot1D> guessIntersect(GeomPlane plane, Evaluatable1D eval) {
1313
1314        //  Find multiple intersection points (roots) by sub-dividing curve into pieces and
1315        //  bracketing any roots that occur in any of those segments.
1316        List<BracketRoot1D> brackets;
1317        try {
1318
1319            brackets = BracketRoot1D.findBrackets(requireNonNull(eval), 0, 1, getNumPointsPerSegment() * 2 + 1);
1320
1321        } catch (RootException err) {
1322            Logger.getLogger(AbstractCurve.class.getName()).log(Level.WARNING,
1323                    "guessIntersect(plane, eval) failed.", err);
1324            return FastTable.newInstance();
1325        }
1326
1327        return brackets;
1328    }
1329
1330    /**
1331     * Return the intersection between an infinite line and this curve.
1332     *
1333     * @param L0   A point on the line (origin of the line).
1334     * @param Ldir The direction vector for the line (does not have to be a unit vector).
1335     * @param tol  Tolerance (in physical space) to refine the point positions to.
1336     * @return A PointString containing zero or more subrange points made by the
1337     *         intersection of this curve with the specified line. If no intersection is
1338     *         found, an empty PointString is returned.
1339     */
1340    @Override
1341    public PointString<SubrangePoint> intersect(GeomPoint L0, GeomVector Ldir, Parameter<Length> tol) {
1342        requireNonNull(tol);
1343        
1344        //  First check to see if the line even comes near the curve.
1345        if (!GeomUtil.lineBoundsIntersect(L0, Ldir, this, null))
1346            return PointString.newInstance();
1347
1348        //  Make sure the line and curve have the same dimensions.
1349        int numDims = getPhyDimension();
1350        if (L0.getPhyDimension() > numDims)
1351            throw new IllegalArgumentException(MessageFormat.format(
1352                    RESOURCES.getString("sameOrFewerDimErr"), "line", "curve"));
1353        else if (L0.getPhyDimension() < numDims)
1354            L0 = L0.toDimension(numDims);
1355        if (Ldir.getPhyDimension() > numDims)
1356            throw new IllegalArgumentException(MessageFormat.format(
1357                    RESOURCES.getString("sameOrFewerDimErr"), "line", "curve"));
1358        else if (Ldir.getPhyDimension() < numDims)
1359            Ldir = Ldir.toDimension(numDims);
1360
1361        //  Create a line-curve intersector object.
1362        LineCurveIntersect intersector = LineCurveIntersect.newInstance(this, tol, L0, Ldir);
1363
1364        //  Intersect the line and curve (intersections are added to "output").
1365        PointString<SubrangePoint> output = intersector.intersect();
1366
1367        //  Clean up.
1368        LineCurveIntersect.recycle(intersector);
1369
1370        return output;
1371    }
1372
1373    /**
1374     * Return the intersection between a line segment and this curve.
1375     *
1376     * @param line The line segment to intersect.
1377     * @param tol  Tolerance (in physical space) to refine the point positions to.
1378     * @return A list containing two PointString instances each containing zero or more
1379     *         subrange points, on this curve and the input line segment respectively,
1380     *         made by the intersection of this curve with the specified line segment. If
1381     *         no intersections are found a list of two empty PointStrings are returned.
1382     */
1383    @Override
1384    public GeomList<PointString<SubrangePoint>> intersect(LineSegment line, Parameter<Length> tol) {
1385        requireNonNull(line);
1386        requireNonNull(tol);
1387        
1388        //  First check to see if the line segment even comes near the curve.
1389        GeomList<PointString<SubrangePoint>> output = GeomList.newInstance();
1390        output.add(PointString.newInstance());
1391        output.add(PointString.newInstance());
1392        if (!GeomUtil.lineSegBoundsIntersect(line.getStart(), line.getEnd(), this))
1393            return output;
1394
1395        //  Check for a degenerate curve case.
1396        if (isDegenerate(tol)) {
1397            //  Just see if the start point on this curve intersects that curve (is on that curve).
1398            SubrangePoint p1 = getPoint(0);
1399            Point p1r = p1.copyToReal();
1400
1401            //  Find the closest point on the surface.
1402            SubrangePoint p2 = line.getClosest(p1r, GTOL);
1403
1404            if (p1r.distance(p2).isLessThan(tol)) {
1405                //  Add the points to the lists of intersections being collected.
1406                output.get(0).add(p1);
1407                output.get(1).add(p2);
1408            }
1409
1410            return output;
1411        }
1412
1413        //  Make sure the line and curve have the same dimensions.
1414        int numDims = getPhyDimension();
1415        if (line.getPhyDimension() > numDims)
1416            throw new IllegalArgumentException(MessageFormat.format(
1417                    RESOURCES.getString("sameOrFewerDimErr"), "line", "curve"));
1418        else if (line.getPhyDimension() < numDims)
1419            line = line.toDimension(numDims);
1420
1421        //  Create a line segment-curve intersector object.
1422        LineSegCurveIntersect intersector = LineSegCurveIntersect.newInstance(this, tol, line, output);
1423
1424        //  Intersect the line segment and curve (intersections are added to "output").
1425        output = intersector.intersect();
1426
1427        //  Clean up.
1428        LineSegCurveIntersect.recycle(intersector);
1429
1430        return output;
1431    }
1432
1433    /**
1434     * Return the intersection between another curve and this curve.
1435     *
1436     * @param curve The curve to intersect with this curve.
1437     * @param tol   Tolerance (in physical space) to refine the point positions to.
1438     * @return A list containing two PointString instances each containing zero or more
1439     *         subrange points, on this curve and the input curve respectively, made by
1440     *         the intersection of this curve with the specified curve. If no
1441     *         intersections are found a list of two empty PointStrings are returned.
1442     */
1443    @Override
1444    public GeomList<PointString<SubrangePoint>> intersect(Curve curve, Parameter<Length> tol) {
1445        requireNonNull(curve);
1446        requireNonNull(tol);
1447        
1448        //  First check to see if the line segment even comes near the curve.
1449        GeomList<PointString<SubrangePoint>> output = GeomList.newInstance();
1450        output.add(PointString.newInstance());
1451        output.add(PointString.newInstance());
1452        if (!GeomUtil.boundsIntersect(curve, this))
1453            return output;
1454
1455        //  Check for a degenerate curve case.
1456        if (isDegenerate(tol)) {
1457            //  Just see if the start point on this curve intersects that curve (is on that curve).
1458            SubrangePoint p1 = getPoint(0);
1459            Point p1r = p1.copyToReal();
1460
1461            //  Find the closest point on the surface.
1462            SubrangePoint p2 = curve.getClosest(p1r, GTOL);
1463
1464            if (p1r.distance(p2).isLessThan(tol)) {
1465                //  Add the points to the lists of intersections being collected.
1466                output.get(0).add(p1);
1467                output.get(1).add(p2);
1468            }
1469
1470            return output;
1471        }
1472
1473        //  Create a curve-curve intersector object.
1474        CurveCurveIntersect intersector = CurveCurveIntersect.newInstance(this, tol, curve, output);
1475
1476        //  Intersect that curve and this curve (intersections are added to "output").
1477        output = intersector.intersect();
1478
1479        //  Clean up.
1480        CurveCurveIntersect.recycle(intersector);
1481
1482        return output;
1483    }
1484
1485    /**
1486     * Return the intersections between a surface and this curve.
1487     *
1488     * @param surface The surface to intersect with this curve.
1489     * @param tol     Tolerance (in physical space) to refine the point positions to.
1490     * @return A list containing two PointString instances each containing zero or more
1491     *         subrange points, on this curve and the input surface respectively, made by
1492     *         the intersection of this curve with the specified surface. If no
1493     *         intersections are found a list of two empty PointStrings are returned.
1494     */
1495    @Override
1496    public GeomList<PointString<SubrangePoint>> intersect(Surface surface, Parameter<Length> tol) {
1497        requireNonNull(surface);
1498        requireNonNull(tol);
1499        
1500        PointString<SubrangePoint> thisOutput = PointString.newInstance();
1501        PointString<SubrangePoint> thatOutput = PointString.newInstance();
1502        GeomList<PointString<SubrangePoint>> output = GeomList.newInstance();
1503        output.add(thisOutput);
1504        output.add(thatOutput);
1505
1506        //  First check to see if the curve even comes near the surface.
1507        if (!GeomUtil.boundsIntersect(this, surface))
1508            return output;
1509
1510        //  Make sure the line and curve have the same dimensions.
1511        Surface origSrf = surface;
1512        int numDims = getPhyDimension();
1513        if (surface.getPhyDimension() > numDims)
1514            throw new IllegalArgumentException(MessageFormat.format(
1515                    RESOURCES.getString("sameOrFewerDimErr"), "surface", "curve"));
1516        else if (surface.getPhyDimension() < numDims)
1517            surface = surface.toDimension(numDims);
1518
1519        //  Check for a degenerate curve case.
1520        if (isDegenerate(tol)) {
1521            //  Just see if the start point on the curve intersects the surface (is on the surface).
1522            SubrangePoint p1 = getPoint(0);
1523            Point p1r = p1.copyToReal();
1524
1525            //  Find the closest point on the surface.
1526            SubrangePoint p2 = origSrf.getClosest(p1r, GTOL);
1527
1528            if (p1r.distance(p2).isLessThan(tol)) {
1529                //  Add the points to the lists of intersections being collected.
1530                thisOutput.add(p1);
1531                thatOutput.add(p2);
1532            }
1533
1534            return output;
1535        }
1536
1537        //  Convert the inputs to the units of this surface.
1538        Unit<Length> unit = getUnit();
1539        surface = surface.to(unit);
1540        tol = tol.to(unit);
1541
1542        //  Determine the parameter tolerance to use.
1543        double arcLength = this.getArcLength(GTOL * 100).getValue();
1544        double ptol = tol.getValue() / arcLength;   //  Tolerance for root finder in parameter space.
1545        if (ptol < GTOL)
1546            ptol = GTOL;
1547
1548        //  Find potential intersection points (roots) by sub-dividing curve into pieces and
1549        //  bracketing any roots that occur in any of those segments.  Then, a 1D root finder
1550        //  is used to find the bracketed roots.
1551        CurveSrfIntersectFun eval = CurveSrfIntersectFun.newInstance(this, surface, GTOL * 1000);
1552        try {
1553
1554            //  Guess brackets surrounding the intersection points (if any).
1555            List<BracketRoot1D> brackets = guessIntersect(surface, eval);
1556
1557            //  Increase the tolerance for actual root finding.
1558            eval.tol = GTOL;
1559
1560            //  Loop over all the brackets found.
1561            for (BracketRoot1D bracket : brackets) {
1562
1563                //  Find the intersections using the evaluator function.
1564                eval.firstPass = true;
1565                double pos = Roots.findRoot1D(eval, bracket, ptol);
1566                SubrangePoint p1 = SubrangePoint.newInstance(this, Point.valueOf(pos));
1567                Point p1r = p1.copyToReal();
1568
1569                //  Find the closest point on the surface.
1570                SubrangePoint p2 = origSrf.getClosest(p1r, eval._nearS, eval._nearT, GTOL);
1571
1572                if (p1r.distance(p2).isLessThan(tol)) {
1573                    //  Add the points to the lists of intersections being collected.
1574                    thisOutput.add(p1);
1575                    thatOutput.add(p2);
1576                }
1577            }
1578
1579        } catch (RootException e) {
1580            //e.printStackTrace();
1581
1582        } finally {
1583            CurveSrfIntersectFun.recycle(eval);
1584        }
1585
1586        return output;
1587    }
1588
1589    /**
1590     * Return a list of brackets that surround possible surface intersection points. Each
1591     * of the brackets will be refined using a root solver in order to find the actual
1592     * intersection points.
1593     * <p>
1594     * The default implementation sub-divides the curve into 17 equal-sized segments and
1595     * looks for a single bracket inside each segment. Sub-classes should override this
1596     * method if they can provide a more efficient/robust method to locate surface
1597     * intersection points.
1598     * </p>
1599     *
1600     * @param surface The surface to find the intersection points on this curve to/from.
1601     * @param eval    A function that calculates the <code>n dot (p(s) - pp)</code> on the
1602     *                curve.
1603     * @return A list of brackets that each surround a potential intersection point.
1604     *         Returns an empty list if there are no intersections found.
1605     */
1606    protected List<BracketRoot1D> guessIntersect(Surface surface, Evaluatable1D eval) {
1607
1608        //  Find multiple intersection points (roots) by sub-dividing curve into pieces and
1609        //  bracketing any roots that occur in any of those segments.
1610        List<BracketRoot1D> brackets;
1611        try {
1612
1613            brackets = BracketRoot1D.findBrackets(requireNonNull(eval), 0, 1, getNumPointsPerSegment() * 2 + 1);
1614
1615        } catch (RootException err) {
1616            Logger.getLogger(AbstractCurve.class.getName()).log(Level.WARNING,
1617                    "guessIntersect(surface, eval) failed.", err);
1618            return FastTable.newInstance();
1619        }
1620
1621        return brackets;
1622    }
1623
1624    /**
1625     * Returns the closest point on this curve to the specified point.
1626     *
1627     * @param point The point to find the closest point on this curve to.
1628     * @param tol   Fractional tolerance (in parameter space) to refine the point position
1629     *              to.
1630     * @return The point on this curve that is closest to the specified point.
1631     */
1632    @Override
1633    public SubrangePoint getClosest(GeomPoint point, double tol) {
1634        return getClosestFarthest(point, tol, true);
1635    }
1636
1637    /**
1638     * Returns the closest point on this curve to the specified point near the specified
1639     * parametric position. This may not return the global closest point!
1640     *
1641     * @param point The point to find the closest point on this curve to.
1642     * @param near  The parametric position where the search for the closest point should
1643     *              be started.
1644     * @param tol   Fractional tolerance (in parameter space) to refine the point position
1645     *              to.
1646     * @return The point on this curve that is closest to the specified point near the
1647     *         specified parametric position.
1648     */
1649    @Override
1650    public SubrangePoint getClosest(GeomPoint point, double near, double tol) {
1651        return getClosestFarthestNear(point, near, tol, true);
1652    }
1653
1654    /**
1655     * Returns the closest points on this curve to the specified list of points.
1656     *
1657     * @param points The list of points to find the closest points on this curve to.
1658     * @param tol    Fractional tolerance (in parameter space) to refine the point
1659     *               positions to.
1660     * @return The list of {@link SubrangePoint} objects on this curve that are closest to
1661     *         the specified points.
1662     */
1663    @Override
1664    public PointString<SubrangePoint> getClosest(List<? extends GeomPoint> points, double tol) {
1665        return getClosestFarthest(points, tol, true);
1666    }
1667
1668    /**
1669     * Returns the farthest point on this curve from the specified point.
1670     *
1671     * @param point The point to find the farthest point on this curve from.
1672     * @param tol   Fractional tolerance (in parameter space) to refine the point position
1673     *              to.
1674     * @return The {@link SubrangePoint} on this curve that is farthest from the specified
1675     *         point.
1676     */
1677    @Override
1678    public SubrangePoint getFarthest(GeomPoint point, double tol) {
1679        return getClosestFarthest(point, tol, false);
1680    }
1681
1682    /**
1683     * Returns the farthest point on this curve from the specified point near the
1684     * specified parametric position. This may not return the global farthest point!
1685     *
1686     * @param point The point to find the farthest point on this curve from.
1687     * @param near  The parametric position where the search for the farthest point should
1688     *              be started.
1689     * @param tol   Fractional tolerance (in parameter space) to refine the point position
1690     *              to.
1691     * @return The {@link SubrangePoint} on this curve that is farthest from the specified
1692     *         point.
1693     */
1694    @Override
1695    public SubrangePoint getFarthest(GeomPoint point, double near, double tol) {
1696        return getClosestFarthestNear(point, near, tol, false);
1697    }
1698
1699    /**
1700     * Returns the farthest points on this curve to the specified list of points.
1701     *
1702     * @param points The list of points to find the farthest points on this curve to.
1703     * @param tol    Fractional tolerance (in parameter space) to refine the point
1704     *               positions to.
1705     * @return The list of {@link SubrangePoint} objects on this curve that are farthest
1706     *         from the specified points.
1707     */
1708    @Override
1709    public PointString<SubrangePoint> getFarthest(List<? extends GeomPoint> points, double tol) {
1710        return getClosestFarthest(points, tol, false);
1711    }
1712
1713    /**
1714     * Returns the closest/farthest points on this curve to the specified point.
1715     *
1716     * @param point   The point to find the closest/farthest point on this curve from.
1717     * @param tol     Fractional tolerance (in parameter space) to refine the point
1718     *                position to.
1719     * @param closest Set to <code>true</code> to return the closest point or
1720     *                <code>false</code> to return the farthest point.
1721     * @return The {@link SubrangePoint} on this curve that is closest/farthest from the
1722     *         specified point.
1723     */
1724    protected SubrangePoint getClosestFarthest(GeomPoint point, double tol, boolean closest) {
1725
1726        double sout = 0;
1727        StackContext.enter();
1728        try {
1729            //  Make sure the point and curve have the same dimensions.
1730            int numDims = getPhyDimension();
1731            if (point.getPhyDimension() > numDims) {
1732                throw new IllegalArgumentException(MessageFormat.format(
1733                        RESOURCES.getString("sameOrFewerDimErr"), "point", "curve"));
1734            } else if (point.getPhyDimension() < numDims) {
1735                point = point.toDimension(numDims);
1736            }
1737
1738            //  Create a function that calculates the dot product of (p-q) and ps.  When this
1739            //  is zero, a point on the curve that forms a vector perpendicular to the curve has been found.
1740            PerpPointEvaluatable perpPointEval = PerpPointEvaluatable.newInstance(this, point);
1741
1742            //  Find brackets surrounding possible closest/farthest points.
1743            List<BracketRoot1D> brackets = guessClosestFarthest(point, closest, perpPointEval);
1744
1745            //  Find the closest/farthest point.
1746            SubrangePoint output = getClosestFarthest(point, tol, closest, perpPointEval, brackets);
1747            sout = output.getParPosition().getValue(0);
1748
1749        } finally {
1750            StackContext.exit();
1751        }
1752
1753        return getPoint(sout);
1754    }
1755
1756    /**
1757     * Returns the closest/farthest points on this curve to the specified list of points.
1758     *
1759     * @param point   The point to find the closest/farthest point on this curve from.
1760     * @param near    The parametric position where the search for the farthest point
1761     *                should be started.
1762     * @param tol     Fractional tolerance (in parameter space) to refine the point
1763     *                position to.
1764     * @param closest Set to <code>true</code> to return the closest point or
1765     *                <code>false</code> to return the farthest point.
1766     * @return The {@link SubrangePoint} on this curve that is closest/farthest from the
1767     *         specified point.
1768     */
1769    private SubrangePoint getClosestFarthestNear(GeomPoint point, double near, double tol, boolean closest) {
1770        validateParameter(near);
1771        
1772        double sout = 0;
1773        StackContext.enter();
1774        try {
1775            //  Make sure the point and curve have the same dimensions.
1776            int numDims = getPhyDimension();
1777            if (point.getPhyDimension() > numDims) {
1778                throw new IllegalArgumentException(MessageFormat.format(
1779                        RESOURCES.getString("sameOrFewerDimErr"), "point", "curve"));
1780            } else if (point.getPhyDimension() < numDims) {
1781                point = point.toDimension(numDims);
1782            }
1783
1784            //  Create a function that calculates the dot product of (p-q) and ps.  When this
1785            //  is zero, a point on the curve that forms a vector perpendicular to the curve has been found.
1786            PerpPointEvaluatable perpPointEval = PerpPointEvaluatable.newInstance(this, point);
1787
1788            //  Find the closet/farthest point near the specified parametric position.
1789            try {
1790                List<BracketRoot1D> brackets = bracketNear(perpPointEval, near);
1791                SubrangePoint output = getClosestFarthest(point, tol, closest, perpPointEval, brackets);
1792                sout = output.getParPosition().getValue(0);
1793                
1794            } catch (RootException err) {
1795                Logger.getLogger(AbstractCurve.class.getName()).log(Level.WARNING,
1796                        "Failed to find closest/farthest point near existing point.", err);
1797                sout = 0;
1798            }
1799            
1800        } finally {
1801            StackContext.exit();
1802        }
1803        
1804        return getPoint(sout);
1805    }
1806
1807    /**
1808     * Returns the closest/farthest points on this curve to the specified list of points.
1809     *
1810     * @param points  The list of points to find the closest/farthest points on this curve
1811     *                to.
1812     * @param tol     Fractional tolerance (in parameter space) to refine the point
1813     *                positions to.
1814     * @param closest Set to <code>true</code> to return the closest point or
1815     *                <code>false</code> to return the farthest point.
1816     * @return The list of {@link SubrangePoint} objects on this curve that are
1817     *         closest/farthest to the specified points.
1818     */
1819    @SuppressWarnings("null")
1820    private PointString<SubrangePoint> getClosestFarthest(List<? extends GeomPoint> points, double tol, boolean closest) {
1821        PointString<SubrangePoint> output = PointString.newInstance();
1822
1823        //  Loop over all the target points in the supplied list and find the closest curve point for each one.
1824        SubrangePoint p_old = null;
1825        int numDims = getPhyDimension();
1826        int size = points.size();
1827        for (int i = 0; i < size; ++i) {
1828            //  Get the target point.
1829            GeomPoint q = points.get(i);
1830
1831            //  Make sure the point and curve have the same dimensions.
1832            if (q.getPhyDimension() > numDims)
1833                throw new IllegalArgumentException(MessageFormat.format(
1834                        RESOURCES.getString("sameOrFewerDimErr"), "points", "curve"));
1835            else if (q.getPhyDimension() < numDims)
1836                q = q.toDimension(numDims);
1837
1838            //  Create a function that calculates the dot product of (p(s)-q) and ps(s).  When this
1839            //  is zero, a point on the curve that forms a vector with the target perpendicular to the curve has been found.
1840            PerpPointEvaluatable perpPointEval = PerpPointEvaluatable.newInstance(this, q);
1841
1842            //  Find multiple closest/farthest points (roots) by sub-dividing curve into pieces and
1843            //  bracketing any roots that occur in any of those segments.
1844            List<BracketRoot1D> brackets;
1845            try {
1846                if (i == 0)
1847                    //  Look across the entire curve.
1848                    brackets = guessClosestFarthest(q, closest, perpPointEval);
1849                else
1850                    //  Look for brackets near the previous point.
1851                    brackets = bracketNear(perpPointEval, p_old.getParPosition().getValue(0));
1852            } catch (RootException err) {
1853                Logger.getLogger(AbstractCurve.class.getName()).log(Level.WARNING,
1854                        "Could not find closest/farthest point.", err);
1855                return output;
1856            }
1857
1858            //  Now use a root finder to actually find the closest/farthest point on the curve.
1859            SubrangePoint p = getClosestFarthest(q, tol, closest, perpPointEval, brackets);
1860            output.add(p);
1861            p_old = p;
1862
1863            //  Recycle temporary objects.
1864            PerpPointEvaluatable.recycle(perpPointEval);
1865        }
1866
1867        return output;
1868    }
1869
1870    /**
1871     * Return a list of brackets that surround possible closest/farthest points. Each of
1872     * the brackets will be refined using a root solver in order to find the actual
1873     * closest or farthest points.
1874     * <p>
1875     * The default implementation sub-divides the curve into 21 equal-sized segments and
1876     * looks for a single bracket inside each segment. Sub-classes should override this
1877     * method if they can provide a more efficient/robust method to locate
1878     * closest/farthest points.
1879     * </p>
1880     *
1881     * @param point   The point to find the closest/farthest point on this curve to/from
1882     *                (guaranteed to be the same physical dimension before this method is
1883     *                called).
1884     * @param closest Set to <code>true</code> to return the closest point or
1885     *                <code>false</code> to return the farthest point.
1886     * @param eval    A function that calculates the dot product of (p(s)-q) and ps(s) on
1887     *                the curve.
1888     * @return A list of brackets that each surround a potential closest/farthest point.
1889     */
1890    protected List<BracketRoot1D> guessClosestFarthest(GeomPoint point, boolean closest,
1891            Evaluatable1D eval) {
1892
1893        //  Find multiple closest/farthest points (roots) by sub-dividing curve into pieces and
1894        //  bracketing any roots that occur in any of those segments.
1895        List<BracketRoot1D> brackets;
1896        try {
1897
1898            brackets = BracketRoot1D.findBrackets(requireNonNull(eval), 0, 1, 21);
1899
1900        } catch (RootException err) {
1901            Logger.getLogger(AbstractCurve.class.getName()).log(Level.WARNING,
1902                    "Failed to guess closest/farthest point.", err);
1903            return FastTable.newInstance();
1904        }
1905
1906        return brackets;
1907    }
1908
1909    /**
1910     * Look for and bracket roots at parametric positions near the specified parametric
1911     * position value.
1912     *
1913     * @param eval A function that calculates the dot product of (p(s)-q) and ps(s) on the
1914     *             curve.
1915     * @param near The parametric position on the curve to begin search for roots near.
1916     * @return A list of brackets of roots to the perpPointEval function each of which
1917     *         could be the closest/farthest.
1918     * @throws jahuwaldt.tools.math.RootException if a bracket could not be found.
1919     */
1920    protected List<BracketRoot1D> bracketNear(Evaluatable1D eval, double near) throws RootException {
1921        requireNonNull(eval);
1922        List<BracketRoot1D> brackets = null;
1923
1924        //  Look near the previous parametric position by expanding the upper and lower limits until
1925        //  at least one bracket is found.
1926        double old_sm = near;
1927        double old_sp = near - 0.01;
1928        if (old_sp < 0)
1929            old_sp = 0;
1930        do {
1931            double sm = old_sm - 0.1;
1932            if (sm < 0)
1933                sm = 0;
1934            double sp = old_sp + 0.1;
1935            if (sp > 1)
1936                sp = 1;
1937            //  Look on either side of near.
1938            brackets = BracketRoot1D.findBrackets(eval, sm, sp, 10);
1939            old_sm = sm;
1940            old_sp = sp;
1941        } while (brackets.size() < 1 && !(MathTools.isApproxEqual(old_sp, 1) && MathTools.isApproxZero(old_sm)));
1942
1943        return brackets;
1944    }
1945
1946    /**
1947     * Returns the closest point on this curve, p, to the specified point, q.
1948     *
1949     * @param point         The point, q, to find the closest point on this curve to
1950     *                      (guaranteed to be the same physical dimension before this
1951     *                      method is called).
1952     * @param tol           Fractional tolerance (in parameter space) to refine the point
1953     *                      position to.
1954     * @param closest       Set to <code>true</code> to return the closest point or
1955     *                      <code>false</code> to return the farthest point.
1956     * @param perpPointEval A function that calculates the dot product of (p(s)-q) and
1957     *                      ps(s) on the curve.
1958     * @param brackets      A list of brackets of roots to the perpPointEval function each
1959     *                      of which could be the closest/farthest.
1960     * @return The {@link SubrangePoint} on this curve that is closest to the specified
1961     *         point.
1962     */
1963    private SubrangePoint getClosestFarthest(GeomPoint point, double tol, boolean closest, PerpPointEvaluatable perpPointEval, List<BracketRoot1D> brackets) {
1964        
1965        double sout = 0;
1966        StackContext.enter();
1967        try {
1968            //  Create a list of potential closest points and add the ends of the curve to it.
1969            FastTable<SubrangePoint> potentials = FastTable.newInstance();
1970            potentials.add(SubrangePoint.newInstance(this, 0));
1971            potentials.add(SubrangePoint.newInstance(this, 1));
1972
1973            try {
1974                //  Loop over all the brackets found.
1975                for (BracketRoot1D bracket : brackets) {
1976                    //  Find the perpendicular point on the curve using the evaluator function.
1977                    perpPointEval.firstPass = true;
1978                    double s = Roots.findRoot1D(perpPointEval, bracket, tol);
1979
1980                    //  Add the point to the list of potential closest/farthest points being collected.
1981                    potentials.add(SubrangePoint.newInstance(this, s));
1982                }
1983
1984            } catch (RootException err) {
1985                Logger.getLogger(AbstractCurve.class.getName()).log(Level.WARNING,
1986                        "Failed to find perpendicular point along curve.", err);
1987
1988            } finally {
1989                if (potentials.isEmpty()) {
1990                    sout = 0;
1991
1992                } else {
1993                    //  Loop over the potentials and find the one that is the closest/farthest.
1994                    SubrangePoint output = potentials.get(0);
1995                    double d2Opt = point.distanceSqValue(output);
1996                    int size = potentials.size();
1997                    for (int i = 1; i < size; ++i) {
1998                        SubrangePoint p = potentials.get(i);
1999                        double dist2 = point.distanceSqValue(p);
2000                        if (closest) {
2001                            if (dist2 < d2Opt) {
2002                                d2Opt = dist2;
2003                                output = p;
2004                            }
2005                        } else {
2006                            if (dist2 > d2Opt) {
2007                                d2Opt = dist2;
2008                                output = p;
2009                            }
2010                        }
2011                    }
2012                    sout = output.getParPosition().getValue(0);
2013                }
2014            }
2015        } finally {
2016            StackContext.exit();
2017        }
2018       
2019        return getPoint(sout);
2020    }
2021
2022    /**
2023     * Returns the closest points (giving the minimum distance) between this curve and the
2024     * specified curve. If the specified curve intersects this curve, then the 1st
2025     * intersection found is returned (not all the intersections are found and there is no
2026     * guarantee about which intersection will be returned). If the specified curve is
2027     * parallel to this curve (all points are equidistant away), then any point between
2028     * the two curves could be returned as the "closest".
2029     *
2030     * @param curve The curve to find the closest points between this curve and that one.
2031     * @param tol   Fractional tolerance (in parameter space) to refine the point position
2032     *              to.
2033     * @return A list containing two SubrangePoint objects that represent the closest
2034     *         point on this curve (index 0) and the closest point on the specified curve
2035     *         (index 1).
2036     */
2037    @Override
2038    public PointString<SubrangePoint> getClosest(Curve curve, double tol) {
2039
2040        //  Make sure that curve and this curve have the same dimensions.
2041        int numDims = getPhyDimension();
2042        if (curve.getPhyDimension() > numDims)
2043            throw new IllegalArgumentException(MessageFormat.format(
2044                    RESOURCES.getString("sameOrFewerDimErr"), "input curve", "curve"));
2045        else if (curve.getPhyDimension() < numDims)
2046            curve = curve.toDimension(numDims);
2047
2048        //  Create a function that calculates the dot product of (p-q) and ps where "q" is the closest
2049        //  point on the target curve.  When this is zero, a point on the curve that forms a
2050        //  vector perpendicular to the curve has been found.
2051        PerpPointParGeomEvaluatable perpPointEval = PerpPointParGeomEvaluatable.newInstance(this, curve, tol);
2052
2053        //  Find brackets surrounding possible closest/farthest points.
2054        List<BracketRoot1D> brackets = guessClosest(curve, perpPointEval);
2055
2056        //  Find the closest point.
2057        PointString<SubrangePoint> output = getClosest(curve, tol, perpPointEval, brackets);
2058
2059        //  Recycle temporary objects.
2060        PerpPointParGeomEvaluatable.recycle(perpPointEval);
2061
2062        return output;
2063    }
2064
2065    /**
2066     * Return a list of brackets that surround possible points on this curve that are
2067     * closest to the target parametric geometry object. Each of the brackets will be
2068     * refined using a root solver in order to find the actual closest points.
2069      * <p>
2070     * The default implementation sub-divides the curve into 21 equal-sized segments and
2071     * looks for a single bracket inside each segment. Sub-classes should override this
2072     * method if they can provide a more efficient/robust method to locate closest
2073     * points.
2074     * </p>
2075     *
2076     * @param parGeom The target parametric geometry to find the closest point on this
2077     *                curve to (guaranteed to be the same physical dimension before this
2078     *                method is called).
2079     * @param eval    A function that calculates the dot product of (p(s)-q) and ps(s) on
2080     *                the curve where "q" is the closest point on the target parametric
2081     *                geometry object.
2082     * @return A list of brackets that each surround a potential closest point on this
2083     *         curve.
2084     */
2085    protected List<BracketRoot1D> guessClosest(ParametricGeometry parGeom, PerpPointParGeomEvaluatable eval) {
2086
2087        //  Find multiple closest/farthest points (roots) by sub-dividing curve into pieces and
2088        //  bracketing any roots that occur in any of those segments.
2089        List<BracketRoot1D> brackets;
2090        try {
2091
2092            brackets = BracketRoot1D.findBrackets(requireNonNull(eval), 0, 1, 21);
2093
2094        } catch (RootException err) {
2095            Logger.getLogger(AbstractCurve.class.getName()).log(Level.WARNING,
2096                    "Failed to bracket closest points.", err);
2097            return FastTable.newInstance();
2098        }
2099
2100        return brackets;
2101    }
2102
2103    /**
2104     * Returns a list containing the closest point on this curve, p, to the closest point
2105     * on the target parametric geometry, q.
2106     *
2107     * @param parGeom       The parametric geometry to find the closest point on this
2108     *                      curve to (guaranteed to be the same physical dimension before
2109     *                      this method is called).
2110     * @param tol           Fractional tolerance (in parameter space) to refine the point
2111     *                      position to.
2112     * @param perpPointEval A function that calculates the dot product of (p(s)-q) and
2113     *                      ps(s) on the curve where "q" is the closest point on the
2114     *                      target parametric geometry object.
2115     * @param brackets      A list of brackets of roots to the perpPointEval function each
2116     *                      of which could be the closest/farthest.
2117     * @return A list containing SubrangePoint objects representing the closest points on
2118     *         this curve (index 0) and the target parametric geometry (index 1).
2119     */
2120    private PointString<SubrangePoint> getClosest(ParametricGeometry parGeom, double tol, PerpPointParGeomEvaluatable perpPointEval,
2121            List<BracketRoot1D> brackets) {
2122        GeomPoint thisParPos, thatParPos;
2123
2124        StackContext.enter();
2125        try {
2126            SubrangePoint thisOutput, thatOutput;
2127            
2128            //  Create a list of potential closest points and add the ends of this curve to it.
2129            FastTable<SubrangePoint> potentials = FastTable.newInstance();
2130            potentials.add(SubrangePoint.newInstance(this, 0));
2131            potentials.add(SubrangePoint.newInstance(this, 1));
2132
2133            try {
2134                //  Loop over all the brackets found.
2135                for (BracketRoot1D bracket : brackets) {
2136                    //  Find the perpendicular point on this curve using the evaluator function.
2137                    double s = Roots.findRoot1D(perpPointEval, bracket, tol);
2138
2139                    //  Add the point to the list of potential closest/farthest points being collected.
2140                    potentials.add(SubrangePoint.newInstance(this, s));
2141                }
2142
2143            } catch (RootException err) {
2144                Logger.getLogger(AbstractCurve.class.getName()).log(Level.WARNING,
2145                        "Failed to find perpendicular points on curve.", err);
2146
2147            } finally {
2148                if (potentials.isEmpty()) {
2149                    thisOutput = getPoint(0);
2150                    thatOutput = parGeom.getClosest(thisOutput, tol);
2151
2152                } else {
2153                    //  Loop over the potentials and find the one that is the closest.
2154                    thisOutput = potentials.get(0);
2155                    thatOutput = parGeom.getClosest(thisOutput, tol);
2156                    double dOpt = thatOutput.distanceValue(thisOutput);
2157                    int size = potentials.size();
2158                    for (int i = 1; i < size; ++i) {
2159                        SubrangePoint p = potentials.get(i);
2160                        SubrangePoint q = parGeom.getClosest(p, tol);
2161                        double dist = q.distanceValue(p);
2162                        if (dist < dOpt) {
2163                            dOpt = dist;
2164                            thisOutput = p;
2165                            thatOutput = q;
2166                        }
2167                    }
2168                }
2169            }
2170            
2171            //  Copy out the parametric positions of the closest points.
2172            thisParPos = StackContext.outerCopy(thisOutput.getParPosition().immutable());
2173            thatParPos = StackContext.outerCopy(thatOutput.getParPosition().immutable());
2174            
2175        } finally {
2176            StackContext.exit();
2177        }
2178
2179        //  Turn parametric positions back into actual subrange points.
2180        PointString<SubrangePoint> output = PointString.newInstance();
2181        output.add(SubrangePoint.newInstance(this, thisParPos));
2182        output.add(SubrangePoint.newInstance(parGeom, thatParPos));
2183
2184        return output;
2185    }
2186
2187    /**
2188     * Returns the closest points (giving the minimum distance) between this curve and the
2189     * specified surface. If this curve intersects the specified surface, then the 1st
2190     * intersection found is returned (not all the intersections are found and there is no
2191     * guarantee about which intersection will be returned). If this curve is parallel to
2192     * the specified surface (all closest points are equidistant away), then any point on
2193     * this curve and the closest point on the target surface will be returned as
2194     * "closest".
2195     *
2196     * @param surface The surface to find the closest points between this curve and that
2197     *                surface.
2198     * @param tol     Fractional tolerance (in parameter space) to refine the point
2199     *                position to.
2200     * @return A list containing two SubrangePoint objects that represent the closest
2201     *         point on this curve (index 0) and the closest point on the specified
2202     *         surface (index 1).
2203     */
2204    @Override
2205    public PointString<SubrangePoint> getClosest(Surface surface, double tol) {
2206
2207        //  Make sure that surface and this curve have the same dimensions.
2208        int numDims = getPhyDimension();
2209        if (surface.getPhyDimension() > numDims)
2210            throw new IllegalArgumentException(MessageFormat.format(
2211                    RESOURCES.getString("sameOrFewerDimErr"), "surface", "curve"));
2212        else if (surface.getPhyDimension() < numDims)
2213            surface = surface.toDimension(numDims);
2214
2215        //  Create a function that calculates the dot product of (p-q) and ps where "q" is the closest
2216        //  point on the target surface.  When this is zero, a point on the curve that forms a
2217        //  vector perpendicular to the curve has been found.
2218        PerpPointParGeomEvaluatable perpPointEval = PerpPointParGeomEvaluatable.newInstance(this, surface, tol);
2219
2220        //  Find brackets surrounding possible closest/farthest points.
2221        List<BracketRoot1D> brackets = guessClosest(surface, perpPointEval);
2222
2223        //  Find the closest point.
2224        PointString<SubrangePoint> output = getClosest(surface, tol, perpPointEval, brackets);
2225
2226        //  Recycle temporary objects.
2227        PerpPointParGeomEvaluatable.recycle(perpPointEval);
2228
2229        return output;
2230    }
2231
2232    /**
2233     * Returns the closest points (giving the minimum distance) between this curve and the
2234     * specified plane. If this curve intersects the specified plane, then the 1st
2235     * intersection found is returned (not all the intersections are found and there is no
2236     * guarantee about which intersection will be returned). If this curve is parallel to
2237     * the specified plane (all closest points are equidistant away), then any point on
2238     * this curve and the closest point on the target plane will be returned as "closest".
2239     *
2240     * @param plane The plane to find the closest points between this curve and that
2241     *              plane.
2242     * @param tol   Fractional tolerance (in parameter space) to refine the point position
2243     *              to.
2244     * @return The closest point found on this curve to the specified plane.
2245     */
2246    @Override
2247    public SubrangePoint getClosest(GeomPlane plane, double tol) {
2248
2249        //  Make sure that the input plane and this curve have the same dimensions.
2250        int numDims = getPhyDimension();
2251        if (plane.getPhyDimension() > numDims)
2252            throw new IllegalArgumentException(MessageFormat.format(
2253                    RESOURCES.getString("sameOrFewerDimErr"), "plane", "curve"));
2254        else if (plane.getPhyDimension() < numDims)
2255            plane = plane.toDimension(numDims);
2256
2257        //  Create a function that calculates the dot product of (p-q) and ps where "q" is the closest
2258        //  point on the target plane.  When this is zero, a point on the curve that forms a
2259        //  vector perpendicular to the curve has been found.
2260        PerpPointPlaneEvaluatable perpPointEval = PerpPointPlaneEvaluatable.newInstance(this, plane);
2261
2262        //  Find brackets surrounding possible closest/farthest points.
2263        List<BracketRoot1D> brackets = guessClosest(perpPointEval);
2264
2265        //  Find the closest point.
2266        SubrangePoint output = getClosest(plane, tol, perpPointEval, brackets);
2267
2268        //  Recycle temporary objects.
2269        PerpPointPlaneEvaluatable.recycle(perpPointEval);
2270
2271        return output;
2272    }
2273
2274    /**
2275     * Return a list of brackets that surround possible points on this curve that are
2276     * closest to the target geometry object. Each of the brackets will be refined using a
2277     * root solver in order to find the actual closest points.
2278     * <p>
2279     * The default implementation sub-divides the curve into 21 equal-sized segments and
2280     * looks for a single bracket inside each segment. Sub-classes should override this
2281     * method if they can provide a more efficient/robust method to locate closest
2282     * points.
2283     * </p>
2284     *
2285     * @param eval A function that calculates the dot product of (p(s)-q) and ps(s) on the
2286     *             curve where "q" is the closest point on the target geometry object.
2287     * @return A list of brackets that each surround a potential closest point on this
2288     *         curve.
2289     */
2290    protected List<BracketRoot1D> guessClosest(PerpPointPlaneEvaluatable eval) {
2291
2292        //  Find multiple closest/farthest points (roots) by sub-dividing curve into pieces and
2293        //  bracketing any roots that occur in any of those segments.
2294        List<BracketRoot1D> brackets;
2295        try {
2296            brackets = BracketRoot1D.findBrackets(requireNonNull(eval), 0, 1, 21);
2297
2298        } catch (RootException err) {
2299            Logger.getLogger(AbstractCurve.class.getName()).log(Level.WARNING,
2300                    "Failed to bracket closest points", err);
2301            return FastTable.newInstance();
2302        }
2303
2304        return brackets;
2305    }
2306
2307    /**
2308     * Returns a list containing the closest point on this curve, p, to the closest point
2309     * on the target plane, q.
2310     *
2311     * @param plane         The plane to find the closest point on this curve to
2312     *                      (guaranteed to be the same physical dimension before this
2313     *                      method is called).
2314     * @param tol           Fractional tolerance (in parameter space) to refine the point
2315     *                      position to.
2316     * @param perpPointEval A function that calculates the dot product of (p(s)-q) and
2317     *                      ps(s) on the curve where "q" is the closest point on the
2318     *                      target plane object.
2319     * @param brackets      A list of brackets of roots to the perpPointEval function each
2320     *                      of which could be the closest/farthest.
2321     * @return The closest point found on this curve to the specified plane.
2322     */
2323    private SubrangePoint getClosest(GeomPlane plane, double tol, PerpPointPlaneEvaluatable perpPointEval,
2324            List<BracketRoot1D> brackets) {
2325        double sOutput;
2326        
2327        StackContext.enter();
2328        try {
2329            SubrangePoint thisOutput;
2330
2331            //  Create a list of potential closest points and add the ends of this curve to it.
2332            FastTable<SubrangePoint> potentials = FastTable.newInstance();
2333            potentials.add(SubrangePoint.newInstance(this, 0));
2334            potentials.add(SubrangePoint.newInstance(this, 1));
2335
2336            //  Loop over all the brackets found.
2337            for (BracketRoot1D bracket : brackets) {
2338                try {
2339                    //  Find the perpendicular point on this curve using the evaluator function.
2340                    double s = Roots.findRoot1D(perpPointEval, bracket, tol);
2341
2342                    //  Add the point to the list of potential closest/farthest points being collected.
2343                    potentials.add(SubrangePoint.newInstance(this, s));
2344
2345                } catch (RootException e) {
2346                    //  Just ignore.
2347                    //e.printStackTrace();
2348                }
2349            }
2350
2351            if (potentials.isEmpty()) {
2352                thisOutput = getPoint(0);
2353
2354            } else {
2355                //  Loop over the potentials and find the one that is the closest.
2356                thisOutput = potentials.get(0);
2357                GeomPoint q = plane.getClosest(thisOutput);
2358                double dOpt = q.distanceValue(thisOutput);
2359                int size = potentials.size();
2360                for (int i = 1; i < size; ++i) {
2361                    SubrangePoint p = potentials.get(i);
2362                    q = plane.getClosest(p);
2363                    double dist = q.distanceValue(p);
2364                    if (dist < dOpt) {
2365                        dOpt = dist;
2366                        thisOutput = p;
2367                    }
2368                }
2369            }
2370
2371            //  Get the parametric position on this curve of the closest point.
2372            sOutput = thisOutput.getParPosition().getValue(0);
2373
2374        } finally {
2375            StackContext.exit();
2376        }
2377        
2378        return getPoint(sOutput);
2379    }
2380
2381    /**
2382     * Returns the most extreme point, either minimum or maximum, in the specified
2383     * coordinate direction on this curve. This is more accurate than
2384     * <code>getBoundsMax</code> &amp; <code>getBoundsMin</code>, but also typically takes
2385     * longer to compute.
2386     *
2387     * @param dim An index indicating the dimension to find the min/max point for (0=X,
2388     *            1=Y, 2=Z, etc).
2389     * @param max Set to <code>true</code> to return the maximum value, <code>false</code>
2390     *            to return the minimum.
2391     * @param tol Fractional tolerance (in parameter space) to refine the min/max point
2392     *            position to.
2393     * @return The point found on this curve that is the min or max in the specified
2394     *         coordinate direction.
2395     * @see #getBoundsMin
2396     * @see #getBoundsMax
2397     */
2398    @Override
2399    public SubrangePoint getLimitPoint(int dim, boolean max, double tol) {
2400        //  Check the input dimension for consistancy.
2401        int thisDim = getPhyDimension();
2402        if (dim < 0 || dim >= thisDim)
2403            throw new DimensionException(MessageFormat.format(
2404                    RESOURCES.getString("inputDimOutOfRange"), "curve"));
2405
2406        double sout = 0;
2407        StackContext.enter();
2408        try {
2409            //  Get the bounds of the curve.
2410            Point boundsMax = getBoundsMax();
2411            Point boundsMin = getBoundsMin();
2412            Parameter<Length> range = boundsMax.get(dim).minus(boundsMin.get(dim));
2413
2414            if (thisDim > 2) {
2415                //  We have at least 3 dimensions.
2416
2417                //  Create a vector pointing in the indicated direction.
2418                MutablePoint p = MutablePoint.newInstance(thisDim);
2419                p.set(dim, Parameter.valueOf(1, SI.METER));
2420                Vector<Dimensionless> uv = Vector.valueOf(p).toUnitVector();
2421
2422                //  Bias the bounds off 1% in the specified dimension
2423                if (max) {
2424                    p.set(boundsMax);
2425                    p.set(dim, p.get(dim).plus(range.times(0.01)));
2426                } else {
2427                    p.set(boundsMin);
2428                    p.set(dim, p.get(dim).minus(range.times(0.01)));
2429                }
2430
2431                //  Create a plane 1% of the diagonal off the bounds.
2432                Plane plane = Plane.valueOf(uv, p.immutable());
2433
2434                //  Return the closest point to the specified plane.
2435                SubrangePoint pnt = getClosest(plane, tol);
2436                sout = pnt.getParPosition().getValue(0);
2437
2438            } else if (thisDim == 2) {
2439                //  We have a 2D curve.
2440
2441                //  Create two points biased off the bounds in the specified direction.
2442                MutablePoint p0, p1;
2443                if (max) {
2444                    p0 = MutablePoint.valueOf(boundsMax);
2445                    p0.set(dim, p0.get(dim).plus(range.times(0.01)));
2446                    p1 = MutablePoint.valueOf(boundsMin);
2447                    p1.set(dim, p0.get(dim));
2448                } else {
2449                    p0 = MutablePoint.valueOf(boundsMin);
2450                    p0.set(dim, p0.get(dim).minus(range.times(0.01)));
2451                    p1 = MutablePoint.valueOf(boundsMax);
2452                    p1.set(dim, p0.get(dim));
2453                }
2454
2455                //  Create a line on the desired side of the curve.
2456                LineSeg seg = LineSeg.valueOf(p0.immutable(), p1.immutable());
2457
2458                //  Find the closest point between this curve the line segment.
2459                PointString<SubrangePoint> str = getClosest(seg, tol);
2460                SubrangePoint p = str.getFirst();
2461                sout = p.getParPosition().getValue(0);
2462
2463            } else {
2464                //  We have a 1 dimensional curve.
2465
2466                //  Bias the bounds off 1% in the specified dimension
2467                Point pp;
2468                if (max) {
2469                    pp = boundsMax.plus(range.times(0.01));
2470                } else {
2471                    pp = boundsMin.minus(range.times(0.01));
2472                }
2473
2474                //  Return the closest point to the specified point.
2475                sout = getClosest(pp, tol).getParPosition().getValue(0);
2476            }
2477        } finally {
2478            StackContext.exit();
2479        }
2480
2481        return getPoint(sout);
2482    }
2483
2484    /**
2485     * Return the signed area of the region enclosed or subtended by a planar curve
2486     * relative to the specified reference point. This method assumes that the curve is
2487     * planar and will give incorrect results if the curve is not planar. Also, the input
2488     * reference point is assumed to lie in the plane of the curve, it if is not, this
2489     * will give incorrect results. If the curve is closed, then the reference point is
2490     * arbitrary (as long as it is in the plane of the curve). If the curve is open, then
2491     * the area returned is that subtended by the curve with straight lines from each end
2492     * of the curve to the input point.
2493     *
2494     * @param refPnt The reference point used for computing the enclosed or subtended area
2495     *               of the curve.
2496     * @param eps    The desired fractional accuracy of the area calculated.
2497     * @return The area enclosed or subtended by this curve (assuming that this curve is
2498     *         planar and that the reference point lies in the plane of this curve).
2499     */
2500    @Override
2501    public Parameter<Area> getEnclosedArea(GeomPoint refPnt, double eps) {
2502        Unit<Area> areaUnit = getUnit().pow(2).asType(Area.class);
2503        
2504        //  If we only have 1 dimension, then there can be no area.
2505        int numDims = getPhyDimension();
2506        if (numDims == 1)
2507            return Parameter.ZERO_AREA.to(areaUnit);
2508
2509        double area = 0;
2510        StackContext.enter();
2511        try {
2512            //  Make sure the point and curve have the same dimensions.
2513            if (refPnt.getPhyDimension() > numDims) {
2514                throw new IllegalArgumentException(MessageFormat.format(
2515                        RESOURCES.getString("sameOrFewerDimErr"), "refPnt", "curve"));
2516            } else if (refPnt.getPhyDimension() < numDims) {
2517                refPnt = refPnt.toDimension(numDims);
2518            }
2519
2520            //  Convert the reference point to the same units as this curve.
2521            refPnt = refPnt.to(getUnit());
2522
2523            EnclAreaEvaluatable areaEvaluator = EnclAreaEvaluatable.newInstance(this, refPnt);
2524            try {
2525                //  Integrate to find the subtended area:   A = int_0^1{ p0p(t) x p'(t) ds}
2526                area = Quadrature.adaptLobatto(areaEvaluator, 0, 1, eps) * 0.5;
2527
2528            } catch (RootException | IntegratorException e) {
2529                //  Fall back on Gaussian Quadrature.
2530                Logger.getLogger(AbstractCurve.class.getName()).log(Level.WARNING,
2531                        "Trying Gaussian Quadrature in getEnclosedArea()...");
2532                try {
2533                    area = Quadrature.gaussLegendre_Wx1N40(areaEvaluator, 0, 1) * 0.5;
2534
2535                } catch (RootException err) {
2536                    Logger.getLogger(AbstractCurve.class.getName()).log(Level.WARNING,
2537                            "Could not calculate enclosed area; using 0.", err);
2538                    area = 0;
2539                }
2540            }
2541        } finally {
2542            StackContext.exit();
2543        }
2544
2545        return Parameter.valueOf(area, areaUnit);
2546    }
2547
2548    /**
2549     * Return <code>true</code> if this curve is degenerate (i.e.: has length less than
2550     * the specified tolerance).
2551     *
2552     * @param tol The tolerance for determining if this curve is degenerate. May not be
2553     *            null.
2554     */
2555    @Override
2556    public boolean isDegenerate(Parameter<Length> tol) {
2557        requireNonNull(tol);
2558        Parameter<Length> cLength = getArcLength(GTOL * 100);
2559        return cLength.compareTo(tol) <= 0;
2560    }
2561
2562    /**
2563     * Return <code>true</code> if this curve is planar or <code>false</code> if it is
2564     * not.
2565     *
2566     * @param tol The geometric tolerance to use in determining if the curve is planar.
2567     */
2568    @Override
2569    public boolean isPlanar(Parameter<Length> tol) {
2570        //  If the curve is less than 3D, then it must be planar.
2571        int numDims = getPhyDimension();
2572        if (numDims <= 2)
2573            return true;
2574
2575        //  Is this curve degenerate?
2576        if (isDegenerate(tol))
2577            return true;
2578
2579        //  Is the curve linear?
2580        if (isLine(tol))
2581            return true;
2582
2583        StackContext.enter();
2584        try {
2585            //  Create a list of segment starting par. positions (and the last end position, 1).
2586            FastTable<Double> segParams = getSegmentParameters();
2587
2588            //  Form a plane from three points on the curve.
2589            Point p0 = getRealPoint(0);
2590            GeomVector<Dimensionless> nhat = getPlaneNormal(p0);
2591
2592            // Number of points extracted on each segment to determine type.
2593            int nPnts = getNumPointsPerSegment();
2594            int nPntsm1 = nPnts - 1;
2595
2596            //  Loop over all the segments.
2597            int numSegs = segParams.size() - 1;
2598            for (int seg = 0; seg < numSegs; ++seg) {
2599
2600                //  Sample at 2*(p+1) points on each segment.
2601                double sSeg = segParams.get(seg);
2602                double sNext = segParams.get(seg + 1);
2603                double sRange = sNext - sSeg;
2604                for (int i = 0; i < nPntsm1; ++i) {
2605                    double s = sSeg + (i * sRange) / nPntsm1;
2606                    Point pi = getRealPoint(s);
2607
2608                    //  Find deviation from a plane for each point by projecting a vector from p0 to pi
2609                    //  onto the normal for the arbitrarily chosen points on the curve.
2610                    //  If the projection is anything other than zero, the points are not
2611                    //  in the same plane.
2612                    Vector<Length> p1pi = pi.minus(p0).toGeomVector();
2613                    if (!p1pi.mag().isApproxZero()) {
2614                        Parameter proj = p1pi.dot(nhat);
2615                        if (proj.isLargerThan(tol))
2616                            return false;
2617                    }
2618                }
2619            }
2620
2621            //  Must be planar.
2622            return true;
2623
2624        } finally {
2625            StackContext.exit();
2626        }
2627    }
2628
2629    /**
2630     * Return a normal vector for the curve from an input point and 2 arbitrarily chosen
2631     * points on a > 2D curve. If the input point is in the plane of the curve, then this
2632     * will return the normal vector for the plane of the curve.
2633     */
2634    private Vector<Dimensionless> getPlaneNormal(GeomPoint p0) {
2635        StackContext.enter();
2636        try {
2637            Vector p0v = p0.toGeomVector();
2638
2639            Point p1 = getRealPoint(0.38196601125);
2640            Point p2 = getRealPoint(0.61803398875);
2641            Vector<Length> p1v = p1.toGeomVector();
2642            Vector<Length> p2v = p2.toGeomVector();
2643            Vector<Length> v1 = p1v.minus(p0v);
2644            Vector<Length> v2 = p2v.minus(p0v);
2645            Vector n = v1.cross(v2);
2646            Vector<Dimensionless> nhat;
2647            if (n.magValue() > GTOL) {
2648                nhat = n.toUnitVector();
2649                nhat.setOrigin(Point.newInstance(getPhyDimension()));
2650            } else {
2651                //  Try a different set of points on the curve.
2652                p1 = getRealPoint(0.25);
2653                p2 = getRealPoint(0.5);
2654                p1v = p1.toGeomVector();
2655                p2v = p2.toGeomVector();
2656                v1 = p1v.minus(p0v);
2657                v2 = p2v.minus(p0v);
2658                n = v1.cross(v2);
2659                if (n.magValue() > GTOL) {
2660                    nhat = n.toUnitVector();
2661                    nhat.setOrigin(Point.newInstance(getPhyDimension()));
2662                } else {
2663                    //  Points from input curve form a straight line.
2664                    nhat = GeomUtil.calcYHat(p1v.minus(p2v).toUnitVector());
2665                    nhat = v2.cross(nhat).toUnitVector();
2666                }
2667            }
2668
2669            return StackContext.outerCopy(nhat);
2670
2671        } finally {
2672            StackContext.exit();
2673        }
2674    }
2675
2676    /**
2677     * Returns <code>true</code> if this curve is a line to within the specified
2678     * tolerance.
2679     *
2680     * @param tol The tolerance for determining if this curve is a line.
2681     */
2682    @Override
2683    public boolean isLine(Parameter<Length> tol) {
2684        //  Reference:  Piegl, L.A., Tiller, W., "Computing Offsets of NURBS Curves and Surfaces",
2685        //              Computer Aided Design 31, 1999, pgs. 147-156.
2686        requireNonNull(tol);
2687        
2688        StackContext.enter();
2689        try {
2690            //  Extract the curve end points.
2691            Point p0 = getRealPoint(0), p1 = getRealPoint(1);
2692
2693            //  If the end points are coincident, then it can't be a line.
2694            Vector<Length> tv = p1.minus(p0).toGeomVector();
2695            if (!tv.mag().isGreaterThan(tol))
2696                return false;
2697
2698            //  Get a unit vector for the potential line.
2699            Vector<Dimensionless> uv = tv.toUnitVector();
2700
2701            //  Create a list of segment starting par. positions (and the last end position, 1).
2702            FastTable<Double> bParams = getSegmentParameters();
2703
2704            // Number of points extracted on each segment to determine type.
2705            int nPnts = getNumPointsPerSegment();
2706            int nPntsm1 = nPnts - 1;
2707
2708            //  Loop over all the segments.
2709            Parameter c2 = tv.dot(tv);
2710            int numSegs = bParams.size() - 1;
2711            for (int seg = 0; seg < numSegs; ++seg) {
2712
2713                //  Sample at nPnts points on each segment.
2714                double sSeg = bParams.get(seg);
2715                double sNext = bParams.get(seg + 1);
2716                double sRange = sNext - sSeg;
2717                for (int i = 0; i < nPntsm1; ++i) {
2718                    double s = sSeg + (i * sRange) / nPntsm1;
2719                    Point pi = getRealPoint(s);
2720
2721                    //  Check distance of this point from infinite line p0-uv.
2722                    if (GeomUtil.pointLineDistance(pi, p0, uv).isGreaterThan(tol))
2723                        return false;
2724
2725                    //  See if the point falls in the segment p0 to p1.
2726                    Vector<Length> w = pi.minus(p0).toGeomVector();
2727                    Parameter c1 = w.dot(tv);
2728                    if (c1.getValue() < 0)
2729                        return false;
2730                    if (c2.isLessThan(c1))
2731                        return false;
2732                }
2733
2734            }
2735
2736            //  Must be a straight line.
2737            return true;
2738
2739        } finally {
2740            StackContext.exit();
2741        }
2742    }
2743
2744    /**
2745     * Returns <code>true</code> if this curve is a circular arc to within the specified
2746     * tolerance.
2747     *
2748     * @param tol The tolerance for determining if this curve is a circular arc.
2749     */
2750    @Override
2751    public boolean isCircular(Parameter<Length> tol) {
2752        //  Reference:  Piegl, L.A., Tiller, W., "Computing Offsets of NURBS Curves and Surfaces",
2753        //              Computer Aided Design 31, 1999, pgs. 147-156.
2754
2755        //  Can't be degenerate.
2756        if (isDegenerate(tol))
2757            return false;
2758
2759        //  Must be planar.
2760        if (!isPlanar(tol))
2761            return false;
2762
2763        StackContext.enter();
2764        try {
2765            //  Create a list of segment starting par. positions (and the last end position, 1).
2766            FastTable<Double> bParams = getSegmentParameters();
2767
2768            //  Get the curvature and principal normal at the start of the curve.
2769            Parameter k0 = getCurvature(0);
2770            if (k0.isApproxZero())
2771                return false;
2772            Vector n0 = getPrincipalNormal(0);
2773
2774            //  Compute the center of the circle.
2775            Point pc0 = Point.valueOf(n0.divide(k0)).plus(n0.getOrigin());
2776
2777            // Number of points extracted on each segment to determine type.
2778            int nPnts = getNumPointsPerSegment();
2779            int nPntsm1 = nPnts - 1;
2780
2781            //  Loop over all the segments.
2782            int numSegs = bParams.size() - 1;
2783            for (int seg = 0; seg < numSegs; ++seg) {
2784
2785                //  Sample at nPnts points on each segment.
2786                double sSeg = bParams.get(seg);
2787                double sNext = bParams.get(seg + 1);
2788                double sRange = sNext - sSeg;
2789                for (int i = 0; i < nPntsm1; ++i) {
2790                    double s = sSeg + (i * sRange) / nPntsm1;
2791
2792                    //  Get the local curvature and principal normal vector.
2793                    Parameter k = getCurvature(s);
2794                    if (k.isApproxZero())
2795                        return false;
2796                    Vector n = getPrincipalNormal(s);
2797
2798                    //  Compute the local center of curvature.
2799                    Point pc = Point.valueOf(n.divide(k)).plus(n.getOrigin());
2800                    Parameter<Length> dist = pc0.distance(pc);
2801                    if (dist.isGreaterThan(tol))
2802                        return false;
2803                }
2804            }
2805
2806            //  The curve must be circular.
2807            return true;
2808
2809        } finally {
2810            StackContext.exit();
2811        }
2812    }
2813
2814    /**
2815     * Return a list containing parameters at the start of logical segments of this curve.
2816     * The first element in this list must always be 0.0 and the last element 1.0. These
2817     * should be good places to break the curve up into pieces for analysis, but otherwise
2818     * are arbitrary. Subclasses should override this method to provide better segment
2819     * starting parameters.
2820     *
2821     * The default implementation always splits the curve into two pieces and returns the
2822     * following parameters [0, 0.50, 1].
2823     *
2824     * @return A list containing parameters at the start of logical segments of this
2825     *         curve.
2826     */
2827    protected FastTable<Double> getSegmentParameters() {
2828        FastTable<Double> knots = FastTable.newInstance();
2829        knots.add(0.0);
2830        knots.add(0.5);
2831        knots.add(1.0);
2832        return knots;
2833    }
2834
2835    /**
2836     * Return the number of points per segment that should be used when analyzing curve
2837     * segments by certain methods. Subclasses should override this method to provide a
2838     * more appropriate number of points per segment.
2839     *
2840     * The default implementation always returns 8.
2841     *
2842     * @return The number of points per segment that should be used when analyzing curve
2843     *         segments by certain methods.
2844     */
2845    protected int getNumPointsPerSegment() {
2846        return 8;   //  2*(3 + 1)
2847    }
2848
2849    /**
2850     * Converts the input parametric position on a curve segment with the specified range
2851     * of parametric positions into a parametric position on the parent curve.
2852     *
2853     * @param s  The parametric position on the curve segment to be converted.
2854     * @param s0 The starting parametric position of the curve segment on the parent
2855     *           curve.
2856     * @param s1 The ending parametric position of the curve segment on the parent curve.
2857     * @return The input segment parametric position converted to the parent curve.
2858     */
2859    protected static double segmentPos2Parent(double s, double s0, double s1) {
2860        double s_out = s0 + (s1 - s0) * s;
2861        if (s_out > 1) {            //  Deal with roundoff issues.
2862            s_out = 1;
2863        } else if (s_out < 0) {
2864            s_out = 0;
2865        }
2866        return s_out;
2867    }
2868
2869    /**
2870     * An Evaulatable1D that returns the magnitude of the tangent vector at the specified
2871     * parametric location. This is used to find the arc length of a curve. This is
2872     * <code>f(x)</code> in:
2873     * <code>L = int_0^1 {f(x)dx} = int_0^1 { sqrt(ps(s) DOT ps(s)) ds }</code>.
2874     */
2875    private static class ArcLengthEvaluatable extends AbstractEvaluatable1D {
2876
2877        private Curve _curve;
2878
2879        public static ArcLengthEvaluatable newInstance(Curve curve) {
2880            if (curve instanceof GeomTransform)
2881                curve = curve.copyToReal();
2882            ArcLengthEvaluatable o = FACTORY.object();
2883            o._curve = curve;
2884            return o;
2885        }
2886
2887        @Override
2888        public double function(double s) throws RootException {
2889            StackContext.enter();
2890            try {
2891                Vector<Length> ps = _curve.getSDerivative(s, 1);        //  The dimensional tangent vector.
2892                double psv = ps.normValue();                            //  sqrt(ps DOT ps) = |ps|
2893                return psv;
2894            } finally {
2895                StackContext.exit();
2896            }
2897        }
2898
2899        public static void recycle(ArcLengthEvaluatable instance) {
2900            FACTORY.recycle(instance);
2901        }
2902
2903        private ArcLengthEvaluatable() { }
2904
2905        private static final ObjectFactory<ArcLengthEvaluatable> FACTORY = new ObjectFactory<ArcLengthEvaluatable>() {
2906            @Override
2907            protected ArcLengthEvaluatable create() {
2908                return new ArcLengthEvaluatable();
2909            }
2910        };
2911
2912    }
2913
2914    /**
2915     * An Evaulatable1D that is used to find a point on a curve that has a tangent vector
2916     * that points to the input point.
2917     */
2918    private static class TangentPointEvaluatable extends AbstractEvaluatable1D {
2919
2920        private Curve _curve;
2921        private Point _q;
2922
2923        public static TangentPointEvaluatable newInstance(Curve curve, GeomPoint q) {
2924            if (curve instanceof GeomTransform)
2925                curve = curve.copyToReal();
2926            TangentPointEvaluatable o = FACTORY.object();
2927            o._curve = curve;
2928            o._q = q.immutable();
2929            return o;
2930        }
2931
2932        /**
2933         * Calculates: 1 - ((q - p(s))/|q - p(s)|) DOT t(s) where t(s) is the tangent
2934         * vector on the curve at s, p(s) is the point on the curve at s and q is the
2935         * point we are finding a tangent on the curve to.
2936         */
2937        @Override
2938        public double function(double s) throws RootException {
2939            StackContext.enter();
2940            try {
2941
2942                Vector<Dimensionless> t = _curve.getTangent(s);
2943                Point p = _curve.getRealPoint(s);
2944                Vector pq = _q.minus(p).toGeomVector();
2945                double tDOTpq = 0;
2946                if (!pq.norm().isApproxZero()) {
2947                    pq = pq.toUnitVector();
2948                    tDOTpq = t.dot(pq).getValue();
2949                }
2950                return 1.0 - MathLib.abs(tDOTpq);
2951
2952            } finally {
2953                StackContext.exit();
2954            }
2955        }
2956
2957        public static void recycle(TangentPointEvaluatable instance) {
2958            FACTORY.recycle(instance);
2959        }
2960
2961        private TangentPointEvaluatable() { }
2962
2963        private static final ObjectFactory<TangentPointEvaluatable> FACTORY = new ObjectFactory<TangentPointEvaluatable>() {
2964            @Override
2965            protected TangentPointEvaluatable create() {
2966                return new TangentPointEvaluatable();
2967            }
2968        };
2969    }
2970
2971    /**
2972     * An Evaulatable1D that is used to find a point on a curve that has the specified arc
2973     * length.
2974     */
2975    private static class ArcPosEvaluatable extends AbstractEvaluatable1D {
2976
2977        private Curve _curve;
2978        private double _eps;
2979        public boolean firstPass;
2980
2981        /**
2982         * Return a new ArcPosEvaluatable function.
2983         *
2984         * @param curve        The curve to find the arc-length along.
2985         * @param targetLength The target length, in meters, to search for.
2986         * @param eps          The desired fractional accuracy on the arc-length.
2987         */
2988        @SuppressWarnings("null")
2989        public static ArcPosEvaluatable newInstance(Curve curve, double targetLength, double eps) {
2990            if (curve instanceof GeomTransform)
2991                curve = curve.copyToReal();
2992
2993            ArcPosEvaluatable o = FACTORY.object();
2994            o._curve = curve.to(SI.METER);  //  Convert everything to meters.
2995            o._eps = eps;
2996            o.setInitialValue(targetLength);
2997            o.firstPass = false;
2998
2999            return o;
3000        }
3001
3002        /**
3003         * Calculates: arcLength(s) - targetArcLength. When this is zero, then the arc
3004         * length position has been found.
3005         */
3006        @Override
3007        public double function(double s) throws RootException {
3008            //  Don't compute anything if the 1st pass flag is set.
3009            if (firstPass) {
3010                firstPass = false;
3011                return 0;
3012            }
3013            StackContext.enter();
3014            try {
3015
3016                double arcLength = _curve.getArcLength(0, s, _eps).getValue();
3017                return arcLength - getInitialValue();
3018
3019            } finally {
3020                StackContext.exit();
3021            }
3022        }
3023
3024        public static void recycle(ArcPosEvaluatable instance) {
3025            FACTORY.recycle(instance);
3026        }
3027
3028        private ArcPosEvaluatable() { }
3029
3030        private static final ObjectFactory<ArcPosEvaluatable> FACTORY = new ObjectFactory<ArcPosEvaluatable>() {
3031            @Override
3032            protected ArcPosEvaluatable create() {
3033                return new ArcPosEvaluatable();
3034            }
3035        };
3036    }
3037
3038    /**
3039     * An Evaulatable1D that is used to find the intersection of a curve and a plane.
3040     */
3041    private static class PlaneIntersectEvaluatable implements Evaluatable1D {
3042
3043        private Curve _curve;
3044        private GeomPlane _plane;
3045        public boolean firstPass;
3046
3047        public static PlaneIntersectEvaluatable newInstance(Curve curve, GeomPlane plane) {
3048            if (curve instanceof GeomTransform)
3049                curve = curve.copyToReal();
3050
3051            //  Create the object.
3052            PlaneIntersectEvaluatable o = FACTORY.object();
3053            o._curve = curve;
3054            o._plane = plane;
3055            o.firstPass = false;
3056
3057            return o;
3058        }
3059
3060        /**
3061         * Calculates: <code>n dot (p(s) - pp)</code> which is the signed distance from
3062         * the point to the plane. When this is zero, then an intersection has been found.
3063         */
3064        @Override
3065        public double function(double s) throws RootException {
3066            //  Don't compute anything if the 1st pass flag is set.
3067            if (firstPass) {
3068                return 0;
3069            }
3070
3071            StackContext.enter();
3072            try {
3073
3074                Point ps = _curve.getRealPoint(s);
3075                GeomPoint pp = _plane.getRefPoint();
3076                Vector<Length> psmpp = ps.minus(pp).toGeomVector();
3077                GeomVector n = _plane.getNormal();
3078
3079                Parameter d = n.dot(psmpp);
3080
3081                return d.getValue();
3082
3083            } finally {
3084                StackContext.exit();
3085            }
3086        }
3087
3088        /**
3089         * Calculates: <code>d(n dot (p(s) - pp))/ds</code> which is the slope of the
3090         * signed distance from the point to the plane wrt to s.
3091         */
3092        @Override
3093        public double derivative(double s) throws RootException {
3094            //  Don't compute anything if the 1st pass flag is set.
3095            if (firstPass) {
3096                firstPass = false;
3097                return 0;
3098            }
3099
3100            StackContext.enter();
3101            try {
3102
3103                Vector dpds = _curve.getSDerivative(s, 1);
3104                GeomVector n = _plane.getNormal();
3105
3106                Parameter d = n.dot(dpds);
3107
3108                return d.getValue();
3109
3110            } finally {
3111                StackContext.exit();
3112            }
3113        }
3114
3115        public static void recycle(PlaneIntersectEvaluatable instance) {
3116            FACTORY.recycle(instance);
3117        }
3118
3119        private PlaneIntersectEvaluatable() { }
3120
3121        private static final ObjectFactory<PlaneIntersectEvaluatable> FACTORY = new ObjectFactory<PlaneIntersectEvaluatable>() {
3122            @Override
3123            protected PlaneIntersectEvaluatable create() {
3124                return new PlaneIntersectEvaluatable();
3125            }
3126        };
3127    }
3128
3129    /**
3130     * An Evaulatable1D that is used to find a parametric position on a curve that forms a
3131     * vector with the target point that is perpendicular to the curves tangent vector at
3132     * that position. Such a perpendicular vector indicates a point that is locally either
3133     * the closest or furthest away from the target point. This is then used with a 1D
3134     * root finder to isolate all the perpendicular points on the curve.
3135     */
3136    protected static class PerpPointEvaluatable implements Evaluatable1D {
3137
3138        private Curve _curve;
3139        private Vector<Length> _point;
3140        public boolean firstPass;
3141        private Vector<Length> _pmq;
3142        private List<Vector<Length>> _ders;
3143
3144        public static PerpPointEvaluatable newInstance(Curve curve, GeomPoint point) {
3145            if (curve instanceof GeomTransform)
3146                curve = curve.copyToReal();
3147
3148            PerpPointEvaluatable o = FACTORY.object();
3149            o._curve = curve;
3150            o._point = point.toGeomVector();
3151            o.firstPass = false;
3152
3153            return o;
3154        }
3155
3156        /**
3157         * Calculates (p(s) - q) DOT ps(s) (dot product between the vector formed by the
3158         * target point and a point on the curve at "s" and the tangent vector at "s".
3159         * Where this is zero, a perpendicular point has been found. A perpendicular point
3160         * is either the closest or furthest away from the target point in a local region.
3161         *
3162         * @param s The parametric position on the curve being evaluated.
3163         */
3164        @Override
3165        public double function(double s) {
3166            //  Don't compute anything if the 1st pass flag is set.
3167            if (firstPass) {
3168                return 0;
3169            }
3170
3171            //  Get derivatives up to 2nd order.  These are also used in derivative().
3172            _ders = _curve.getSDerivatives(s, 2);
3173
3174            //  Calculate: p(s) - q.  This is also used in derivative().
3175            Vector<Length> p = _ders.get(0);
3176            _pmq = p.minus(_point);
3177            if (_pmq.mag().isApproxZero())
3178                return 0;
3179
3180            Vector ps = _ders.get(1);
3181            Parameter cosa = _pmq.dot(ps);  //  Projections of (p-q) onto tangent vector.
3182
3183            return cosa.getValue();
3184        }
3185
3186        /**
3187         * Calculates: <code>d((p(s) - q) DOT ps(s))/ds</code> which is the slope of the
3188         * projection of the target point to curve point vector onto the tangent vector.
3189         *
3190         * @param s The parametric position on the curve being evaluated.
3191         */
3192        @Override
3193        public double derivative(double s) {
3194            //  Don't compute anything if the 1st pass flag is set.
3195            if (firstPass) {
3196                firstPass = false;
3197                return 0;
3198            }
3199
3200            /*  d(fg) = df*g + f*dg
3201             d(p(s) - q)/ds = ps(s);    d(ps(s))/ds = pss(s)
3202             d((p(s) - q) DOT ps(s))/ds = ps(s) DOT ps(s) + (p(s) - q) DOT pss(s)
3203             ((p(s) - q) DOT ps(s))/ds = |ps(s)|^2 + (p(s) - q) DOT pss(s)
3204             Derivatives and p(s) - q calculated in "function()".
3205             */
3206            Vector ps = _ders.get(1);
3207            Vector pss = _ders.get(2);
3208            FastTable.recycle((FastTable)_ders);
3209
3210            Parameter value = ps.dot(ps).plus(_pmq.dot(pss));
3211
3212            Vector.recycle(_pmq);
3213            Vector.recycle(ps);
3214            Vector.recycle(pss);
3215
3216            return value.getValue();
3217        }
3218
3219        public static void recycle(PerpPointEvaluatable instance) {
3220            FACTORY.recycle(instance);
3221        }
3222
3223        private PerpPointEvaluatable() { }
3224
3225        private static final ObjectFactory<PerpPointEvaluatable> FACTORY = new ObjectFactory<PerpPointEvaluatable>() {
3226            @Override
3227            protected PerpPointEvaluatable create() {
3228                return new PerpPointEvaluatable();
3229            }
3230
3231            @Override
3232            protected void cleanup(PerpPointEvaluatable obj) {
3233                obj._curve = null;
3234                obj._point = null;
3235                obj._pmq = null;
3236                obj._ders = null;
3237            }
3238        };
3239
3240    }
3241
3242    /**
3243     * An Evaulatable1D that is used to find a parametric position on a curve that forms a
3244     * vector with the target point (which is the closest point on the target curve or
3245     * surface) that is perpendicular to the curve's tangent vector at that position. Such
3246     * a perpendicular vector indicates a point that is locally either the closest or
3247     * furthest away from the target point. This is then used with a 1D root finder to
3248     * isolate all the potential nearest/farthest points on the two curves.
3249     */
3250    protected static class PerpPointParGeomEvaluatable extends AbstractEvaluatable1D {
3251
3252        private Curve _thisCurve;
3253        private ParametricGeometry _thatParGeom = null;
3254        private double _tol;
3255        public boolean firstPass;
3256
3257        public static PerpPointParGeomEvaluatable newInstance(Curve thisCurve, ParametricGeometry thatParGeom, double tol) {
3258            if (thisCurve instanceof GeomTransform)
3259                thisCurve = thisCurve.copyToReal();
3260
3261            PerpPointParGeomEvaluatable o = FACTORY.object();
3262            o._thisCurve = thisCurve;
3263            o._thatParGeom = thatParGeom;
3264            o._tol = tol;
3265            o.firstPass = false;
3266
3267            return o;
3268        }
3269
3270        /**
3271         * Calculates (p(s) - q) DOT ps(s) (dot product between the vector formed by the
3272         * target point (closest point on target parametric geometry) and a point on this
3273         * curve at "s" and the tangent vector at "s". Where this is zero, a perpendicular
3274         * point has been found. A perpendicular point is either the closest or furthest
3275         * away from the target point in a local region.
3276         *
3277         * @param s The parametric position on the curve being evaluated.
3278         * @return The dot product between the vector formed by the target point and a
3279         *         point on this curve.
3280         */
3281        @Override
3282        public double function(double s) {
3283            //  Don't compute anything if the 1st pass flag is set.
3284            if (firstPass) {
3285                firstPass = false;
3286                return 0;
3287            }
3288
3289            StackContext.enter();
3290            try {
3291                //  Get derivatives up to 2nd order.
3292                FastTable<Vector> ders = (FastTable)_thisCurve.getSDerivatives(s, 2);
3293
3294                //  Extract the point on this curve at "s".
3295                Vector<Length> p = ders.get(0);
3296                Point pnt = Point.valueOf(p);
3297
3298                //  Get the closest point on the target parametric geometry.
3299                Vector<Length> q = _thatParGeom.getClosest(pnt, _tol).toGeomVector();
3300
3301                //  Calculate: p(s) - q.
3302                Vector<Length> pmq = p.minus(q);
3303                if (pmq.mag().isApproxZero()) {
3304                    return 0;
3305                }
3306
3307                Vector ps = ders.get(1);
3308                Parameter cosa = pmq.dot(ps);   //  Projection of (p-q) onto tangent vector.
3309
3310                return cosa.getValue();
3311
3312            } finally {
3313                StackContext.exit();
3314            }
3315        }
3316
3317        public static void recycle(PerpPointParGeomEvaluatable instance) {
3318            FACTORY.recycle(instance);
3319        }
3320
3321        private PerpPointParGeomEvaluatable() { }
3322
3323        private static final ObjectFactory<PerpPointParGeomEvaluatable> FACTORY = new ObjectFactory<PerpPointParGeomEvaluatable>() {
3324            @Override
3325            protected PerpPointParGeomEvaluatable create() {
3326                return new PerpPointParGeomEvaluatable();
3327            }
3328        };
3329
3330    }
3331
3332    /**
3333     * An Evaulatable1D that is used to find a parametric position on a curve that forms a
3334     * vector with the target point (which is the closest point on the target plane) that
3335     * is perpendicular to the curve's tangent vector at that position. Such a
3336     * perpendicular vector indicates a point that is locally either the closest or
3337     * furthest away from the target point. This is then used with a 1D root finder to
3338     * isolate all the potential nearest/farthest points on the curve to the plane.
3339     */
3340    protected static class PerpPointPlaneEvaluatable extends AbstractEvaluatable1D {
3341
3342        private Curve _thisCurve;
3343        private GeomPlane _thatPlane = null;
3344        public boolean firstPass;
3345
3346        public static PerpPointPlaneEvaluatable newInstance(Curve thisCurve, GeomPlane thatPlane) {
3347            if (thisCurve instanceof GeomTransform)
3348                thisCurve = thisCurve.copyToReal();
3349
3350            PerpPointPlaneEvaluatable o = FACTORY.object();
3351            o._thisCurve = thisCurve;
3352            o._thatPlane = thatPlane;
3353            o.firstPass = false;
3354
3355            return o;
3356        }
3357
3358        /**
3359         * Calculates (p(s) - q) DOT ps(s) (dot product between the vector formed by the
3360         * target point (closest point on target plane) and a point on this curve at "s"
3361         * and the tangent vector at "s". Where this is zero, a perpendicular point has
3362         * been found. A perpendicular point is either the closest or furthest away from
3363         * the target point in a local region.
3364         *
3365         * @param s The parametric position on this curve being evaluated.
3366         * @return The dot product between the vector formed by the target point and a
3367         *         point on this curve
3368         */
3369        @Override
3370        public double function(double s) {
3371            //  Don't compute anything if the 1st pass flag is set.
3372            if (firstPass) {
3373                firstPass = false;
3374                return 0;
3375            }
3376
3377            StackContext.enter();
3378            try {
3379                //  Get derivatives up to 2nd order.
3380                FastTable<Vector> ders = (FastTable)_thisCurve.getSDerivatives(s, 2);
3381
3382                //  Extract the point on this curve at "s".
3383                Vector<Length> p = ders.get(0);
3384                Point pnt = Point.valueOf(p);
3385
3386                //  Get the closest point on the target plane geometry.
3387                Vector<Length> q = _thatPlane.getClosest(pnt).toGeomVector();
3388
3389                //  Calculate: p(s) - q.
3390                Vector<Length> pmq = p.minus(q);
3391                if (pmq.mag().isApproxZero()) {
3392                    return 0;
3393                }
3394
3395                Vector ps = ders.get(1);
3396                Parameter cosa = pmq.dot(ps);   //  Projection of (p-q) onto tangent vector.
3397
3398                return cosa.getValue();
3399
3400            } finally {
3401                StackContext.exit();
3402            }
3403        }
3404
3405        public static void recycle(PerpPointPlaneEvaluatable instance) {
3406            FACTORY.recycle(instance);
3407        }
3408
3409        private PerpPointPlaneEvaluatable() { }
3410
3411        private static final ObjectFactory<PerpPointPlaneEvaluatable> FACTORY = new ObjectFactory<PerpPointPlaneEvaluatable>() {
3412            @Override
3413            protected PerpPointPlaneEvaluatable create() {
3414                return new PerpPointPlaneEvaluatable();
3415            }
3416        };
3417
3418    }
3419
3420    /**
3421     * A class used to intersect an infinite line and a curve.
3422     */
3423    private static class LineCurveIntersect {
3424
3425        private static final int MAXDEPTH = 20;
3426        private Curve _thisCurve;
3427        private Parameter<Length> _tol;
3428        private Point _L0;
3429        private Vector<Length> _Ldir;
3430        private FastTable<Double> _sParams;
3431
3432        @SuppressWarnings("null")
3433        public static LineCurveIntersect newInstance(Curve thisCurve,
3434                Parameter<Length> tol, GeomPoint L0, GeomVector<Length> Ldir) {
3435
3436            //  The magnitude of Ldir doesn't matter, but it can not be "Dimensionless".
3437            if (Ldir.getUnit().equals(Dimensionless.UNIT)) {
3438                Ldir = Vector.valueOf(Ldir.toFloat64Vector(), SI.METER);
3439            }
3440
3441            if (thisCurve instanceof GeomTransform)
3442                thisCurve = thisCurve.copyToReal();
3443
3444            //  Create an instance of this class and fill in the class data.
3445            LineCurveIntersect o = FACTORY.object();
3446            o._thisCurve = thisCurve;
3447            o._tol = tol.to(thisCurve.getUnit());
3448            o._L0 = L0.immutable().to(thisCurve.getUnit());
3449            o._Ldir = Ldir.immutable();
3450            o._sParams = FastTable.newInstance();
3451
3452            return o;
3453        }
3454
3455        /**
3456         * Method that actually carries out the intersection adding the intersection
3457         * points to the list provided in the constructor.
3458         *
3459         * @return A PointString containing zero or more subrange points made by the
3460         *         intersection of the curve with the specified line. If no intersection
3461         *         is found, an empty PointString is returned.
3462         */
3463        public PointString<SubrangePoint> intersect() {
3464            divideAndConquerLine(1, _thisCurve, 0, 1);
3465
3466            PointString<SubrangePoint> output = PointString.newInstance();
3467            int size = _sParams.size();
3468            for (int i = 0; i < size; ++i) {
3469                double s = _sParams.get(i);
3470                output.add(SubrangePoint.newInstance(_thisCurve, s));
3471            }
3472            FastTable.recycle(_sParams);
3473            _sParams = null;
3474
3475            return output;
3476        }
3477
3478        /**
3479         * Uses a recursive "Divide and Conquer" approach to intersecting an infinite line
3480         * with a curve. On each call, the following happens:
3481         * <pre>
3482         *      1) The curve is checked to see if it is approx. a straight line.
3483         *          1b) If it is, then a line-line intersection is performed, the
3484         *              intersection point added to the output list and the method exited.
3485         *      2) The curve is divided into two halves.
3486         *          2b) Each half is tested to see if it's bounding box is intersected by the line.
3487         *              If it is, then this method is called with that segment recursively.
3488         *              Otherwise, the method is exited.
3489         * </pre>
3490         *
3491         * A number of class variables are used to pass information to this recursion:
3492         * <pre>
3493         *      _tol  The tolerance to use when determining if the curve is locally linear.
3494         *      _L0   A point on the line being intersected with this curve (origin of the line).
3495         *      _Ldir   A vector indicating the direction of the line.
3496         *      _sParams A list used to collect the output subrange intersection point parametric locations.
3497         * </pre>
3498         *
3499         * @param crv The curve segment being checked for intersection with the line.
3500         * @param s0  The parametric position of the start of this curve segment on the
3501         *            master curve.
3502         * @param s1  The parametric position of the end of this curve segment on the
3503         *            master curve.
3504         */
3505        private void divideAndConquerLine(int depth, Curve crv, double s0, double s1) {
3506            StackContext.enter();
3507            try {
3508                //  Check to see if the input curve is within tolerance of a straight line.
3509                Point p0 = crv.getRealPoint(0);
3510                Point p1 = crv.getRealPoint(1);
3511                Point pAvg = p0.plus(p1).divide(2);
3512                Point pm = crv.getRealPoint(0.5);
3513                if (depth > MAXDEPTH || pAvg.distance(pm).isLessThan(_tol)) {
3514                    //  The input curve segment is nearly a straight line.
3515
3516                    //  Intersect the input line and the line approximating the curve.
3517                    Vector<Length> crvDir = p1.minus(p0).toGeomVector();
3518                    MutablePoint sLn = MutablePoint.newInstance(2);
3519                    int numDims = p0.getPhyDimension();
3520                    Unit<Length> unit = p0.getUnit();
3521                    MutablePoint pi1 = MutablePoint.newInstance(numDims, unit);
3522                    MutablePoint pi2 = MutablePoint.newInstance(numDims, unit);
3523                    GeomUtil.lineLineIntersect(_L0, _Ldir, p0, crvDir, _tol, sLn, pi1, pi2);
3524
3525                    if (pi1.distance(pi2).isLessThan(_tol)) {
3526                        //  Add an intersection point to the list.
3527                        double sln = sLn.getValue(1);
3528                        //  Convert from segment par. pos. to full curve pos.
3529                        double s = segmentPos2Parent(sln, s0, s1);
3530                        _sParams.add(s);
3531                    }
3532
3533                } else {
3534                    //  Sub-divide the curve into two segments.
3535                    GeomList<Curve> segments = crv.splitAt(0.5);
3536
3537                    //  Check for possible intersections on the left segment.
3538                    Curve seg = segments.get(0);
3539                    if (GeomUtil.lineBoundsIntersect(_L0, _Ldir, seg, null)) {
3540                        //  May be an intersection.
3541                        double sl = s0;
3542                        double sh = (s0 + s1) / 2.0;
3543
3544                        //  Recurse down to a finer level of detail.
3545                        divideAndConquerLine(depth + 1, seg, sl, sh);
3546                    }
3547
3548                    //  Check for possible intersections on the right segment.
3549                    seg = segments.get(1);
3550                    if (GeomUtil.lineBoundsIntersect(_L0, _Ldir, seg, null)) {
3551                        //  May be an intersection.
3552                        double sl = (s0 + s1) / 2.0;
3553                        double sh = s1;
3554
3555                        //  Recurse down to a finer level of detail.
3556                        divideAndConquerLine(depth + 1, seg, sl, sh);
3557                    }
3558                }
3559
3560            } finally {
3561                StackContext.exit();
3562            }
3563        }
3564
3565        public static void recycle(LineCurveIntersect instance) {
3566            FACTORY.recycle(instance);
3567        }
3568
3569        private LineCurveIntersect() { }
3570
3571        private static final ObjectFactory<LineCurveIntersect> FACTORY = new ObjectFactory<LineCurveIntersect>() {
3572            @Override
3573            protected LineCurveIntersect create() {
3574                return new LineCurveIntersect();
3575            }
3576
3577            @Override
3578            protected void cleanup(LineCurveIntersect obj) {
3579                obj._thisCurve = null;
3580                obj._tol = null;
3581                obj._L0 = null;
3582                obj._Ldir = null;
3583                obj._sParams = null;
3584            }
3585        };
3586
3587    }
3588
3589    /**
3590     * A class used to intersect an line segment and this curve.
3591     */
3592    private static class LineSegCurveIntersect {
3593
3594        private static final int MAXDEPTH = 20;
3595        private Curve _thisCurve;
3596        private Parameter<Length> _tol;
3597        private LineSegment _line;
3598        private Point _p0;
3599        private Point _p1;
3600        private Vector<Length> _Lu;
3601        private GeomList<PointString<SubrangePoint>> _output;
3602        private FastTable<Double> _crvS, _lineS;
3603
3604        @SuppressWarnings("null")
3605        public static LineSegCurveIntersect newInstance(Curve thisCurve, Parameter<Length> tol,
3606                LineSegment line, GeomList<PointString<SubrangePoint>> output) {
3607
3608            if (thisCurve instanceof GeomTransform)
3609                thisCurve = thisCurve.copyToReal();
3610            if (line instanceof GeomTransform)
3611                line = line.copyToReal();
3612            Unit<Length> unit = thisCurve.getUnit();
3613
3614            //  Create an instance of this class and fill in the class variables.
3615            LineSegCurveIntersect o = FACTORY.object();
3616            o._thisCurve = thisCurve;
3617            o._tol = tol.to(unit);
3618            o._line = line.to(unit);
3619            o._p0 = line.getStart().immutable().to(unit);
3620            o._p1 = line.getEnd().immutable().to(unit);
3621            o._Lu = line.getDerivativeVector().immutable().to(unit);
3622            o._output = output;
3623            o._crvS = FastTable.newInstance();
3624            o._lineS = FastTable.newInstance();
3625
3626            return o;
3627        }
3628
3629        /**
3630         * Method that actually carries out the intersection adding the intersection
3631         * points to the list of lists provided in the constructor.
3632         *
3633         * @return A list containing two PointString instances each containing zero or
3634         *         more subrange points, on this curve and the input line segment
3635         *         respectively, made by the intersection of this curve with the specified
3636         *         line segment. If no intersections are found a list of two empty
3637         *         PointStrings are returned.
3638         */
3639        public GeomList<PointString<SubrangePoint>> intersect() {
3640            divideAndConquerLineSeg(1, _thisCurve, 0, 1);
3641
3642            //  Convert stored parametric positions to subrange points on the curves.
3643            int size = _crvS.size();
3644            for (int i = 0; i < size; ++i) {
3645                double s = _crvS.get(i);
3646                _output.get(0).add(SubrangePoint.newInstance(_thisCurve, s));
3647                s = _lineS.get(i);
3648                _output.get(1).add(SubrangePoint.newInstance(_line, s));
3649            }
3650            FastTable.recycle(_crvS);
3651            _crvS = null;
3652            FastTable.recycle(_lineS);
3653            _lineS = null;
3654
3655            return _output;
3656        }
3657
3658        /**
3659         * Uses a recursive "Divide and Conquer" approach to intersecting a line segment
3660         * with a curve. On each call, the following happens:
3661         * <pre>
3662         *      1) The curve is checked to see if it is approx. a straight line.
3663         *          1b) If it is, then a line-line intersection is performed, the
3664         *              intersection point added to the output list and the method exited.
3665         *      2) The curve is divided into two halves.
3666         *          2b) Each half is tested to see if it's bounding box is intersected by the line.
3667         *              If it is, then this method is called with that segment recursively.
3668         *              Otherwise, the method is exited.
3669         * </pre>
3670         *
3671         * A number of class variables are used to pass information to this recursion:
3672         * <pre>
3673         *      _tol  The tolerance to use when determining if the curve is locally linear.
3674         *      _line The line segment instance we are intersecting the curve with.
3675         *      _p0   A point at the start of the line being intersected with this curve.
3676         *      _p1   A point at the end of the line being intersected with this curve.
3677         *      _Lu   A vector indicating the direction of the line.
3678         *      _crvS A list used to collect the output subrange intersection point parametric
3679         *              locations on this curve.
3680         *      _lineS A list used to collect the output subrange intersection point parametric
3681         *              locations on the line segment curve.
3682         * </pre>
3683         *
3684         * @param crv The curve segment being checked for intersection with the line.
3685         * @param s0  The parametric position of the start of this curve segment on the
3686         *            master curve.
3687         * @param s1  The parametric position of the end of this curve segment on the
3688         *            master curve.
3689         */
3690        private void divideAndConquerLineSeg(int depth, Curve crv, double s0, double s1) {
3691            StackContext.enter();
3692            try {
3693                //  Check to see if the input curve is within tolerance of a straight line.
3694                Point p0 = crv.getRealPoint(0);
3695                Point p1 = crv.getRealPoint(1);
3696                Point pAvg = p0.plus(p1).divide(2);
3697                Point pm = crv.getRealPoint(0.5);
3698                if (depth > MAXDEPTH || pAvg.distance(pm).isLessThan(_tol)) {
3699                    //  The input curve segment is nearly a straight line.
3700
3701                    //  Intersect the input line and the line approximating the curve.
3702                    Vector<Length> crvDir = p1.minus(p0).toGeomVector();
3703                    MutablePoint sLn = MutablePoint.newInstance(2);
3704                    int numDims = p0.getPhyDimension();
3705                    Unit<Length> unit = p0.getUnit();
3706                    MutablePoint pi1 = MutablePoint.newInstance(numDims, unit);
3707                    MutablePoint pi2 = MutablePoint.newInstance(numDims, unit);
3708                    GeomUtil.lineLineIntersect(_p0, _Lu, p0, crvDir, _tol, sLn, pi1, pi2);
3709
3710                    //  Make sure the intersection points fall on the line segments.
3711                    double sl_1 = sLn.getValue(0);
3712                    if (sl_1 > 1) {
3713                        pi1.set(_line.getEnd());
3714                        sl_1 = 1;
3715                    } else if (sl_1 < 0) {
3716                        pi1.set(_line.getStart());
3717                        sl_1 = 0;
3718                    }
3719                    double sl_2 = sLn.getValue(1);
3720                    if (sl_2 > 1) {
3721                        pi2.set(crv.getRealPoint(1));
3722                        sl_2 = 1;
3723                    } else if (sl_2 < 0) {
3724                        pi2.set(crv.getRealPoint(0));
3725                        sl_2 = 0;
3726                    }
3727                    Parameter<Length> dist = pi1.distance(pi2);
3728                    if (dist.isLessThan(_tol)) {
3729                        //  Add an intersection point on each curve to the output lists.
3730                        _lineS.add(sl_1);
3731                        //  Convert from segment par. pos. to full curve pos.
3732                        double s = segmentPos2Parent(sl_2, s0, s1);
3733                        _crvS.add(s);
3734                    }
3735
3736                } else {
3737                    //  Sub-divide the curve into two segments.
3738                    GeomList<Curve> segments = crv.splitAt(0.5);
3739
3740                    //  Check for possible intersections on the left segment.
3741                    Curve seg = segments.get(0);
3742                    if (GeomUtil.lineSegBoundsIntersect(_p0, _p1, seg)) {
3743                        //  May be an intersection.
3744                        double sl = s0;
3745                        double sh = (s0 + s1) / 2;
3746
3747                        //  Recurse down to a finer level of detail.
3748                        divideAndConquerLineSeg(depth + 1, seg, sl, sh);
3749                    }
3750
3751                    //  Check for possible intersections on the right segment.
3752                    seg = segments.get(1);
3753                    if (GeomUtil.lineSegBoundsIntersect(_p0, _p1, seg)) {
3754                        //  May be an intersection.
3755                        double sl = (s0 + s1) / 2;
3756                        double sh = s1;
3757
3758                        //  Recurse down to a finer level of detail.
3759                        divideAndConquerLineSeg(depth + 1, seg, sl, sh);
3760                    }
3761                }
3762
3763            } finally {
3764                StackContext.exit();
3765            }
3766        }
3767
3768        public static void recycle(LineSegCurveIntersect instance) {
3769            FACTORY.recycle(instance);
3770        }
3771
3772        private LineSegCurveIntersect() { }
3773
3774        private static final ObjectFactory<LineSegCurveIntersect> FACTORY = new ObjectFactory<LineSegCurveIntersect>() {
3775            @Override
3776            protected LineSegCurveIntersect create() {
3777                return new LineSegCurveIntersect();
3778            }
3779
3780            @Override
3781            protected void cleanup(LineSegCurveIntersect obj) {
3782                obj._thisCurve = null;
3783                obj._tol = null;
3784                obj._line = null;
3785                obj._p0 = null;
3786                obj._p1 = null;
3787                obj._Lu = null;
3788                obj._output = null;
3789                obj._crvS = null;
3790                obj._lineS = null;
3791            }
3792        };
3793    }
3794
3795    /**
3796     * A class used to intersect an line segment and this curve.
3797     */
3798    private static class CurveCurveIntersect {
3799
3800        private static final int MAXDEPTH = 20;
3801        private AbstractCurve _thisCurve;
3802        private Parameter<Length> _tol;
3803        private Curve _thatCurve;
3804        private GeomList<PointString<SubrangePoint>> _output;
3805        private FastTable<Double> _thisS, _thatS;
3806
3807        public static CurveCurveIntersect newInstance(AbstractCurve thisCurve, Parameter<Length> tol,
3808                Curve thatCurve, GeomList<PointString<SubrangePoint>> output) {
3809
3810            //  Make sure the line and curve have the same dimensions.
3811            int numDims = thisCurve.getPhyDimension();
3812            if (thatCurve.getPhyDimension() > numDims)
3813                throw new IllegalArgumentException(MessageFormat.format(
3814                        RESOURCES.getString("sameOrFewerDimErr"), "target curve", "curve"));
3815            else if (thatCurve.getPhyDimension() < numDims)
3816                thatCurve = thatCurve.toDimension(numDims);
3817
3818            if (thisCurve instanceof GeomTransform)
3819                thisCurve = (AbstractCurve)thisCurve.copyToReal();
3820            if (thatCurve instanceof GeomTransform)
3821                thatCurve = thatCurve.copyToReal();
3822
3823            //  Create an instance of this class and fill in the class variables.
3824            CurveCurveIntersect o = FACTORY.object();
3825            o._thisCurve = thisCurve;
3826            o._tol = tol;
3827            o._thatCurve = thatCurve;
3828            o._output = output;
3829            o._thisS = FastTable.newInstance();
3830            o._thatS = FastTable.newInstance();
3831
3832            return o;
3833        }
3834
3835        /**
3836         * Method that actually carries out the intersection adding the intersection
3837         * points to the list of lists provided in the constructor.
3838         *
3839         * @return A list containing two PointString instances each containing zero or
3840         *         more subrange points, on this curve and the input line segment
3841         *         respectively, made by the intersection of this curve with the specified
3842         *         line segment. If no intersections are found a list of two empty
3843         *         PointStrings are returned.
3844         */
3845        public GeomList<PointString<SubrangePoint>> intersect() {
3846            divideAndConquerCurveSeg(1, _thisCurve, 0, 1, _thatCurve, 0, 1);
3847
3848            //  Convert stored parametric positions to subrange points on the curves.
3849            int size = _thisS.size();
3850            for (int i = 0; i < size; ++i) {
3851                double s = _thisS.get(i);
3852                _output.get(0).add(SubrangePoint.newInstance(_thisCurve, s));
3853                s = _thatS.get(i);
3854                _output.get(1).add(SubrangePoint.newInstance(_thatCurve, s));
3855            }
3856            FastTable.recycle(_thisS);
3857            _thisS = null;
3858            FastTable.recycle(_thatS);
3859            _thatS = null;
3860
3861            return _output;
3862        }
3863
3864        /**
3865         * Uses a recursive "Divide and Conquer" approach to intersecting a curve segment
3866         * with another curve. On each call, the following happens:
3867         * <pre>
3868         *      1) The curve is checked to see if it is approx. a straight line.
3869         *          1b) If it is, then a line-line intersection is performed, the
3870         *              intersection point added to the output list and the method exited.
3871         *      2) The curve is divided into two halves.
3872         *          2b) Each half is tested to see if it's bounding box is intersected by the line.
3873         *              If it is, then this method is called with that segment recursively.
3874         *              Otherwise, the method is exited.
3875         * </pre>
3876         *
3877         * A number of class variables are used to pass information to this recursion:
3878         * <pre>
3879         *      _tol  The tolerance to use when determining if the curve is locally linear.
3880         *      _line The line segment instance we are intersecting the curve with.
3881         *      _p0   A point at the start of the line being intersected with this curve.
3882         *      _p1   A point at the end of the line being intersected with this curve.
3883         *      _Lu   A vector indicating the direction of the line.
3884         *      _output A list used to collect the output subrange intersection points.
3885         * </pre>
3886         *
3887         * @param crv1 The 1st curve segment being checked for intersection with "that"
3888         *             curve.
3889         * @param s0_1 The parametric position of the start of the 1st curve segment on
3890         *             the master "this" curve.
3891         * @param s1_1 The parametric position of the end of the 1st curve segment on the
3892         *             master "this" curve.
3893         * @param crv2 The 2nd curve segment being checked for intersection with "this"
3894         *             curve.
3895         * @param s0_2 The parametric position of the start of the 2nd curve segment on
3896         *             the master "that" curve.
3897         * @param s1_2 The parametric position of the end of the 2nd curve segment on the
3898         *             master "that" curve.
3899         */
3900        private void divideAndConquerCurveSeg(int depth, Curve crv1, double s0_1, double s1_1,
3901                Curve crv2, double s0_2, double s1_2) {
3902
3903            StackContext.enter();
3904            try {
3905                //  Check to see if the 1st input curve is within tolerance of a straight line.
3906                Point p0 = crv1.getRealPoint(0);
3907                Point p1 = crv1.getRealPoint(1);
3908                Point pAvg = p0.plus(p1).divide(2);
3909                Point pm = crv1.getRealPoint(0.5);
3910                if (depth > MAXDEPTH || pAvg.distance(pm).isLessThan(_tol)) {
3911                    //  The 1st curve segment is nearly a straight line.
3912                    LineSeg line1 = LineSeg.valueOf(crv1);
3913
3914                    // Intersect the line segment with crv2.
3915                    GeomList<PointString<SubrangePoint>> strs = crv2.intersect(line1, _tol);
3916                    if (!strs.get(0).isEmpty()) {
3917                        //  Intersections were found.
3918                        //  Loop over the intersection points.
3919                        int numPnts = strs.get(0).size();
3920                        for (int i = 0; i < numPnts; ++i) {
3921                            SubrangePoint pnt1 = strs.get(1).get(i);
3922                            SubrangePoint pnt2 = strs.get(0).get(i);
3923
3924                            //  Convert from segment par. pos. to full curve pos.
3925                            double s1 = pnt1.getParPosition().getValue(0);
3926                            double s = segmentPos2Parent(s1, s0_1, s1_1);
3927                            _thisS.add(s);
3928
3929                            double s2 = pnt2.getParPosition().getValue(0);
3930                            s = segmentPos2Parent(s2, s0_2, s1_2);
3931                            _thatS.add(s);
3932                        }
3933                    }
3934
3935                    return;
3936                }
3937
3938                //  Check to see if the 2nd input curve is within tolerance of a straight line.
3939                p0 = crv2.getRealPoint(0);
3940                p1 = crv2.getRealPoint(1);
3941                pAvg = p0.plus(p1).divide(2);
3942                pm = crv2.getRealPoint(0.5);
3943                if (pAvg.distance(pm).isLessThan(_tol)) {
3944                    //  The 2nd curve segment is nearly a straight line.
3945                    LineSeg line2 = LineSeg.valueOf(crv2);
3946
3947                    // Intersect the line segment with crv1.
3948                    GeomList<PointString<SubrangePoint>> strs = crv1.intersect(line2, _tol);
3949                    if (!strs.get(0).isEmpty()) {
3950                        //  Intersections were found.
3951                        //  Loop over the intersection points.
3952                        int numPnts = strs.get(0).size();
3953                        for (int i = 0; i < numPnts; ++i) {
3954                            SubrangePoint pnt1 = strs.get(0).get(i);
3955                            SubrangePoint pnt2 = strs.get(1).get(i);
3956
3957                            //  Convert from segment par. pos. to full curve pos.
3958                            double s1 = pnt1.getParPosition().getValue(0);
3959                            double s = segmentPos2Parent(s1, s0_1, s1_1);
3960                            _thisS.add(s);
3961
3962                            double s2 = pnt2.getParPosition().getValue(0);
3963                            s = segmentPos2Parent(s2, s0_2, s1_2);
3964                            _thatS.add(s);
3965                        }
3966                    }
3967
3968                } else {
3969                    //  Sub-divide the input curves into two segments each.
3970                    GeomList<Curve> crv1Segs = crv1.splitAt(0.5);
3971                    GeomList<Curve> crv2Segs = crv2.splitAt(0.5);
3972
3973                    //  Check for possible intersections on the left segment of curve 1.
3974                    Curve seg1 = crv1Segs.get(0);
3975                    Curve seg2 = crv2Segs.get(0);
3976                    if (GeomUtil.boundsIntersect(seg1, seg2)) {
3977                        //  May be an intersection.
3978                        double sl_1 = s0_1;
3979                        double sh_1 = (s0_1 + s1_1) / 2;
3980                        double sl_2 = s0_2;
3981                        double sh_2 = (s0_2 + s1_2) / 2;
3982
3983                        //  Recurse down to a finer level of detail.
3984                        divideAndConquerCurveSeg(depth + 1, seg1, sl_1, sh_1, seg2, sl_2, sh_2);
3985                    }
3986
3987                    seg2 = crv2Segs.get(1);
3988                    if (GeomUtil.boundsIntersect(seg1, seg2)) {
3989                        //  May be an intersection.
3990                        double sl_1 = s0_1;
3991                        double sh_1 = (s0_1 + s1_1) / 2;
3992                        double sl_2 = (s0_2 + s1_2) / 2;
3993                        double sh_2 = s1_2;
3994
3995                        //  Recurse down to a finer level of detail.
3996                        divideAndConquerCurveSeg(depth + 1, seg1, sl_1, sh_1, seg2, sl_2, sh_2);
3997                    }
3998
3999                    //  Check for possible intersections on the right segment of curve 1.
4000                    seg1 = crv1Segs.get(1);
4001                    seg2 = crv2Segs.get(0);
4002                    if (GeomUtil.boundsIntersect(seg1, seg2)) {
4003                        //  May be an intersection.
4004                        double sl_1 = (s0_1 + s1_1) / 2;
4005                        double sh_1 = s1_1;
4006                        double sl_2 = s0_2;
4007                        double sh_2 = (s0_2 + s1_2) / 2;
4008
4009                        //  Recurse down to a finer level of detail.
4010                        divideAndConquerCurveSeg(depth + 1, seg1, sl_1, sh_1, seg2, sl_2, sh_2);
4011                    }
4012
4013                    seg2 = crv2Segs.get(1);
4014                    if (GeomUtil.boundsIntersect(seg1, seg2)) {
4015                        //  May be an intersection.
4016                        double sl_1 = (s0_1 + s1_1) / 2;
4017                        double sh_1 = s1_1;
4018                        double sl_2 = (s0_2 + s1_2) / 2;
4019                        double sh_2 = s1_2;
4020
4021                        //  Recurse down to a finer level of detail.
4022                        divideAndConquerCurveSeg(depth + 1, seg1, sl_1, sh_1, seg2, sl_2, sh_2);
4023                    }
4024
4025                }
4026
4027            } finally {
4028                StackContext.exit();
4029            }
4030        }
4031
4032        public static void recycle(CurveCurveIntersect instance) {
4033            FACTORY.recycle(instance);
4034        }
4035
4036        private CurveCurveIntersect() { }
4037
4038        private static final ObjectFactory<CurveCurveIntersect> FACTORY = new ObjectFactory<CurveCurveIntersect>() {
4039            @Override
4040            protected CurveCurveIntersect create() {
4041                return new CurveCurveIntersect();
4042            }
4043
4044            @Override
4045            protected void cleanup(CurveCurveIntersect obj) {
4046                obj._thisCurve = null;
4047                obj._tol = null;
4048                obj._output = null;
4049            }
4050        };
4051    }
4052
4053    /**
4054     * An Evaluatable1D function that calculates the signed distance between a point along
4055     * a curve and a point on a surface. This is used by a 1D root finder to drive that
4056     * distance to zero.
4057     */
4058    private static class CurveSrfIntersectFun extends AbstractEvaluatable1D {
4059
4060        private Surface _srf;
4061        private Curve _thisCurve;
4062        public double tol;
4063        public double _nearS, _nearT;
4064        public boolean firstPass;
4065
4066        public static CurveSrfIntersectFun newInstance(Curve thisCurve, Surface surface, double tol) {
4067
4068            if (surface instanceof GeomTransform)
4069                surface = (Surface)surface.copyToReal();
4070            if (thisCurve instanceof GeomTransform)
4071                thisCurve = thisCurve.copyToReal();
4072
4073            CurveSrfIntersectFun o = FACTORY.object();
4074            o._srf = surface;
4075            o._thisCurve = thisCurve;
4076            o._nearS = -1;
4077            o._nearT = 0;
4078            o.tol = tol;
4079            o.firstPass = false;
4080
4081            return o;
4082        }
4083
4084        /**
4085         * Function that calculates the signed distance between a point on an curve, Q(u),
4086         * and the closest point on a parametric surface, P(s,t).
4087         *
4088         * @param u Parametric distance along the curve.
4089         * @return The signed distance between the point on the curve at u and the closest
4090         *         point on the surface.
4091         */
4092        @Override
4093        public double function(double u) throws RootException {
4094            //  Don't compute anything if the 1st pass flag is set.
4095            if (firstPass) {
4096                firstPass = false;
4097                return 0;
4098            }
4099
4100            StackContext.enter();
4101            try {
4102                //  Compute the point on the curve at u.
4103                Point Qu = _thisCurve.getRealPoint(u);
4104
4105                //  Find the closest point on the surface near "Qu".
4106                SubrangePoint Pst;
4107                if (_nearS < 0)
4108                    Pst = _srf.getClosest(Qu, tol);
4109                else
4110                    Pst = _srf.getClosest(Qu, _nearS, _nearT, tol);
4111                GeomPoint st = Pst.getParPosition();
4112                _nearS = st.getValue(0);
4113                _nearT = st.getValue(1);
4114
4115                //  Find the signed distance between Pst and Qu.
4116                Vector<Dimensionless> nhat = _srf.getNormal(_nearS, _nearT);
4117                Vector<Length> dP = Qu.toGeomVector().minus(Pst.toGeomVector());
4118                double dPsign = dP.dot(nhat).getValue();    //  Projection of dP onto nhat.
4119                dPsign = dPsign >= 0 ? 1 : -1;
4120                double d = dP.normValue();                  //  Unsigned distance.
4121                d = d * dPsign;                             //  Signed distance.
4122
4123                return d;
4124
4125            } finally {
4126                StackContext.exit();
4127            }
4128        }
4129
4130        public static void recycle(CurveSrfIntersectFun instance) {
4131            FACTORY.recycle(instance);
4132        }
4133
4134        private CurveSrfIntersectFun() { }
4135
4136        @SuppressWarnings("unchecked")
4137        private static final ObjectFactory<CurveSrfIntersectFun> FACTORY = new ObjectFactory<CurveSrfIntersectFun>() {
4138            @Override
4139            protected CurveSrfIntersectFun create() {
4140                return new CurveSrfIntersectFun();
4141            }
4142
4143            @Override
4144            protected void cleanup(CurveSrfIntersectFun obj) {
4145                obj._srf = null;
4146                obj._thisCurve = null;
4147            }
4148        };
4149    }
4150
4151    /**
4152     * An Evaulatable1D that returns the increment in area subtended by a curve tangent
4153     * vector at the specified parametric location. This is used to find the enclosed area
4154     * of a curve. This is <code>f(s)</code> in:
4155     * <code>A = int_0^1 { f(s) ds } = int_0^1 { ((p(s)-p0) x p'(s)) ds }</code>.
4156     */
4157    private static class EnclAreaEvaluatable extends AbstractEvaluatable1D {
4158
4159        private AbstractCurve _curve;
4160        private Vector<Dimensionless> _nhat;
4161        private Vector<Length> _p0;
4162        private boolean _is2D;
4163        private boolean _isZero;
4164        public int iterations;
4165
4166        @SuppressWarnings("null")
4167        public static EnclAreaEvaluatable newInstance(AbstractCurve curve, GeomPoint p0) {
4168            if (curve instanceof GeomTransform)
4169                curve = (AbstractCurve)curve.copyToReal();
4170            if (p0 instanceof GeomTransform)
4171                p0 = p0.copyToReal();
4172
4173            EnclAreaEvaluatable o = FACTORY.object();
4174            o._curve = curve;
4175            o._p0 = p0.toGeomVector();
4176            o._is2D = curve.getPhyDimension() == 2;
4177            o._isZero = false;
4178
4179            if (!o._is2D) {
4180                if (curve.isLine(Parameter.valueOf(1e-10, SI.METER))) {
4181                    //  For straight line curves, form normal from curve end points.
4182                    Vector<Length> v1 = curve.getRealPoint(0).toGeomVector().minus(o._p0);
4183                    Vector<Length> v2 = curve.getRealPoint(1).toGeomVector().minus(o._p0);
4184                    Vector n = v1.cross(v2);
4185                    if (n.magValue() <= 1e-10) {
4186                        //  The target point is colinear with the line segment.  The enclosed area is always zero.
4187                        o._isZero = true;
4188                    } else
4189                        o._nhat = n.toUnitVector();
4190
4191                } else
4192                    o._nhat = curve.getPlaneNormal(p0);
4193            }
4194            o.iterations = 0;
4195
4196            return o;
4197        }
4198
4199        @Override
4200        public double function(double s) throws RootException {
4201            ++iterations;
4202
4203            //  If the area must be zero, then don't compute anything.
4204            if (_isZero)
4205                return 0;
4206
4207            StackContext.enter();
4208            try {
4209                //  Extract the curve point (p) and velocity (p') at s.
4210                List<Vector<Length>> ders = _curve.getSDerivatives(s, 1);
4211                Vector<Length> p = ders.get(0);
4212                Vector<Length> pp = ders.get(1);
4213
4214                //  Form the p0p vector.
4215                Vector<Length> p0p = p.minus(_p0);
4216
4217                //  Calculate fs = p0p(s) x p'(s)
4218                double fs;
4219                if (_is2D) {
4220                    //  In 2D u x v = det(u, v) = ux*vy - uy*vx
4221                    fs = p0p.getValue(0) * pp.getValue(1) - p0p.getValue(1) * pp.getValue(0);
4222
4223                } else {
4224                    Vector p0pxpp = p0p.cross(pp);
4225
4226                    //  Determine the sign by projecting the p0p x pp vector onto the
4227                    //  plane normal vector for the curve.
4228                    fs = p0pxpp.dot(_nhat).getValue();
4229                }
4230
4231                return fs;
4232
4233            } finally {
4234                StackContext.exit();
4235            }
4236        }
4237
4238        public static void recycle(EnclAreaEvaluatable instance) {
4239            FACTORY.recycle(instance);
4240        }
4241
4242        private EnclAreaEvaluatable() { }
4243
4244        private static final ObjectFactory<EnclAreaEvaluatable> FACTORY = new ObjectFactory<EnclAreaEvaluatable>() {
4245            @Override
4246            protected EnclAreaEvaluatable create() {
4247                return new EnclAreaEvaluatable();
4248            }
4249
4250            @Override
4251            protected void cleanup(EnclAreaEvaluatable obj) {
4252                obj._curve = null;
4253                obj._nhat = null;
4254                obj._p0 = null;
4255            }
4256        };
4257
4258    }
4259
4260}