001/**
002 * NurbsSurface -- Implementation and interface in common to all NURBS surfaces.
003 *
004 * Copyright (C) 2010-2017, by Joseph A. Huwaldt. All rights reserved.
005 *
006 * This library is free software; you can redistribute it and/or modify it under the terms
007 * of the GNU Lesser General Public License as published by the Free Software Foundation;
008 * either version 2.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 Library General Public License for more details.
013 *
014 * You should have received a copy of the GNU Lesser General Public License along with
015 * this program; if not, write to the Free Software Foundation, Inc., 59 Temple Place -
016 * Suite 330, Boston, MA 02111-1307, USA. Or visit: http://www.gnu.org/licenses/lgpl.html
017 *
018 * This source file is based on, but heavily modified from, the LGPL jgeom (Geometry
019 * Library for Java) code by Samuel Gerber, Copyright (C) 2005.
020 */
021package geomss.geom.nurbs;
022
023import geomss.geom.*;
024import jahuwaldt.js.param.Parameter;
025import java.text.MessageFormat;
026import java.util.Collections;
027import java.util.List;
028import static java.util.Objects.isNull;
029import static java.util.Objects.nonNull;
030import static java.util.Objects.requireNonNull;
031import javax.measure.quantity.Dimensionless;
032import javax.measure.quantity.Length;
033import javolution.context.StackContext;
034import javolution.lang.MathLib;
035import javolution.text.Text;
036import javolution.text.TextBuilder;
037import javolution.util.FastMap;
038import javolution.util.FastTable;
039import javolution.xml.XMLFormat;
040import javolution.xml.stream.XMLStreamException;
041import org.jscience.mathematics.number.Float64;
042import org.jscience.mathematics.vector.Float64Vector;
043
044/**
045 * The interface and implementation in common to all NURBS surfaces.
046 *
047 * <p> Modified by: Joseph A. Huwaldt </p>
048 *
049 * @author Joseph A. Huwaldt, Date: June 14, 2010
050 * @version January 31, 2017
051 */
052@SuppressWarnings({"serial", "CloneableImplementsClone"})
053public abstract class NurbsSurface extends AbstractSurface<NurbsSurface> {
054
055    /**
056     * Return a matrix or network of control points for this surface.
057     *
058     * @return the ordered control points
059     */
060    public abstract ControlPointNet getControlPoints();
061
062    /**
063     * Return the control point matrix size in the s-direction (down a column of control
064     * points).
065     *
066     * @return The control point matrix size in the s-direction
067     */
068    public abstract int getNumberOfRows();
069
070    /**
071     * Return the control point matrix size in the t-direction (across the columns of
072     * control points).
073     *
074     * @return The control point matrix size in the t-direction
075     */
076    public abstract int getNumberOfColumns();
077
078    /**
079     * Return the s-direction knot vector of this surface.
080     *
081     * @return The s-knot vector.
082     */
083    public abstract KnotVector getSKnotVector();
084
085    /**
086     * Return the t-direction knot vector of this surface.
087     *
088     * @return The t-knot vector.
089     */
090    public abstract KnotVector getTKnotVector();
091
092    /**
093     * Return the s-degree of the NURBS surface.
094     *
095     * @return s-degree of surface
096     */
097    public abstract int getSDegree();
098
099    /**
100     * Return the t-degree of the NURBS surface.
101     *
102     * @return t-degree of surface
103     */
104    public abstract int getTDegree();
105
106    /**
107     * Return a new surface that is identical to this one but with the transpose of the
108     * parameterization of this surface. The S and T directions will be swapped.
109     *
110     * @return A new surface that is identical to this one but with the transpose of the
111     *         parameterization of this surface.
112     */
113    @Override
114    public NurbsSurface transpose() {
115        //  Transpose the control point network.
116        ControlPointNet oldCPNet = getControlPoints();
117        ControlPointNet cpNet = oldCPNet.transpose();
118
119        //  Reverse and swap the knot vectors.
120        KnotVector kvT = getSKnotVector().reverse();
121        KnotVector kvS = getTKnotVector().reverse();
122
123        //  Construct a new surface.
124        BasicNurbsSurface srf = BasicNurbsSurface.newInstance(cpNet, kvS, kvT);
125        copyState(srf); //  Copy over the super-class state for this surface to the new one.
126
127        return srf;
128    }
129
130    /**
131     * Return a new surface that is identical to this one, but with the S-parameterization
132     * reversed.
133     *
134     * @return A new surface that is identical to this one, but with the
135     *         S-parameterization reversed.
136     * @see #reverseT
137     */
138    @Override
139    public NurbsSurface reverseS() {
140        //  Reverse the control point network rows.
141        ControlPointNet oldCPNet = getControlPoints();
142        ControlPointNet cpNet = oldCPNet.reverseRows();
143
144        //  Reverse the S-knot vector.
145        KnotVector kvS = getSKnotVector().reverse();
146
147        //  Construct a new surface.
148        BasicNurbsSurface srf = BasicNurbsSurface.newInstance(cpNet, kvS, getTKnotVector());
149        copyState(srf); //  Copy over the super-class state for this surface to the new one.
150
151        return srf;
152    }
153
154    /**
155     * Return a new surface that is identical to this one, but with the T-parameterization
156     * reversed.
157     *
158     * @return A new surface that is identical to this one, but with the
159     *         T-parameterization reversed.
160     * @see #reverseS
161     */
162    @Override
163    public NurbsSurface reverseT() {
164        //  Reverse the control point network columns.
165        ControlPointNet oldCPNet = getControlPoints();
166        ControlPointNet cpNet = oldCPNet.reverseColumns();
167
168        //  Reverse the T-knot vector.
169        KnotVector kvT = getTKnotVector().reverse();
170
171        //  Construct a new surface.
172        BasicNurbsSurface srf = BasicNurbsSurface.newInstance(cpNet, getSKnotVector(), kvT);
173        copyState(srf); //  Copy over the super-class state for this surface to the new one.
174
175        return srf;
176    }
177
178    /**
179     * Split this {@link NurbsSurface} at the specified parametric S-position returning a
180     * list containing two new surfaces (a lower surface with smaller S-parametric
181     * positions than "s" and an upper surface with larger S-parametric positions).
182     *
183     * @param s The S-parametric position where this surface should be split (must not be
184     *          0 or 1!).
185     * @return A list containing two surfaces: 0 == the lower surface, 1 == the upper
186     *         surface.
187     */
188    @Override
189    public GeomList<NurbsSurface> splitAtS(double s) {
190        validateParameter(s, 0);
191        if (parNearEnds(s, TOL_ST)) {
192            throw new IllegalArgumentException(
193                    MessageFormat.format(RESOURCES.getString("canNotSplitAtEnds"), "surface"));
194        }
195
196        GeomList<NurbsSurface> output = GeomList.newInstance();     //  Created outside of StackContext.
197        StackContext.enter();
198        try {
199            //  Construct the upper and lower control point networks and S-knot vector
200            //  by splitting each defining S-curve.
201            FastTable<List<ControlPoint>> cpMatrixL = FastTable.newInstance();
202            FastTable<List<ControlPoint>> cpMatrixU = FastTable.newInstance();
203            KnotVector kvSL = getSKnotVector();
204            KnotVector kvSU = kvSL;
205            int numCols = getNumberOfColumns();
206            for (int i = 0; i < numCols; ++i) {
207                BasicNurbsCurve crv = getSCurve(i);
208                GeomList<NurbsCurve> crvList = crv.splitAt(s);
209                cpMatrixL.add(crvList.get(0).getControlPoints());
210                cpMatrixU.add(crvList.get(1).getControlPoints());
211                kvSL = crvList.get(0).getKnotVector();
212                kvSU = crvList.get(1).getKnotVector();
213            }
214
215            //  Create the new control point networks.
216            ControlPointNet cpNetL = ControlPointNet.valueOf(cpMatrixL);
217            ControlPointNet cpNetU = ControlPointNet.valueOf(cpMatrixU);
218
219            //  Construct the new surfaces.
220            KnotVector kvT = getTKnotVector();
221            BasicNurbsSurface srfL = BasicNurbsSurface.newInstance(cpNetL, kvSL, kvT);
222            BasicNurbsSurface srfU = BasicNurbsSurface.newInstance(cpNetU, kvSU, kvT);
223
224            //  Copy surface segments to the outer context.
225            output.add(StackContext.outerCopy(srfL));
226            output.add(StackContext.outerCopy(srfU));
227
228        } finally {
229            StackContext.exit();
230        }
231
232        return output;
233    }
234
235    /**
236     * Split this {@link NurbsSurface} at the specified parametric T-position returning a
237     * list containing two new surfaces (a lower surface with smaller T-parametric
238     * positions than "t" and an upper surface with larger T-parametric positions).
239     *
240     * @param t The T-parametric position where this surface should be split (must not be
241     *          0 or 1!).
242     * @return A list containing two surfaces: 0 == the lower surface, 1 == the upper
243     *         surface.
244     */
245    @Override
246    public GeomList<NurbsSurface> splitAtT(double t) {
247        validateParameter(0, t);
248        if (parNearEnds(t, TOL_ST)) {
249            throw new IllegalArgumentException(
250                    MessageFormat.format(RESOURCES.getString("canNotSplitAtEnds"), "surface"));
251        }
252
253        GeomList<NurbsSurface> output = GeomList.newInstance();     //  Created outside of StackContext.
254        StackContext.enter();
255        try {
256            //  Construct the upper and lower control point networks and T-knot vector
257            //  by splitting each defining S-curve.
258            FastTable<List<ControlPoint>> cpMatrixL = FastTable.newInstance();
259            FastTable<List<ControlPoint>> cpMatrixU = FastTable.newInstance();
260            KnotVector kvTL = getTKnotVector();
261            KnotVector kvTU = kvTL;
262            int numRows = getNumberOfRows();
263            for (int i = 0; i < numRows; ++i) {
264                BasicNurbsCurve crv = getTCurve(i);
265                GeomList<NurbsCurve> crvList = crv.splitAt(t);
266                cpMatrixL.add(crvList.get(0).getControlPoints());
267                cpMatrixU.add(crvList.get(1).getControlPoints());
268                kvTL = crvList.get(0).getKnotVector();
269                kvTU = crvList.get(1).getKnotVector();
270            }
271
272            //  Create the new control point networks.
273            ControlPointNet cpNetL = ControlPointNet.valueOf(cpMatrixL).transpose();
274            ControlPointNet cpNetU = ControlPointNet.valueOf(cpMatrixU).transpose();
275
276            //  Construct the new surfaces.
277            KnotVector kvS = getSKnotVector();
278            BasicNurbsSurface srfL = BasicNurbsSurface.newInstance(cpNetL, kvS, kvTL);
279            BasicNurbsSurface srfU = BasicNurbsSurface.newInstance(cpNetU, kvS, kvTU);
280
281            //  Copy surface segments to the outer context.
282            output.add(StackContext.outerCopy(srfL));
283            output.add(StackContext.outerCopy(srfU));
284
285        } finally {
286            StackContext.exit();
287        }
288
289        return output;
290    }
291
292    /**
293     * Return the T=0 Boundary for this surface as a NURBS curve.
294     *
295     * @return The T=0 Boundary for this surface as a NURBS curve.
296     */
297    @Override
298    public Curve getT0Curve() {
299        return getSCurve(0);
300    }
301
302    /**
303     * Return the T=1 Boundary for this surface as a NURBS curve.
304     *
305     * @return The T=1 Boundary for this surface as a NURBS curve.
306     */
307    @Override
308    public Curve getT1Curve() {
309        Curve crv = getSCurve(getNumberOfColumns() - 1);
310        return crv;
311    }
312
313    /**
314     * Return the S=0 Boundary for this surface as a NURBS curve.
315     *
316     * @return The S=0 Boundary for this surface as a NURBS curve.
317     */
318    @Override
319    public Curve getS0Curve() {
320        return getTCurve(0);
321    }
322
323    /**
324     * Return the S=1 Boundary for this surface as a NURBS curve.
325     *
326     * @return The S=1 Boundary for this surface as a NURBS curve.
327     */
328    @Override
329    public Curve getS1Curve() {
330        Curve crv = getTCurve(getNumberOfRows() - 1);
331        return crv;
332    }
333
334    /**
335     * Return the S-curve with the specified index from this surface as a NURBS curve.
336     *
337     * @param index The T-index for the curve to return from 0 to getNumberOfColumns()-1.
338     * @return The <code>BasicNurbsCurve</code> built from the specified column in the
339     *         ControlPointNet.
340     */
341    public BasicNurbsCurve getSCurve(int index) {
342        ControlPointNet net = getControlPoints();
343        List<ControlPoint> cps = net.getColumn(index);
344        BasicNurbsCurve crv = BasicNurbsCurve.newInstance(cps, getSKnotVector());
345        FastTable.recycle((FastTable)cps);
346        return crv;
347    }
348
349    /**
350     * Return the T-curve with the specified index from this surface as a NURBS curve.
351     *
352     * @param index The S-index for the curve to return from 0 to getNumberOfRows()-1.
353     * @return The <code>BasicNurbsCurve</code> built from the specified row in the
354     *         ControlPointNet.
355     */
356    public BasicNurbsCurve getTCurve(int index) {
357        ControlPointNet net = getControlPoints();
358        List<ControlPoint> cps = net.getRow(index);
359        BasicNurbsCurve crv = BasicNurbsCurve.newInstance(cps, getTKnotVector());
360        FastTable.recycle((FastTable)cps);
361        return crv;
362    }
363
364    /**
365     * Create and return a new {@link BasicNurbsSurface NURBS surface} that is
366     * geometrically identical to this one but with a new knot inserted at the specified
367     * parametric S-location. The parameterization of the new surface is identical to this
368     * one, but the new surface will have a new knot at the specified s-location and a new
369     * set of control points.
370     *
371     * @param s The parametric s-position where the new knot should be inserted from 0 to
372     *          1.
373     * @param r The number of times that a knot should be inserted at the specified
374     *          position.
375     * @return A new NURBS surface that is identical to this one, but with a new knot
376     *         inserted in the specified parametric s-position.
377     */
378    public BasicNurbsSurface insertSKnot(double s, int r) {
379        validateParameter(s,0);
380
381        StackContext.enter();
382        try {
383            KnotVector kv = null;
384            FastTable<List<ControlPoint>> cpNetList = FastTable.newInstance();
385
386            //  Add a knot to each defining column curve.
387            int numCols = getNumberOfColumns();
388            for (int i = 0; i < numCols; ++i) {
389                BasicNurbsCurve c1 = getSCurve(i);
390                BasicNurbsCurve c2 = c1.insertKnot(s, r);
391                kv = c2.getKnotVector();
392                cpNetList.add(c2.getControlPoints());
393            }
394
395            //  Create a new surface from the new knot vector and set of control points.
396            BasicNurbsSurface srf = BasicNurbsSurface.newInstance(ControlPointNet.valueOf(cpNetList), kv, getTKnotVector());
397            return StackContext.outerCopy(srf);
398
399        } finally {
400            StackContext.exit();
401        }
402    }
403
404    /**
405     * Create and return a new {@link BasicNurbsSurface NURBS surface} that is
406     * geometrically identical to this one but with a new knot inserted at the specified
407     * parametric T-location. The parameterization of the new surface is identical to this
408     * one, but the new surface will have a new knot at the specified T-location and a new
409     * set of control points.
410     *
411     * @param t The parametric T-position where the new knot should be inserted from 0 to
412     *          1.
413     * @param r The number of times that a knot should be inserted at the specified
414     *          position.
415     * @return A new NURBS surface that is identical to this one, but with new knots
416     *         inserted in the specified parametric T-position.
417     */
418    public BasicNurbsSurface insertTKnot(double t, int r) {
419        validateParameter(0,t);
420
421        StackContext.enter();
422        try {
423            KnotVector kv = null;
424            FastTable<FastTable<ControlPoint>> cpNetList = FastTable.newInstance();
425            int numCols = getNumberOfColumns() + 1;
426
427            //  Create an empty list of lists of control points.
428            int numRows = getNumberOfRows();
429            for (int i = 0; i < numCols; ++i) {
430                FastTable<ControlPoint> tbl = FastTable.newInstance();
431                tbl.setSize(numRows);
432                cpNetList.add(tbl);
433            }
434
435            //  Add the knot to each defining row curve.
436            for (int j = 0; j < numRows; ++j) {
437                BasicNurbsCurve c1 = getTCurve(j);
438                BasicNurbsCurve c2 = c1.insertKnot(t, r);
439                kv = c2.getKnotVector();
440                List<ControlPoint> cps = c2.getControlPoints();
441                for (int i = 0; i < numCols; ++i) {
442                    cpNetList.get(i).set(j, cps.get(i));
443                }
444            }
445
446            //  Create a new surface from the new knot vector and set of control points.
447            BasicNurbsSurface srf = BasicNurbsSurface.newInstance(ControlPointNet.valueOf(cpNetList), getSKnotVector(), kv);
448            return StackContext.outerCopy(srf);
449
450        } finally {
451            StackContext.exit();
452        }
453    }
454
455    /**
456     * Return a {@link BasicNurbsSurface NURBS surface} that is geometrically identical to
457     * this one but with the specified list of knots merged into it's S-knot vector. The
458     * parameterization of the new surface is identical to this surface, but the new
459     * surface will have any new knots at the specified locations added and a new set of
460     * control points. If the input knot list matches the knot list in the S direction of
461     * this surface, then this surface is returned.
462     *
463     * @param Um The knot vector to merge with this surface's S-knot vector.
464     * @return A NURBS surface that is identical to this one, but with the merger of the
465     *         input knot vector and this surface's S-knot vector.
466     */
467    public NurbsSurface mergeSKnotVector(Float64Vector Um) {
468
469        //  Check to see if any knots need inserting.
470        Float64Vector iKnots = knotsToInsert(getSKnotVector(), Um);
471        if (isNull(iKnots))
472            return this;
473
474        NurbsSurface newSrf = this.refineSKnotVector(iKnots);
475        copyState(newSrf);  //  Copy over the super-class state for this surface to the new one.
476
477        return newSrf;
478    }
479
480    /**
481     * Return a {@link BasicNurbsSurface NURBS surface} that is geometrically identical to
482     * this one but with the specified list of knots merged into it's T-knot vector. The
483     * parameterization of the new surface is identical to this surface, but the new
484     * surface will have any new knots at the specified locations added and a new set of
485     * control points. If the input knot list matches the knot list in the T direction of
486     * this surface, then this surface is returned.
487     *
488     * @param Um The knot vector to merge with this surface's T-knot vector. May not be null.
489     * @return A NURBS surface that is identical to this one, but with the merger of the
490     *         input knot vector and this surface's T knot vector.
491     */
492    public NurbsSurface mergeTKnotVector(Float64Vector Um) {
493
494        //  Check to see if any knots need inserting.
495        Float64Vector iKnots = knotsToInsert(getTKnotVector(), requireNonNull(Um));
496        if (isNull(iKnots))
497            return this;
498
499        NurbsSurface newSrf = this.refineTKnotVector(iKnots);
500        copyState(newSrf);  //  Copy over the super-class state for this surface to the new one.
501
502        return newSrf;
503    }
504
505    private static Float64Vector knotsToInsert(KnotVector U, Float64Vector Um) {
506        StackContext.enter();
507        try {
508            int Usize = U.length();
509            int Umsize = Um.getDimension();
510
511            //  Find the knots that need inserting.
512            FastTable<Float64> Iknots = FastTable.newInstance();
513
514            boolean done = false;
515            int ia = 0, ib = 0;
516            while (!done) {
517                double Umv = Um.getValue(ib);
518                double Uv = U.getValue(ia);
519                if (MathLib.abs(Umv - Uv) <= TOL_ST) {      //  Umv == Uv
520                    ++ib;
521                    ++ia;
522                } else {
523                    Iknots.add(Float64.valueOf(Umv));
524                    ++ib;
525                }
526                done = (ia >= Usize || ib >= Umsize);
527            }
528
529            //  Check to see if any knots need inserting.
530            if (Iknots.isEmpty())
531                return null;
532
533            Float64Vector iKnots = Float64Vector.valueOf(Iknots);
534            return StackContext.outerCopy(iKnots);
535
536        } finally {
537            StackContext.exit();
538        }
539    }
540
541    /**
542     * Create and return a new {@link BasicNurbsSurface NURBS surface} that is
543     * geometrically identical to this one but with the specified list of knots inserted
544     * into it's S-knot vector. The parameterization of the new surface is identical to
545     * this surface, but the new surface will have the new knots at the specified
546     * locations in the S-direction and a new set of control points. This is more
547     * efficient than repeatedly calling insertSKnot() for multiple knot insertions.
548     *
549     * @param newKnots The parametric S-locations where the knots are to be inserted. May
550     *                 not be null.
551     * @return A new NURBS surface that is identical to this one, but with the new knots
552     *         inserted at the specified locations in the parametric S-direction.
553     */
554    public BasicNurbsSurface refineSKnotVector(double[] newKnots) {
555        return refineSKnotVector(Float64Vector.valueOf(requireNonNull(newKnots)));
556    }
557
558    /**
559     * Create and return a new {@link BasicNurbsSurface NURBS surface} that is
560     * geometrically identical to this one but with the specified list of knots inserted
561     * into it's S-knot vector. The parameterization of the new surface is identical to
562     * this surface, but the new surface will have the new knots at the specified
563     * locations in the S-direction and a new set of control points. This is more
564     * efficient than repeatedly calling insertSKnot() for multiple knot insertions.
565     *
566     * @param newKnots The parametric S-locations where the knots are to be inserted. May
567     *                 not be null.
568     * @return A new NURBS surface that is identical to this one, but with the new knots
569     *         inserted at the specified locations in the parametric S-direction.
570     */
571    public BasicNurbsSurface refineSKnotVector(Float64Vector newKnots) {
572
573        //  Check on the validity of the new knots.
574        int r = newKnots.getDimension() - 1;
575        double sOld = 0;
576        for (int i = 0; i <= r; ++i) {
577            double s = newKnots.getValue(i);
578            validateParameter(s,0);
579            if (s < sOld) {
580                throw new IllegalArgumentException(RESOURCES.getString("knotsNotIncreasingErr"));
581            }
582            sOld = s;
583        }
584
585        StackContext.enter();
586        try {
587            KnotVector kv = null;
588            FastTable<List<ControlPoint>> cpNetList = FastTable.newInstance();
589
590            //  Add the new knots to each defining column curve.
591            int numCols = getNumberOfColumns();
592            for (int i = 0; i < numCols; ++i) {
593                BasicNurbsCurve c1 = getSCurve(i);
594                BasicNurbsCurve c2 = c1.refineKnotVector(newKnots);
595                kv = c2.getKnotVector();
596                cpNetList.add(c2.getControlPoints());
597            }
598
599            //  Create a new surface from the new knot vector and set of control points.
600            BasicNurbsSurface srf = BasicNurbsSurface.newInstance(ControlPointNet.valueOf(cpNetList), kv, getTKnotVector());
601            return StackContext.outerCopy(srf);
602
603        } finally {
604            StackContext.exit();
605        }
606    }
607
608    /**
609     * Create and return a new {@link BasicNurbsSurface NURBS surface} that is
610     * geometrically identical to this one but with the specified list of knots inserted
611     * into it's T-knot vector. The parameterization of the new surface is identical to
612     * this surface, but the new surface will have the new knots at the specified
613     * locations in the T-direction and a new set of control points. This is more
614     * efficient than repeatedly calling insertTKnot() for multiple knot insertions.
615     *
616     * @param newKnots The parametric T-locations where the knots are to be inserted. May
617     *                 not be null.
618     * @return A new NURBS surface that is identical to this one, but with the new knots
619     *         inserted at the specified locations in the parametric T-direction.
620     */
621    public BasicNurbsSurface refineTKnotVector(double[] newKnots) {
622        return refineTKnotVector(Float64Vector.valueOf(requireNonNull(newKnots)));
623    }
624
625    /**
626     * Create and return a new {@link BasicNurbsSurface NURBS surface} that is
627     * geometrically identical to this one but with the specified list of knots inserted
628     * into it's T-knot vector. The parameterization of the new surface is identical to
629     * this surface, but the new surface will have the new knots at the specified
630     * locations in the T-direction and a new set of control points. This is more
631     * efficient than repeatedly calling insertTKnot() for multiple knot insertions.
632     *
633     * @param newKnots The parametric T-locations where the knots are to be inserted. May
634     *                 not be null.
635     * @return A new NURBS surface that is identical to this one, but with the new knots
636     *         inserted at the specified locations in the parametric T-direction.
637     */
638    public BasicNurbsSurface refineTKnotVector(Float64Vector newKnots) {
639
640        //  Check on the validity of the new knots.
641        int r = newKnots.getDimension() - 1;
642        double tOld = 0;
643        for (int i = 0; i <= r; ++i) {
644            double t = newKnots.getValue(i);
645            validateParameter(0,t);
646            if (t < tOld) {
647                throw new IllegalArgumentException(RESOURCES.getString("knotsNotIncreasingErr"));
648            }
649            tOld = t;
650        }
651
652        StackContext.enter();
653        try {
654            KnotVector kv = null;
655            FastTable<FastTable<ControlPoint>> cpNetList = FastTable.newInstance();
656            int numCols = getNumberOfColumns() + newKnots.getDimension();
657
658            //  Create an empty list of lists of control points.
659            int numRows = getNumberOfRows();
660            for (int i = 0; i < numCols; ++i) {
661                FastTable<ControlPoint> tbl = FastTable.newInstance();
662                tbl.setSize(numRows);
663                cpNetList.add(tbl);
664            }
665
666            //  Add the knot to each defining row curve.
667            for (int j = 0; j < numRows; ++j) {
668                BasicNurbsCurve c1 = getTCurve(j);
669                BasicNurbsCurve c2 = c1.refineKnotVector(newKnots);
670                kv = c2.getKnotVector();
671                List<ControlPoint> cps = c2.getControlPoints();
672                for (int i = 0; i < numCols; ++i) {
673                    cpNetList.get(i).set(j, cps.get(i));
674                }
675            }
676
677            //  Create a new surface from the new knot vector and set of control points.
678            BasicNurbsSurface srf = BasicNurbsSurface.newInstance(ControlPointNet.valueOf(cpNetList), getSKnotVector(), kv);
679            return StackContext.outerCopy(srf);
680
681        } finally {
682            StackContext.exit();
683        }
684    }
685
686    /**
687     * Attempts to remove the knot in the S-direction with the specified index from this
688     * NURBS surface the specified number of times. If the knot can not be removed without
689     * changing the shape of the surface, then a reference to this surface is returned:
690     * <code>if (srf == srf.removeKnot(...)) { no knots removed }</code>.
691     *
692     * @param index       The index of the knot to be removed (degree + 1 &le; index &lt;
693     *                    knot vector length - degree).
694     * @param numRemovals The number of times that the knot at "index" should be removed
695     *                    (1 &le; num &le; multiplicity of the knot). If numRemovals is
696     *                    &gt; multiplicity of the knot, then numRemovals is set to the
697     *                    multiplicity of the knot.
698     * @param tolerance   The knot(s) specified will be removed as long as the deviation
699     *                    of the surface is everywhere less than this tolerance. May not
700     *                    be null.
701     * @return A new NURBS surface that is identical to this one (to within the specified
702     *         tolerance), but with the specified knot in the S-direction removed or a
703     *         reference to this surface if the knot could not be removed.
704     */
705    public NurbsSurface removeSKnot(int index, int numRemovals, Parameter<Length> tolerance) {
706        requireNonNull(tolerance);
707        
708        KnotVector kv = null;
709        FastTable<List<ControlPoint>> cpNetList = FastTable.newInstance();
710
711        //  Try to remove a knot from each defining column curve.
712        int numCols = getNumberOfColumns();
713        for (int i = 0; i < numCols; ++i) {
714            BasicNurbsCurve c1 = getSCurve(i);
715            BasicNurbsCurve c2 = (BasicNurbsCurve)c1.removeKnot(index, numRemovals, tolerance);
716            if (c2 == c1)
717                break;
718            kv = c2.getKnotVector();
719            cpNetList.add(c2.getControlPoints());
720        }
721
722        //  If we weren't successfull, then just return a reference to this surface.
723        NurbsSurface srf = this;
724
725        //  Were we successfull at removing the knot from all the defining sections?
726        if (cpNetList.size() == numCols) {
727            //  Create a new surface from the new knot vector and set of control points.
728            srf = BasicNurbsSurface.newInstance(ControlPointNet.valueOf(cpNetList), kv, getTKnotVector());
729        }
730
731        return srf;
732    }
733
734    /**
735     * Attempts to remove the knot in the T-direction with the specified index from this
736     * NURBS surface the specified number of times. If the knot can not be removed without
737     * changing the shape of the surface, then a reference to this surface is returned:
738     * <code>if (srf == srf.removeKnot(...)) { no knots removed }</code>.
739     *
740     * @param index       The index of the knot to be removed (degree + 1 &le; index &lt;
741     *                    knot vector length - degree).
742     * @param numRemovals The number of times that the knot at "index" should be removed
743     *                    (1 &le; num &le; multiplicity of the knot). If numRemovals is
744     *                    &gt; multiplicity of the knot, then numRemovals is set to the
745     *                    multiplicity of the knot.
746     * @param tolerance   The knot(s) specified will be removed as long as the deviation
747     *                    of the surface is everywhere less than this tolerance. May not
748     *                    be null.
749     * @return A new NURBS surface that is identical to this one (to within the specified
750     *         tolerance), but with the specified knot in the T-direction removed or a
751     *         reference to this surface if the knot could not be removed.
752     */
753    public NurbsSurface removeTKnot(int index, int numRemovals, Parameter<Length> tolerance) {
754        requireNonNull(tolerance);
755        
756        KnotVector kv = null;
757        FastTable<FastTable<ControlPoint>> cpNetList = FastTable.newInstance();
758        int numCols = getNumberOfColumns() + 1;
759
760        //  Create an empty list of lists of control points.
761        int numRows = getNumberOfRows();
762        for (int i = 0; i < numCols; ++i) {
763            FastTable<ControlPoint> tbl = FastTable.newInstance();
764            tbl.setSize(numRows);
765            cpNetList.add(tbl);
766        }
767
768        //  Try to remove a knot from each defining row curve.
769        boolean successfull = true;
770        for (int j = 0; j < numRows; ++j) {
771            BasicNurbsCurve c1 = getTCurve(j);
772            BasicNurbsCurve c2 = (BasicNurbsCurve)c1.removeKnot(index, numRemovals, tolerance);
773            if (c2 != c1) {
774                kv = c2.getKnotVector();
775                List<ControlPoint> cps = c2.getControlPoints();
776                for (int i = 0; i < numCols; ++i) {
777                    cpNetList.get(i).set(j, cps.get(i));
778                }
779
780            } else {
781                successfull = false;
782                break;
783            }
784
785        }
786
787        //  If we weren't successfull, then just return a reference to this surface.
788        NurbsSurface srf = this;
789
790        //  Were we successfull at removing the knot from all the defining sections?
791        if (successfull) {
792            //  Create a new surface from the new knot vector and set of control points.
793            srf = BasicNurbsSurface.newInstance(ControlPointNet.valueOf(cpNetList), getSKnotVector(), kv);
794        }
795
796        return srf;
797    }
798
799    /**
800     * Guess an initial set of points to use when gridding to tolerance. The initial
801     * points must include each corner of the surface as well as any discontinuities in
802     * the surface. Other than that, they should be distributed as best as possible to
803     * indicate areas of higher curvature, but should be generated as efficiently as
804     * possible.
805     * <p>
806     * This implementation returns the start and end points of each Bezier surface patch.
807     * It also uses the surface boundary curves as guides for inserting additional
808     * interior points.
809     * </p>
810     *
811     * @param tol The maximum distance that a straight line between the final set of
812     *            gridded points may deviate from this surface. May not be null or zero.
813     * @return An array of subrange points to use as a starting point for the grid to
814     *         tolerance algorithms.
815     */
816    @Override
817    protected PointArray<SubrangePoint> guessGrid2TolPoints(Parameter<Length> tol) {
818        requireNonNull(tol);
819        
820        //  Create a list of Bezier segment starting parametric positions ( and the last end position, 1).
821        FastTable<FastTable<Double>> params = getPatchParameters();
822        FastTable<Double> bParamsS = params.get(0);
823        FastTable<Double> bParamsT = params.get(1);
824        FastTable.recycle(params);
825
826        //  If there are only start and end points, insert 0.5 to ensure there is always at least one intermediate point.
827        if (bParamsS.size() < 3) {
828            bParamsS.add(1, 0.5);
829        }
830        if (bParamsT.size() < 3) {
831            bParamsT.add(1, 0.5);
832        }
833
834        if (getSDegree() > 1 || getTDegree() > 1) {
835            StackContext.enter();
836            try {
837                //  Extract edge curves and use their parameterizations as a 1st guess.
838
839                BasicNurbsCurve s0crv = (BasicNurbsCurve)this.getS0Curve();
840                enrichParPosition(bParamsT, s0crv, tol);
841
842                BasicNurbsCurve s1crv = (BasicNurbsCurve)this.getS1Curve();
843                enrichParPosition(bParamsT, s1crv, tol);
844
845                BasicNurbsCurve t0crv = (BasicNurbsCurve)this.getT0Curve();
846                enrichParPosition(bParamsS, t0crv, tol);
847
848                BasicNurbsCurve t1crv = (BasicNurbsCurve)this.getT1Curve();
849                enrichParPosition(bParamsS, t1crv, tol);
850
851            } finally {
852                StackContext.exit();
853            }
854        }
855
856        //  Turn the lists of parametric positions into an array of subrange points.
857        PointArray<SubrangePoint> pntArray = PointArray.newInstance();
858        for (double t : bParamsT) {
859            PointString<SubrangePoint> str = PointString.newInstance();
860            for (double s : bParamsS) {
861                SubrangePoint pnt = SubrangePoint.newInstance(this, Point.valueOf(s, t));
862                str.add(pnt);
863            }
864            pntArray.add(str);
865        }
866
867        FastTable.recycle(bParamsS);
868        FastTable.recycle(bParamsT);
869
870        return pntArray;
871    }
872
873    /**
874     * Enrich an existing list of parametric positions with those extracted from a curve
875     * gridded to the input tolerance.
876     *
877     * @param parPosLst The existing list of parametric positions on a curve.
878     * @param crv       The curve to extract additional parametric positions using the
879     *                  specified tolerance.
880     * @param tol       The tolerance to use when gridding the specified curve.
881     */
882    private static void enrichParPosition(List<Double> parPosLst, Curve crv, Parameter<Length> tol) {
883        //  Grid the curve to tolerance.
884        PointString<SubrangePoint> subPnts = crv.gridToTolerance(tol);
885
886        //  Loop over all the curve points and extract the parametric positions.
887        int numPnts = subPnts.size() - 1;
888        for (int i = 1; i < numPnts; ++i) {
889            SubrangePoint pi = subPnts.get(i);
890            double ti = pi.getParPosition().getValue(0);
891
892            //  Add the curve parametric position to the existing list of parametric positions if it is new.
893            int idx = Collections.binarySearch(parPosLst, ti);  //  Search for the desired value.
894            if (idx < 0) {
895                //  The value is not in the list, so insert it.
896                idx = -(idx + 1);   //  Convert output from binarySearch to insertion point index.
897                parPosLst.add(idx, ti);
898            }
899        }
900    }
901
902    /**
903     * Method that guesses the most likely location for the closest or farthest point on a
904     * surface and returns that location as a 2D point containing the "s" and "t"
905     * parameter values. This is called by getClosest() and getFarthest(). This is
906     * required in order for the root-finding algorithm to reliably refine the closest
907     * point to the correct location.
908     *
909     * @param point   The point to find the closest point on this surface to. May not be null.
910     * @param closest Set to <code>true</code> to return the closest point or
911     *                <code>false</code> to return the farthest point.
912     * @return The 2D parametric point on this surface that is closest to the specified
913     *         point.
914     */
915    @Override
916    protected GeomPoint guessClosestFarthest(GeomPoint point, boolean closest) {
917        requireNonNull(point);
918        //  Inspired by: Ma, Hewitt, "Point inversion and projection for NURBS curve and surface", CAGD, 2003.
919
920        //  Create a list of Bezier segment starting parametric positions ( and the last end position, 1).
921        FastTable<FastTable<Double>> params = getPatchParameters();
922        FastTable<Double> bParamsS = params.get(0);
923        FastTable<Double> bParamsT = params.get(1);
924        FastTable.recycle(params);
925
926        //  Decompose the surface into Bezier patches.
927        GeomList<GeomList<BasicNurbsSurface>> patches = SurfaceUtils.decomposeToBezier(this);
928
929        //  Create a list of potential closest points.
930        FastTable<SubrangePoint> potentials = FastTable.newInstance();
931
932        //  Loop over all the Bezier patches.
933        int column = 0;
934        int ncols = patches.size();
935        for (int i = 0; i < ncols; ++i) {
936            GeomList<BasicNurbsSurface> patchCol = patches.get(i);
937            double t0 = bParamsT.get(column);
938            double tWidth = bParamsT.get(column + 1) - t0;
939
940            int row = 0;
941            int nrows = patchCol.size();
942            for (int j = 0; j < nrows; ++j) {
943                BasicNurbsSurface patch = patchCol.get(j);
944                double s0 = bParamsS.get(row);
945                double sWidth = bParamsS.get(row + 1) - s0;
946
947                //  Determine if this patch is a "valid" patch (simple and convex).
948                ControlPointNet patchCPNet = patch.getControlPoints();
949                if (isValidCPNet(patchCPNet)) {
950                    //  Determine if the boundary curves are closest to the target point.
951                    if (isBezierEdgeCurveNearest(point, patchCPNet)) {
952
953                        //  One of the boundary curves is nearest.
954                        potentials.add(closestFarthestEdgeCurve(patch, s0, t0, sWidth, tWidth, point, closest));
955
956                    } else {
957                        //  An interior point is nearest.
958                        potentials.add(closestFarthestInterior(patch, s0, t0, sWidth, tWidth, point, closest));
959                    }
960
961                } else {
962                    //  Handle invalid patches by just searching for interior points.
963                    potentials.add(closestFarthestInterior(patch, s0, t0, sWidth, tWidth, point, closest));
964                }
965
966                BasicNurbsSurface.recycle(patch);
967                ++row;
968            }
969            GeomList.recycle(patchCol);
970            ++column;
971        }
972        GeomList.recycle(patches);
973        FastTable.recycle(bParamsS);
974        FastTable.recycle(bParamsT);
975
976        //  Find the closest/farthest of all the potential points.
977        SubrangePoint pOpt = closestFarthestInList(potentials, point, closest);
978        FastTable.recycle(potentials);
979
980        //  Extract the parametric position of the closest point.
981        GeomPoint parPos = pOpt.getParPosition();
982
983        //  Return the parametric position of the best guess at the closest/farthest point.
984        return parPos;
985    }
986
987    /**
988     * Method that guesses the most likely location for the closest or farthest point on a
989     * surface and returns that location as a 2D point containing the "s" and "t"
990     * parameter values. This is called by getClosest() and getFarthest(). This is
991     * required in order for the root-finding algorithm to reliably refine the closest
992     * point to the correct location.
993     *
994     * @param plane   The plane to find the closest point on this surface to. May not be null.
995     * @param closest Set to <code>true</code> to return the closest point or
996     *                <code>false</code> to return the farthest point.
997     * @return The 2D parametric point on this surface that is closest to the specified
998     *         plane.
999     */
1000    @Override
1001    protected GeomPoint guessClosestFarthest(GeomPlane plane, boolean closest) {
1002        requireNonNull(plane);
1003        
1004        //  Create a list of Bezier segment starting parametric positions ( and the last end position, 1).
1005        FastTable<FastTable<Double>> params = getPatchParameters();
1006        FastTable<Double> bParamsS = params.get(0);
1007        FastTable<Double> bParamsT = params.get(1);
1008        FastTable.recycle(params);
1009
1010        //  Decompose the surface into Bezier patches.
1011        GeomList<GeomList<BasicNurbsSurface>> patches = SurfaceUtils.decomposeToBezier(this);
1012
1013        //  Create a list of potential closest points.
1014        FastTable<SubrangePoint> potentials = FastTable.newInstance();
1015
1016        //  Loop over all the Bezier patches.
1017        int column = 0;
1018        int ncols = patches.size();
1019        for (int i = 0; i < ncols; ++i) {
1020            GeomList<BasicNurbsSurface> patchCol = patches.get(i);
1021            double t0 = bParamsT.get(column);
1022            double tWidth = bParamsT.get(column + 1) - t0;
1023
1024            int row = 0;
1025            int nrows = patchCol.size();
1026            for (int j = 0; j < nrows; ++j) {
1027                BasicNurbsSurface patch = patchCol.get(j);
1028                double s0 = bParamsS.get(row);
1029                double sWidth = bParamsS.get(row + 1) - s0;
1030
1031                //  Just search for interior points.
1032                potentials.add(closestFarthestInterior(patch, s0, t0, sWidth, tWidth, plane, closest));
1033
1034                BasicNurbsSurface.recycle(patch);
1035                ++row;
1036            }
1037            GeomList.recycle(patchCol);
1038            ++column;
1039        }
1040        GeomList.recycle(patches);
1041        FastTable.recycle(bParamsS);
1042        FastTable.recycle(bParamsT);
1043
1044        //  Find the closest/farthest of all the potential points.
1045        SubrangePoint pOpt = closestFarthestInList(potentials, plane, closest);
1046        FastTable.recycle(potentials);
1047
1048        //  Extract the parametric position of the closest point.
1049        GeomPoint parPos = pOpt.getParPosition();
1050
1051        //  Return the parametric position of the best guess at the closest/farthest point.
1052        return parPos;
1053    }
1054
1055    /**
1056     * Determine if the provided control point network is "valid" (simple and and convex)
1057     * in both the S and T directions.
1058     *
1059     * @param cpNet THe control point network to evaluate.
1060     */
1061    private static boolean isValidCPNet(ControlPointNet cpNet) {
1062        //  Algorithm 2 from Ma, Hewitt, "Point inversion and projection for NURBS curve and surface", CAGD, 2003.
1063
1064        StackContext.enter();
1065        try {
1066            //  Detect every control polygon in S direction
1067            int m = cpNet.getNumberOfColumns();
1068            for (int i = 0; i < m; ++i) {
1069                List<ControlPoint> P = cpNet.getColumn(i);      //  Generate a control polygon
1070
1071                if (!CurveUtils.isSimpleConvexPolygon(P)) {
1072                    //  Control point net is not valid
1073                    return false;
1074                }
1075            }
1076
1077            //  Detect every control polygon in T direction
1078            int n = cpNet.getNumberOfRows();
1079            for (int i = 0; i < n; ++i) {
1080                List<ControlPoint> P = cpNet.getRow(i);         //  Generate a control polygon
1081
1082                if (!CurveUtils.isSimpleConvexPolygon(P)) {
1083                    //  Control point net is not valid
1084                    return false;
1085                }
1086            }
1087
1088            return true;
1089
1090        } finally {
1091            StackContext.exit();
1092        }
1093    }
1094
1095    /**
1096     * Returns <code>true</code> if a boundary curve point is the closest/farthest point
1097     * in the specified Bezier patch to the target point. Returns <code>false</code> if
1098     * the closest/farthest point may be interior to the Bezier patch.
1099     *
1100     * @param point The point to find the closest/farthest point on this surface to/from.
1101     * @param cpNet The control point network for the Bezier patch.
1102     * @return <code>true</code> if a boundary curve point of the specified patch is
1103     *         closest to the target point.
1104     */
1105    private static boolean isBezierEdgeCurveNearest(GeomPoint point, ControlPointNet cpNet) {
1106        //  Algorithm 4 from Ma, Hewitt, "Point inversion and projection for NURBS curve and surface", CAGD, 2003
1107        //  (with logic reversed).
1108        //StackContext.enter();
1109        try {
1110            boolean flag = true;
1111
1112            //  Detect every control polygon in S direction
1113            int m = cpNet.getNumberOfColumns();
1114            for (int i = 0; i < m; ++i) {
1115                List<ControlPoint> P = cpNet.getColumn(i);      //  Generate a control polygon
1116
1117                if (!CurveUtils.isBezierEndPointNearest(point, P)) {
1118                    flag = false;
1119                    break;
1120                }
1121            }
1122            if (flag) {
1123                return true;    //  The nearest point is the point on the boundary curves
1124            }
1125            //  Detect every control polygon in T direction
1126            int n = cpNet.getNumberOfRows();
1127            for (int i = 0; i < n; ++i) {
1128                List<ControlPoint> P = cpNet.getRow(i);         //  Generate a control polygon
1129
1130                if (!CurveUtils.isBezierEndPointNearest(point, P)) {
1131                    flag = false;
1132                    break;
1133                }
1134            }
1135
1136            return flag;
1137
1138        } finally {
1139            //StackContext.exit();
1140        }
1141    }
1142
1143    /**
1144     * Find the closest point among the edge curves of the specified Bezier patch and
1145     * create a subrange point on this surface from that point and return it.
1146     *
1147     * @param patch   The Bezier patch to search the edge curves for.
1148     * @param s0      The parametric S position on this surface of the 0,0 corner of
1149     *                patch.
1150     * @param t0      The parametric T position on this surface of the 0,0 corner of
1151     *                patch.
1152     * @param sWidth  The parametric width of the patch in the S-direction.
1153     * @param tWidth  The parametric width of the patch in the T-direction.
1154     * @param point   The point to find the closest/farthest point on this surface
1155     *                to/from.
1156     * @param closest Set to <code>true</code> to return the closest point or
1157     *                <code>false</code> to return the farthest point.
1158     * @param A       subrange point on this surface corresponding to the closest boundary
1159     *                point on the patch.
1160     */
1161    private SubrangePoint closestFarthestEdgeCurve(BasicNurbsSurface patch,
1162            double s0, double t0, double sWidth, double tWidth, GeomPoint point, boolean closest) {
1163
1164        FastTable<SubrangePoint> potentials = FastTable.newInstance();
1165
1166        if (closest) {
1167            //  Looking for closest point.
1168
1169            //  S=0 side.
1170            Curve c = patch.getS0Curve();
1171            GeomPoint parPos = c.getClosest(point, GTOL).getParPosition();
1172            potentials.add(getPoint(s0, parPos.getValue(0) * tWidth + t0));
1173            BasicNurbsCurve.recycle((BasicNurbsCurve)c);
1174
1175            //  S=1 side.
1176            c = patch.getS1Curve();
1177            parPos = c.getClosest(point, GTOL).getParPosition();
1178            potentials.add(getPoint(sWidth + s0, parPos.getValue(0) * tWidth + t0));
1179            BasicNurbsCurve.recycle((BasicNurbsCurve)c);
1180
1181            //  T=0 side.
1182            c = patch.getT0Curve();
1183            parPos = c.getClosest(point, GTOL).getParPosition();
1184            potentials.add(getPoint(parPos.getValue(0) * sWidth + s0, t0));
1185            BasicNurbsCurve.recycle((BasicNurbsCurve)c);
1186
1187            //  T=1 side.
1188            c = patch.getT1Curve();
1189            parPos = c.getClosest(point, GTOL).getParPosition();
1190            potentials.add(getPoint(parPos.getValue(0) * sWidth + s0, tWidth + t0));
1191            BasicNurbsCurve.recycle((BasicNurbsCurve)c);
1192
1193        } else {
1194            //  Looking for farthest point.
1195
1196            //  S=0 side.
1197            Curve c = patch.getS0Curve();
1198            GeomPoint parPos = c.getFarthest(point, GTOL).getParPosition();
1199            potentials.add(getPoint(s0, parPos.getValue(0) * tWidth + t0));
1200            BasicNurbsCurve.recycle((BasicNurbsCurve)c);
1201
1202            //  S=1 side.
1203            c = patch.getS1Curve();
1204            parPos = c.getFarthest(point, GTOL).getParPosition();
1205            potentials.add(getPoint(sWidth + s0, parPos.getValue(0) * tWidth + t0));
1206            BasicNurbsCurve.recycle((BasicNurbsCurve)c);
1207
1208            //  T=0 side.
1209            c = patch.getT0Curve();
1210            parPos = c.getFarthest(point, GTOL).getParPosition();
1211            potentials.add(getPoint(parPos.getValue(0) * sWidth + s0, t0));
1212            BasicNurbsCurve.recycle((BasicNurbsCurve)c);
1213
1214            //  T=1 side.
1215            c = patch.getT1Curve();
1216            parPos = c.getFarthest(point, GTOL).getParPosition();
1217            potentials.add(getPoint(parPos.getValue(0) * sWidth + s0, tWidth + t0));
1218            BasicNurbsCurve.recycle((BasicNurbsCurve)c);
1219        }
1220
1221        //  Find the closest/farthest of all the potential points.
1222        SubrangePoint pOpt = closestFarthestInList(potentials, point, closest);
1223        FastTable.recycle(potentials);
1224
1225        return pOpt;
1226    }
1227
1228    /**
1229     * Find the closest point on a grid of points that are interior to the specified
1230     * Bezier patch and create a subrange point on this surface from that point and add
1231     * that to the potentials list.
1232     *
1233     * @param patch   The Bezier patch to search.
1234     * @param s0      The parametric S position on this surface of the 0,0 corner of
1235     *                patch.
1236     * @param t0      The parametric T position on this surface of the 0,0 corner of
1237     *                patch.
1238     * @param sWidth  The parametric width of the patch in the S-direction.
1239     * @param tWidth  The parametric width of the patch in the T-direction.
1240     * @param point   The point to find the closest/farthest point on this surface
1241     *                to/from.
1242     * @param closest Set to <code>true</code> to return the closest point or
1243     *                <code>false</code> to return the farthest point.
1244     * @param A       subrange point on this surface corresponding to the closest grid
1245     *                point on the patch.
1246     */
1247    private SubrangePoint closestFarthestInterior(BasicNurbsSurface patch,
1248            double s0, double t0, double sWidth, double tWidth, GeomPoint point, boolean closest) {
1249
1250        //  Find the closest grid point on the patch.
1251        int gridSize = MathLib.max(patch.getSDegree(), patch.getTDegree()) + 2;
1252        SubrangePoint p = getClosestFarthestOnGrid(patch, point, closest, gridSize);
1253
1254        //  Convert to a subrange point on this surface and add to potentials list.
1255        GeomPoint parPos = p.getParPosition();
1256        p = getPoint(parPos.getValue(0) * sWidth + s0, parPos.getValue(1) * tWidth + t0);
1257
1258        return p;
1259    }
1260
1261    /**
1262     * Find the closest point on a grid of points that are interior to the specified
1263     * Bezier patch and create a subrange point on this surface from that point and add
1264     * that to the potentials list.
1265     *
1266     * @param patch   The Bezier patch to search.
1267     * @param s0      The parametric S position on this surface of the 0,0 corner of
1268     *                patch.
1269     * @param t0      The parametric T position on this surface of the 0,0 corner of
1270     *                patch.
1271     * @param sWidth  The parametric width of the patch in the S-direction.
1272     * @param tWidth  The parametric width of the patch in the T-direction.
1273     * @param plane   The plane to find the closest/farthest point on this surface
1274     *                to/from.
1275     * @param closest Set to <code>true</code> to return the closest point or
1276     *                <code>false</code> to return the farthest point.
1277     * @param A       subrange point on this surface corresponding to the closest grid
1278     *                point to the plane on the patch.
1279     */
1280    private SubrangePoint closestFarthestInterior(BasicNurbsSurface patch,
1281            double s0, double t0, double sWidth, double tWidth, GeomPlane plane, boolean closest) {
1282
1283        //  Find the closest grid point on the patch.
1284        int gridSize = MathLib.max(patch.getSDegree(), patch.getTDegree()) + 2;
1285        SubrangePoint p = getClosestFarthestOnGrid(patch, plane, closest, gridSize);
1286
1287        //  Convert to a subrange point on this surface and add to potentials list.
1288        GeomPoint parPos = p.getParPosition();
1289        p = getPoint(parPos.getValue(0) * sWidth + s0, parPos.getValue(1) * tWidth + t0);
1290
1291        return p;
1292    }
1293
1294    /**
1295     * Return the closest/farthest point in the specified list of points to the specified
1296     * point.
1297     *
1298     * @param pointList The list of subrange points to search.
1299     * @param point     The point to find the closest/farthest point on this surface
1300     *                  to/from.
1301     * @param closest   Set to <code>true</code> to return the closest point or
1302     *                  <code>false</code> to return the farthest point.
1303     * @return The closest/farthest point in the list to/from the target point.
1304     */
1305    private static SubrangePoint closestFarthestInList(List<SubrangePoint> pointList, GeomPoint point, boolean closest) {
1306
1307        SubrangePoint pOpt = pointList.get(0);
1308        double dOpt = pOpt.distanceValue(point);
1309        int size = pointList.size();
1310        for (int i = 1; i < size; ++i) {
1311            SubrangePoint p = pointList.get(i);
1312            double dist = p.distanceValue(point);
1313            if (closest) {
1314                if (dist < dOpt) {
1315                    dOpt = dist;
1316                    pOpt = p;
1317                }
1318            } else {
1319                if (dist > dOpt) {
1320                    dOpt = dist;
1321                    pOpt = p;
1322                }
1323            }
1324        }
1325
1326        return pOpt;
1327    }
1328
1329    /**
1330     * Return the closest/farthest point in the specified list of points to the specified
1331     * plane.
1332     *
1333     * @param pointList The list of subrange points to search.
1334     * @param plane     The plane to find the closest/farthest point on this surface
1335     *                  to/from.
1336     * @param closest   Set to <code>true</code> to return the closest point or
1337     *                  <code>false</code> to return the farthest point.
1338     * @return The closest/farthest point in the list to/from the target plane.
1339     */
1340    private static SubrangePoint closestFarthestInList(List<SubrangePoint> pointList, GeomPlane plane, boolean closest) {
1341
1342        SubrangePoint pOpt = pointList.get(0);
1343        double dOpt = pOpt.distanceValue(plane.getClosest(pOpt));
1344        int size = pointList.size();
1345        for (int i = 1; i < size; ++i) {
1346            SubrangePoint p = pointList.get(i);
1347            double dist = p.distanceValue(plane.getClosest(p));
1348            if (closest) {
1349                if (dist < dOpt) {
1350                    dOpt = dist;
1351                    pOpt = p;
1352                }
1353            } else {
1354                if (dist > dOpt) {
1355                    dOpt = dist;
1356                    pOpt = p;
1357                }
1358            }
1359        }
1360
1361        return pOpt;
1362    }
1363
1364    /**
1365     * Returns transformed version of this element. The returned object implements
1366     * {@link geomss.geom.GeomTransform} and contains this element as a child.
1367     *
1368     * @param transform The transformation to apply to this geometry. May not be null.
1369     * @return A new triangle that is identical to this one with the specified
1370     *         transformation applied.
1371     * @throws DimensionException if this point is not 3D.
1372     */
1373    @Override
1374    public NurbsSurfaceTrans getTransformed(GTransform transform) {
1375        return NurbsSurfaceTrans.newInstance(this, requireNonNull(transform));
1376    }
1377
1378    /**
1379     * Return a NURBS surface representation of this surface to within the specified
1380     * tolerance. This implementation returns a reference to this surface itself since it
1381     * is already a NURBS surface and the tolerance parameter is ignored.
1382     *
1383     * @param tol Ignored. May pass <code>null</code>.
1384     * @return A reference to this NURBS surface is returned.
1385     */
1386    @Override
1387    public NurbsSurface toNurbs(Parameter<Length> tol) {
1388        return this;
1389    }
1390
1391    /**
1392     * Return a NURBS surface representation of this surface. An exact representation is
1393     * returned.
1394     *
1395     * @return A reference to this NURBS surface is returned.
1396     */
1397    public NurbsSurface toNurbs() {
1398        return toNurbs(null);
1399    }
1400
1401    /**
1402     * Return <code>true</code> if this NurbsSurface contains valid and finite numerical
1403     * components. A value of <code>false</code> will be returned if any of the control
1404     * point values or knots are NaN or Inf.
1405     *
1406     * @return true if this NurbsSurface contains valid and finite numerical components.
1407     */
1408    @Override
1409    public boolean isValid() {
1410
1411        //  Check the knot vectors.
1412        KnotVector kv = getSKnotVector();
1413        if (!kv.isValid())
1414            return false;
1415        kv = getTKnotVector();
1416        if (!kv.isValid())
1417            return false;
1418
1419        //  Check the control points.
1420        ControlPointNet net = getControlPoints();
1421
1422        return net.isValid();
1423    }
1424
1425    /**
1426     * Return <code>true</code> if this surface is degenerate (i.e.: has area less than
1427     * the specified tolerance squared).
1428     *
1429     * @param tol The tolerance for determining if this surface is degenerate. May not be
1430     *            null.
1431     * @return true if this surface is degenerate.
1432     */
1433    @Override
1434    public boolean isDegenerate(Parameter<Length> tol) {
1435        //  The surface has zero area if it's control point network has zero area.
1436
1437        requireNonNull(tol);
1438        StackContext.enter();
1439        try {
1440            ControlPointNet cpNet = getControlPoints();
1441
1442            int ncols = cpNet.getNumberOfColumns();
1443            for (int j = 0; j < ncols; ++j) {
1444                List<ControlPoint> Pw = cpNet.getColumn(j);
1445                Parameter<Length> length = Parameter.ZERO_LENGTH.to(getUnit());
1446                int size = Pw.size();
1447                Point Pim1 = Pw.get(0).getPoint();
1448                for (int i = 1; i < size; ++i) {
1449                    Point Pi = Pw.get(i).getPoint();
1450                    length = length.plus(Pim1.distance(Pi));
1451                    Pim1 = Pi;
1452                }
1453                //  If any column of the control point network is not zero length,
1454                //  then the surface is not degenerate.
1455                if (length.isGreaterThan(tol))
1456                    return false;
1457            }
1458
1459            return true;
1460
1461        } finally {
1462            StackContext.exit();
1463        }
1464    }
1465
1466    /**
1467     * Return <code>true</code> if this surface is planar or <code>false</code> if it is
1468     * not.
1469     *
1470     * @param tol The geometric tolerance to use in determining if the surface is planar.
1471     *            May not be null.
1472     * @return true if this surface is planar.
1473     */
1474    @Override
1475    public boolean isPlanar(Parameter<Length> tol) {
1476        requireNonNull(tol);
1477        StackContext.enter();
1478        try {
1479            //  If the surface is less than 3D, then it must be planar.
1480            int numDims = getPhyDimension();
1481            if (numDims <= 2)
1482                return true;
1483
1484            //  If the surface is degenerate (zero area; a point or curve), then check edge curves for being planar.
1485            if (isDegenerate(tol)) {
1486                Curve sCrv = getT0Curve();
1487                Curve tCrv = getS0Curve();
1488                if (sCrv.isPlanar(tol) && tCrv.isPlanar(tol)) {
1489                    return true;
1490                }
1491            }
1492
1493            //  Form a plane from the diagonals of the surface's control points.
1494            ControlPointNet cpNet = getControlPoints();
1495            int nSm1 = cpNet.getNumberOfRows() - 1;
1496            int nTm1 = cpNet.getNumberOfColumns() - 1;
1497            Point p0 = cpNet.get(0, 0).getPoint();
1498            Point p1 = cpNet.get(nSm1, nTm1).getPoint();
1499            Vector<Length> d1 = p1.minus(p0).toGeomVector();
1500            p0 = cpNet.get(nSm1, 0).getPoint();
1501            p1 = cpNet.get(0, nTm1).getPoint();
1502            Vector<Length> d2 = p1.minus(p0).toGeomVector();
1503            Vector<Dimensionless> nhat = d1.cross(d2).toUnitVector();
1504            nhat.setOrigin(Point.newInstance(numDims));
1505
1506            //  The curve is planar if the control points are co-planar.
1507            int ncols = cpNet.getNumberOfColumns();
1508            for (int i = 0; i < ncols; ++i) {
1509                List<ControlPoint> col = cpNet.getColumn(i);
1510                int nrows = col.size();
1511                for (int j = 0; j < nrows; ++j) {
1512                    ControlPoint cp = col.get(j);
1513                    Point p = cp.getPoint();
1514                    p1 = GeomUtil.pointPlaneClosest(p, p0, nhat);
1515                    if (p.distance(p1).isGreaterThan(tol))
1516                        return false;
1517                }
1518            }
1519
1520            return true;
1521
1522        } finally {
1523            StackContext.exit();
1524        }
1525    }
1526
1527    /**
1528     * Return a list of lists containing parameters at the start of logical patches of
1529     * this surface in each parametric direction. The first list contains the parameters
1530     * in the "S" direction and the 2nd in the "T" direction. The first element in each
1531     * list must always be 0.0 and the last element 1.0. These should be good places to
1532     * break the surface up into pieces for analysis, but otherwise are arbitrary.
1533     * Subclasses should override this method to provide better patch starting parameters.
1534     *<p>
1535     * This implementation returns the parameters at the start of each Bezier patch of
1536     * this NURBS surface.
1537     *</p>
1538     *
1539     * @return A list of lists containing parameters at the start of logical patches of
1540     *         this surface in each parametric direction.
1541     */
1542    @Override
1543    protected FastTable<FastTable<Double>> getPatchParameters() {
1544        FastTable<FastTable<Double>> lsts = FastTable.newInstance();
1545        lsts.add(SurfaceUtils.getBezierStartParameters(this, 0));
1546        lsts.add(SurfaceUtils.getBezierStartParameters(this, 1));
1547        return lsts;
1548    }
1549
1550    /**
1551     * Return the number of points per patch in the "S" direction that should be used when
1552     * analyzing surface patches by certain methods.
1553     *<p>
1554     * This implementation always returns 2*(p + 1) where p is the degree of the surface
1555     * in the "S" direction.
1556     *</p>
1557     *
1558     * @return The number of points per patch in the "S" direction that should be used
1559     *         when analyzing surface patches.
1560     */
1561    @Override
1562    protected int getNumPointsPerPatchS() {
1563        return 2 * (getSDegree() + 1);
1564    }
1565
1566    /**
1567     * Return the number of points per patch in the "T" direction that should be used when
1568     * analyzing surface patches by certain methods.
1569     *<p>
1570     * This implementation always returns 2*(p + 1) where p is the degree of the surface
1571     * in the "T" direction.
1572     *</p>
1573     *
1574     * @return The number of points per patch in the "T" direction that should be used
1575     *         when analyzing surface patches.
1576     */
1577    @Override
1578    protected int getNumPointsPerPatchT() {
1579        return 2 * (getTDegree() + 1);
1580    }
1581
1582    /**
1583     * Returns the text representation of this geometry element that consists of the name
1584     * followed by the control point values, followed by the knot vectors.
1585     *
1586     * @return the text representation of this geometry element.
1587     */
1588    @Override
1589    public Text toText() {
1590        TextBuilder tmp = TextBuilder.newInstance();
1591        tmp.append('{');
1592        String nameStr = getName();
1593        boolean hasName = nonNull(nameStr);
1594        if (hasName) {
1595            tmp.append(nameStr);
1596            tmp.append(" = {");
1597        }
1598
1599        ControlPointNet cpNet = getControlPoints();
1600        tmp.append(cpNet.toText());
1601        tmp.append(getSKnotVector().toText());
1602        tmp.append(getTKnotVector().toText());
1603
1604        if (hasName) {
1605            tmp.append('}');
1606        }
1607        tmp.append('}');
1608        Text txt = tmp.toText();
1609        TextBuilder.recycle(tmp);
1610        return txt;
1611    }
1612
1613    /**
1614     * Holds the default XML representation for this object.
1615     */
1616    @SuppressWarnings("FieldNameHidesFieldInSuperclass")
1617    protected static final XMLFormat<NurbsSurface> XML = new XMLFormat<NurbsSurface>(NurbsSurface.class) {
1618
1619        @Override
1620        public void read(XMLFormat.InputElement xml, NurbsSurface obj) throws XMLStreamException {
1621            String name = xml.getAttribute("name", (String)null);
1622            if (nonNull(name))
1623                obj.setName(name);
1624            FastMap userData = xml.get("UserData", FastMap.class);
1625            if (nonNull(userData)) {
1626                obj.putAllUserData(userData);
1627                FastMap.recycle(userData);
1628            }
1629        }
1630
1631        @Override
1632        public void write(NurbsSurface obj, XMLFormat.OutputElement xml) throws XMLStreamException {
1633            String name = obj.getName();
1634            if (nonNull(name))
1635                xml.setAttribute("name", name);
1636            FastMap userData = (FastMap)obj.getAllUserData();
1637            if (!userData.isEmpty())
1638                xml.add(userData, "UserData", FastMap.class);
1639            FastMap.recycle(userData);
1640        }
1641    };
1642
1643}