Page MenuHome
Differential D6807 Diff 21646 extern/opennurbs/opennurbs_surface.cpp

Changeset View

Standalone View

View Options

extern/opennurbs/opennurbs_surface.cpp

  • This file was added.
/* $NoKeywords: $ */
/*
//
// Copyright (c) 1993-2012 Robert McNeel & Associates. All rights reserved.
// OpenNURBS, Rhinoceros, and Rhino3D are registered trademarks of Robert
// McNeel & Associates.
//
// THIS SOFTWARE IS PROVIDED "AS IS" WITHOUT EXPRESS OR IMPLIED WARRANTY.
// ALL IMPLIED WARRANTIES OF FITNESS FOR ANY PARTICULAR PURPOSE AND OF
// MERCHANTABILITY ARE HEREBY DISCLAIMED.
//
// For complete openNURBS copyright information see <http://www.opennurbs.org>.
//
////////////////////////////////////////////////////////////////
*/
#include "opennurbs.h"
#if !defined(ON_COMPILING_OPENNURBS)
// This check is included in all opennurbs source .c and .cpp files to insure
// ON_COMPILING_OPENNURBS is defined when opennurbs source is compiled.
// When opennurbs source is being compiled, ON_COMPILING_OPENNURBS is defined
// and the opennurbs .h files alter what is declared and how it is declared.
#error ON_COMPILING_OPENNURBS must be defined when compiling opennurbs
#endif
ON_VIRTUAL_OBJECT_IMPLEMENT(ON_Surface,ON_Geometry,"4ED7D4E1-E947-11d3-BFE5-0010830122F0");
ON_Surface::ON_Surface()
: ON_Geometry()
{}
ON_Surface::ON_Surface(const ON_Surface& src)
: ON_Geometry(src)
{}
unsigned int ON_Surface::SizeOf() const
{
unsigned int sz = ON_Geometry::SizeOf();
sz += (sizeof(*this) - sizeof(ON_Geometry));
// Currently, the size of m_stree is not included
// because this is cached runtime information.
// Applications that care about object size are
// typically storing "inactive" objects for potential
// future use and should call DestroyRuntimeCache(true)
// to remove any runtime cache information.
return sz;
}
ON_Surface& ON_Surface::operator=(const ON_Surface& src)
{
DestroySurfaceTree();
ON_Geometry::operator=(src);
return *this;
}
ON_Surface::~ON_Surface()
{
// Do not call the (virtual) DestroyRuntimeCache or
// DestroySurfaceTree (which calls DestroyRuntimeCache()
// because it opens the potential for crashes in a
// "dirty" destructors of classes derived from ON_Surface
// that to not use DestroyRuntimeCache() in their
// destructors and to not set deleted pointers to zero.
}
ON_Surface* ON_Surface::DuplicateSurface() const
{
return Duplicate();
}
ON::object_type ON_Surface::ObjectType() const
{
return ON::surface_object;
}
bool ON_Surface::GetDomain( int dir, double* t0, double* t1 ) const
{
ON_Interval d = Domain(dir);
if ( t0 ) *t0 = d[0];
if ( t1 ) *t1 = d[1];
return d.IsIncreasing();
}
bool ON_Surface::GetSurfaceSize(
double* width,
double* height
) const
{
if ( width )
*width = 0.0;
if ( height )
*height = 0.0;
return false;
}
bool ON_Surface::SetDomain( int dir, ON_Interval domain )
{
return ( dir >= 0
&& dir <= 1
&& domain.IsIncreasing()
&& SetDomain( dir, domain[0], domain[1] )) ? true : false;
}
bool ON_Surface::SetDomain(
int, // 0 sets first parameter's domain, 1 gets second parameter's domain
double, double
)
{
// TODO make this pure virutual when all the source code is available
return false;
}
//////////
// If t is in the domain of the surface, GetSpanVectorIndex() returns the
// span vector index "i" such that span_vector[i] <= t <= span_vector[i+1].
// The "side" parameter determines which span is selected when t is at the
// end of a span.
//
//virtual
bool ON_Surface::GetSpanVectorIndex(
int dir, // 0 gets first parameter's domain, 1 gets second parameter's domain
double t, // [IN] t = evaluation parameter
int side, // [IN] side 0 = default, -1 = from below, +1 = from above
int* span_vector_i, // [OUT] span vector index
ON_Interval* span_domain // [OUT] domain of the span containing "t"
) const
{
bool rc = false;
int i;
int span_count = SpanCount(dir);
if ( span_count > 0 ) {
double* span_vector = (double*)onmalloc((span_count+1)*sizeof(span_vector[0]));
rc = GetSpanVector( dir, span_vector );
if (rc) {
i = ON_NurbsSpanIndex( 2, span_count, span_vector, t, side, 0 );
if ( i >= 0 && i <= span_count ) {
if ( span_vector_i )
*span_vector_i = i;
if ( span_domain )
span_domain->Set( span_vector[i], span_vector[i+1] );
}
else
rc = false;
}
onfree(span_vector);
}
return rc;
}
bool ON_Surface::GetParameterTolerance( // returns tminus < tplus: parameters tminus <= s <= tplus
int dir,
double t, // t = parameter in domain
double* tminus, // tminus
double* tplus // tplus
) const
{
bool rc = false;
ON_Interval d = Domain( dir );
if ( d.IsIncreasing() )
rc = ON_GetParameterTolerance( d.Min(), d.Max(), t, tminus, tplus );
return rc;
}
ON_Surface::ISO
ON_Surface::IsIsoparametric( const ON_Curve& curve, const ON_Interval* subdomain ) const
{
ISO iso = not_iso;
if ( subdomain )
{
ON_Interval cdom = curve.Domain();
double t0 = cdom.NormalizedParameterAt(subdomain->Min());
double t1 = cdom.NormalizedParameterAt(subdomain->Max());
if ( t0 < t1-ON_SQRT_EPSILON )
{
if ( (t0 > ON_SQRT_EPSILON && t0 < 1.0-ON_SQRT_EPSILON) || (t1 > ON_SQRT_EPSILON && t1 < 1.0-ON_SQRT_EPSILON) )
{
cdom.Intersection(*subdomain);
if ( cdom.IsIncreasing() )
{
ON_NurbsCurve nurbs_curve;
if ( curve.GetNurbForm( nurbs_curve, 0.0,&cdom) )
{
return IsIsoparametric( nurbs_curve, 0 );
}
}
}
}
}
ON_BoundingBox bbox;
double tolerance = 0.0;
const int dim = curve.Dimension();
if ( (dim == 2 || dim==3) && curve.GetBoundingBox(bbox) )
{
iso = IsIsoparametric( bbox );
switch (iso) {
case x_iso:
case W_iso:
case E_iso:
// make sure curve is a (nearly) vertical line
// and weed out vertical scribbles
tolerance = bbox.m_max.x - bbox.m_min.x;
if ( tolerance < ON_ZERO_TOLERANCE && ON_ZERO_TOLERANCE*1024.0 <= (bbox.m_max.y-bbox.m_min.y) )
{
// 26 March 2007 Dale Lear
// If tolerance is tiny, then use ON_ZERO_TOLERANCE
// This fixes cases where iso curves where not getting
// the correct flag because tol=1e-16 and the closest
// point to line had calculation errors of 1e-15.
tolerance = ON_ZERO_TOLERANCE;
}
if ( !curve.IsLinear( tolerance ) )
iso = not_iso;
break;
case y_iso:
case S_iso:
case N_iso:
// make sure curve is a (nearly) horizontal line
// and weed out horizontal scribbles
tolerance = bbox.m_max.y - bbox.m_min.y;
if ( tolerance < ON_ZERO_TOLERANCE && ON_ZERO_TOLERANCE*1024.0 <= (bbox.m_max.x-bbox.m_min.x) )
{
// 26 March 2007 Dale Lear
// If tolerance is tiny, then use ON_ZERO_TOLERANCE
// This fixes cases where iso curves where not getting
// the correct flag because tol=1e-16 and the closest
// point to line had calculation errors of 1e-15.
tolerance = ON_ZERO_TOLERANCE;
}
if ( !curve.IsLinear( tolerance ) )
iso = not_iso;
break;
default:
// nothing here
break;
}
}
return iso;
}
ON_Surface::ISO
ON_Surface::IsIsoparametric( const ON_BoundingBox& bbox ) const
{
ISO iso = not_iso;
if ( bbox.m_min.z == bbox.m_max.z ) {
const double ds = bbox.m_max.x - bbox.m_min.x;
const double dt = bbox.m_max.y - bbox.m_min.y;
double a, b, s0, s1, t0, t1;
ON_Interval d = Domain(0);
s0 = d.Min();
s1 = d.Max();
d = Domain(1);
t0 = d.Min();
t1 = d.Max();
double stol = (s1-s0)/32.0;
double ttol = (t1-t0)/32.0;
if ( s0 < s1 && t0 < t1 && ( ds <= stol || dt <= ttol) )
{
if ( ds*(t1-t0) <= dt*(s1-s0) )
{
// check for s = constant iso
if ( bbox.m_max.x <= s0+stol )
{
// check for west side iso
GetParameterTolerance( 0, s0, &a, &b);
if ( a <= bbox.m_min.x && bbox.m_max.x <= b )
iso = W_iso;
}
else if ( bbox.m_min.x >= s1-stol )
{
// check for east side iso
GetParameterTolerance( 0, s1, &a, &b);
if ( a <= bbox.m_min.x && bbox.m_max.x <= b )
iso = E_iso;
}
if ( iso == not_iso && (s0 < bbox.m_max.x || bbox.m_min.x < s1) )
{
// check for interior "u = constant" iso
GetParameterTolerance( 0, 0.5*(bbox.m_min.x+bbox.m_max.x), &a, &b);
if ( a <= bbox.m_min.x && bbox.m_max.x <= b )
iso = x_iso;
}
}
else
{
// check for t = constant iso
if ( bbox.m_max.y <= t0+ttol )
{
// check for south side iso
GetParameterTolerance( 1, t0, &a, &b);
if ( a < bbox.m_min.y && bbox.m_max.y <= b )
iso = S_iso;
}
else if ( bbox.m_min.y >= t1-ttol )
{
// check for north side iso
GetParameterTolerance( 1, t1, &a, &b);
if ( a < bbox.m_min.y && bbox.m_max.y <= b )
iso = N_iso;
}
if ( iso == not_iso && (t0 < bbox.m_max.x || bbox.m_min.x < t1) )
{
// check for interior "t = constant" iso
GetParameterTolerance( 1, 0.5*(bbox.m_min.y+bbox.m_max.y), &a, &b);
if ( a < bbox.m_min.y && bbox.m_max.y <= b )
iso = y_iso;
}
}
}
}
return iso;
}
bool ON_Surface::IsPlanar( ON_Plane* plane, double tolerance ) const
{
return false;
}
bool
ON_Surface::IsClosed(int dir) const
{
ON_Interval d = Domain(dir);
if ( d.IsIncreasing() && Dimension() <= 3 ) {
const int span_count = SpanCount(dir?0:1);
const int span_degree = Degree(dir?0:1);
if ( span_count > 0 && span_degree > 0 )
{
ON_SimpleArray<double> s(span_count+1);
s.SetCount(span_count+1);
int n = 2*span_degree+1;
double delta = 1.0/n;
ON_3dPoint P(ON_3dPoint::Origin);
ON_3dPoint Q(ON_3dPoint::Origin);
int hintP[2] = {0,0};
int hintQ[2] = {0,0};
double *u0, *u1, *v0, *v1;
double t;
ON_Interval sp;
if ( dir ) {
v0 = &d.m_t[0];
v1 = &d.m_t[1];
u0 = &t;
u1 = &t;
}
else {
u0 = &d.m_t[0];
u1 = &d.m_t[1];
v0 = &t;
v1 = &t;
}
if ( GetSpanVector( dir?0:1, s.Array() ) )
{
int span_index, i;
for ( span_index = 0; span_index < span_count; span_index++ ) {
sp.Set(s[span_index],s[span_index+1]);
for ( i = 0; i < n; i++ ) {
t = sp.ParameterAt(i*delta);
if ( !Evaluate( *u0, *v0, 1, 3, P, 0, hintP ) )
return false;
if ( !Evaluate( *u1, *v1, 2, 3, Q, 0, hintQ ) )
return false;
if ( false == ON_PointsAreCoincident( 3, 0, &P.x, &Q.x ) )
return false;
}
}
return true;
}
}
}
return false;
}
bool ON_Surface::IsPeriodic(int dir) const
{
return false;
}
bool ON_Surface::GetNextDiscontinuity(
int dir,
ON::continuity c,
double t0,
double t1,
double* t,
int* hint,
int* dtype,
double cos_angle_tolerance,
double curvature_tolerance
) const
{
// 28 Jan 2005 - untested code
// this function must be overridden by surface objects that
// can have parametric discontinuities on the interior of the curve.
bool rc = false;
// using tmp_dtype means I don't have to keep checking for a null dtype;
int tmp_dtype = 0;
if ( !dtype )
dtype = &tmp_dtype;
*dtype = 0;
if ( t0 != t1 )
{
bool bTestC0 = false;
bool bTestD1 = false;
bool bTestD2 = false;
bool bTestT = false;
bool bTestK = false;
switch(c)
{
case ON::continuity::C0_locus_continuous:
bTestC0 = true;
break;
case ON::continuity::C1_locus_continuous:
bTestC0 = true;
bTestD1 = true;
break;
case ON::continuity::C2_locus_continuous:
bTestC0 = true;
bTestD1 = true;
bTestD2 = true;
break;
case ON::continuity::G1_locus_continuous:
bTestC0 = true;
bTestT = true;
break;
case ON::continuity::G2_locus_continuous:
bTestC0 = true;
bTestT = true;
bTestK = true;
break;
default:
// other values ignored on purpose.
break;
}
if ( bTestC0 )
{
// 20 March 2003 Dale Lear:
// Have to look for locus discontinuities at ends.
// Must test both ends becuase t0 > t1 is valid input.
// In particular, for ON_CurveProxy::GetNextDiscontinuity()
// to work correctly on reversed "real" curves, the
// t0 > t1 must work right.
int hinta[2], hintb[2], span_index, j;
ON_Interval domain = Domain(dir);
ON_Interval span_domain;
ON_2dPoint st0, st1;
ON_3dVector Da[6], Db[6];
ON_3dVector& D1a = Da[1+dir];
ON_3dVector& D1b = Db[1+dir];
ON_3dVector& D2a = Da[3+2*dir];
ON_3dVector& D2b = Db[3+2*dir];
if ( t0 < domain[1] && t1 >= domain[1] )
t1 = domain[1];
else if ( t0 > domain[0] && t1 <= domain[0] )
t1 = domain[0];
if ( (t0 < domain[1] && t1 >= domain[1]) || (t0 > domain[0] && t1 <= domain[0]) )
{
if ( IsClosed(dir) )
{
int span_count = SpanCount(1-dir);
double* span_vector = (span_count>0) ? ((double*)onmalloc((span_count+1)*sizeof(*span_vector))) : 0;
if (!GetSpanVector(1-dir,span_vector))
span_count = 0;
st0[dir] = domain[0];
st1[dir] = domain[1];
for ( span_index = 0; span_index < span_count && 1 != *dtype; span_index++ )
{
span_domain.Set(span_vector[span_index],span_vector[span_index+1]);
for ( j = (span_index?1:0); j <= 2 && 1 != *dtype; j++ )
{
st0[1-dir] = span_domain.ParameterAt(0.5*j);
st1[1-dir] = st0[1-dir];
if ( bTestD1 || bTestT )
{
// need to check locus continuity at start/end of closed surface.
if ( Evaluate(st0.x,st0.y,2,3,&Da[0].x,1,hinta)
&& Evaluate(st1.x,st1.y,2,3,&Db[0].x,2,hintb) )
{
if ( bTestD1 )
{
if ( !(D1a-D1b).IsTiny(D1b.MaximumCoordinate()*ON_SQRT_EPSILON ) )
{
if ( dtype )
*dtype = 1;
*t = t1;
rc = true;
}
else if ( bTestD2 && !(D2a-D2b).IsTiny(D2b.MaximumCoordinate()*ON_SQRT_EPSILON) )
{
if ( dtype )
*dtype = 2;
*t = t1;
rc = true;
}
}
else if ( bTestT )
{
ON_3dVector Ta, Tb, Ka, Kb;
ON_EvCurvature( D1a, D2a, Ta, Ka );
ON_EvCurvature( D1b, D2b, Tb, Kb );
if ( Ta*Tb < cos_angle_tolerance )
{
if ( dtype )
*dtype = 1;
*t = t1;
rc = true;
}
else if ( bTestK && (Ka-Kb).Length() > curvature_tolerance )
{
if ( dtype )
*dtype = 2;
*t = t1;
rc = true;
}
}
}
}
}
}
if ( span_vector)
{
onfree(span_vector);
}
}
else
{
// open curves are not locus continuous at ends.
*dtype = 0; // locus C0 discontinuity
*t = t1;
rc = true;
}
}
}
}
return rc;
}
void ON_IsG1Closed(const ON_Surface& Srf, bool closed[2])
{
ON_Interval dom[2];
double t;
dom[0] = Srf.Domain(0);
dom[1] = Srf.Domain(1);
closed[0] = Srf.IsClosed(0) && !Srf.GetNextDiscontinuity(0, ON::continuity::G1_locus_continuous, dom[0][0], dom[0][1], &t);
closed[1] = Srf.IsClosed(1) && !Srf.GetNextDiscontinuity(1, ON::continuity::G1_locus_continuous, dom[1][0], dom[1][1], &t);
}
static
bool PrincipalCurvaturesAreContinuous(
bool bSmoothTest,
double k1a, double k2a, // side "a" principal curvatures
double k1b, double k2b, // side "b" principal curvatures
double curvature_tolerance
)
{
ON_3dVector K[2];
K[0].x = k1a;
K[0].y = 0.0;
K[0].z = 0.0;
K[1].x = k1b;
K[1].y = 0.0;
K[1].z = 0.0;
// compare the first principal curvatures
bool rc = ( bSmoothTest )
? ON_IsGsmoothCurvatureContinuous(K[0],K[1],0.0,curvature_tolerance)
: ON_IsG2CurvatureContinuous(K[0],K[1],0.0,curvature_tolerance);
if ( rc )
{
// compare the second principal curvatures
K[0].x = k2a;
K[1].x = k2b;
rc = ( bSmoothTest )
? ON_IsGsmoothCurvatureContinuous(K[0],K[1],0.0,curvature_tolerance)
: ON_IsG2CurvatureContinuous(K[0],K[1],0.0,curvature_tolerance);
}
return rc;
}
bool ON_Surface::IsContinuous(
ON::continuity desired_continuity,
double s,
double t,
int* hint, // default = nullptr,
double point_tolerance, // default=ON_ZERO_TOLERANCE
double d1_tolerance, // default==ON_ZERO_TOLERANCE
double d2_tolerance, // default==ON_ZERO_TOLERANCE
double cos_angle_tolerance, // default==ON_DEFAULT_ANGLE_TOLERANCE_COSINE
double curvature_tolerance // default==ON_SQRT_EPSILON
) const
{
int qi, span_count[2];
span_count[0] = SpanCount(0);
span_count[1] = SpanCount(1);
if ( span_count[0] <= 1 && span_count[1] <= 1 )
return true;
ON_3dPoint P[4];
ON_3dVector Ds[4], Dt[4], Dss[4], Dst[4], Dtt[4], N[4], K1[4], K2[4];
double gauss[4], mean[4], kappa1[4], kappa2[4], sq[4], tq[4];
switch ( desired_continuity )
{
case ON::continuity::C0_locus_continuous:
case ON::continuity::C1_locus_continuous:
case ON::continuity::C2_locus_continuous:
case ON::continuity::G1_locus_continuous:
case ON::continuity::G2_locus_continuous:
{
// 7 April 2005 Dale Lear
// This locus continuity test was added. Prior to
// this time, this function ignored the word "locus".
// The reason for the change is that Chuck's filleting code
// needs to query the continuity at the seams of closed surfaces.
// See ON::continuity comments. The different sq[] values
// must NOT be used when s == Domain(0)[0] and must always
// be used when s == Domain(0)[1]. In particular, if a surface
// is not closed in the "s" direction and s == Domain(1)[1], then
// the answer to any locus query is false.
ON_Interval d = Domain(0);
if ( s == d[1] )
{
sq[0] = sq[1] = d[0];
sq[2] = sq[3] = d[1];
}
else
{
sq[0] = sq[1] = sq[2] = sq[3] = s;
}
d = Domain(1);
// See ON::continuity comments. The different tq[] values
// must NOT be used when t == Domain(1)[0] and must always
// be used when t == Domain(1)[1]. In particular, if a surface
// is not closed in the "t" direction and t == Domain(1)[1], then
// the answer to any locus query is false.
if ( t == d[1] )
{
tq[0] = tq[3] = d[0];
tq[1] = tq[2] = d[1];
}
else
{
tq[0] = tq[1] = tq[2] = tq[3] = t;
}
}
break;
default:
sq[0] = sq[1] = sq[2] = sq[3] = s;
tq[0] = tq[1] = tq[2] = tq[3] = t;
break;
}
desired_continuity = ON::ParametricContinuity((int)desired_continuity);
// this is slow and uses evaluation
// virtual overrides on curve classes that can have multiple spans
// are much faster because the avoid evaluation
switch ( desired_continuity )
{
case ON::continuity::C0_continuous:
for ( qi = 0; qi < 4; qi++ )
{
if ( !EvPoint( sq[qi], tq[qi], P[qi], qi+1 ) )
return false;
if ( qi )
{
if ( !(P[qi]-P[qi-1]).IsTiny(point_tolerance) )
return false;
}
}
if ( !(P[3]-P[0]).IsTiny(point_tolerance) )
return false;
break;
case ON::continuity::C1_continuous:
for ( qi = 0; qi < 4; qi++ )
{
if ( !Ev1Der( sq[qi], tq[qi], P[qi], Ds[qi], Dt[qi], qi+1, hint ) )
return false;
if ( qi )
{
if ( !(P[qi]-P[qi-1]).IsTiny(point_tolerance) )
return false;
if ( !(Ds[qi]-Ds[qi-1]).IsTiny(d1_tolerance) )
return false;
if ( !(Dt[qi]-Dt[qi-1]).IsTiny(d1_tolerance) )
return false;
}
}
if ( !(P[3]-P[0]).IsTiny(point_tolerance) )
return false;
if ( !(Ds[3]-Ds[0]).IsTiny(d1_tolerance) )
return false;
if ( !(Dt[3]-Dt[0]).IsTiny(d1_tolerance) )
return false;
break;
case ON::continuity::C2_continuous:
for ( qi = 0; qi < 4; qi++ )
{
if ( !Ev2Der( sq[qi], tq[qi], P[qi], Ds[qi], Dt[qi],
Dss[qi], Dst[qi], Dtt[qi],
qi+1, hint ) )
return false;
if ( qi )
{
if ( !(P[qi]-P[qi-1]).IsTiny(point_tolerance) )
return false;
if ( !(Ds[qi]-Ds[qi-1]).IsTiny(d1_tolerance) )
return false;
if ( !(Dt[qi]-Dt[qi-1]).IsTiny(d1_tolerance) )
return false;
if ( !(Dss[qi]-Dss[qi-1]).IsTiny(d2_tolerance) )
return false;
if ( !(Dst[qi]-Dst[qi-1]).IsTiny(d2_tolerance) )
return false;
if ( !(Dtt[qi]-Dtt[qi-1]).IsTiny(d2_tolerance) )
return false;
}
}
if ( !(P[3]-P[0]).IsTiny(point_tolerance) )
return false;
if ( !(Ds[3]-Ds[0]).IsTiny(d1_tolerance) )
return false;
if ( !(Dt[3]-Dt[0]).IsTiny(d1_tolerance) )
return false;
if ( !(Dss[3]-Dss[0]).IsTiny(d2_tolerance) )
return false;
if ( !(Dst[3]-Dst[0]).IsTiny(d2_tolerance) )
return false;
if ( !(Dtt[3]-Dtt[0]).IsTiny(d2_tolerance) )
return false;
break;
case ON::continuity::G1_continuous:
for ( qi = 0; qi < 4; qi++ )
{
if ( !EvNormal( sq[qi], tq[qi], P[qi], N[qi], qi+1 ) )
return false;
if ( qi )
{
if ( !(P[qi]-P[qi-1]).IsTiny(point_tolerance) )
return false;
if ( N[qi]*N[qi-1] < cos_angle_tolerance )
return false;
}
}
if ( !(P[3]-P[0]).IsTiny(point_tolerance) )
return false;
if ( N[3]*N[0] < cos_angle_tolerance )
return false;
break;
case ON::continuity::G2_continuous:
case ON::continuity::Gsmooth_continuous:
{
bool bSmoothCon = (ON::continuity::Gsmooth_continuous == desired_continuity);
for ( qi = 0; qi < 4; qi++ )
{
if ( !Ev2Der( sq[qi], tq[qi], P[qi], Ds[qi], Dt[qi],
Dss[qi], Dst[qi], Dtt[qi],
qi+1, hint ) )
return false;
if (!ON_EvNormal(qi+1, Ds[qi], Dt[qi], Dss[qi], Dst[qi], Dtt[qi], N[qi]))
return false;
ON_EvPrincipalCurvatures( Ds[qi], Dt[qi], Dss[qi], Dst[qi], Dtt[qi], N[qi],
&gauss[qi], &mean[qi], &kappa1[qi], &kappa2[qi],
K1[qi], K2[qi] );
if ( qi )
{
if ( !(P[qi]-P[qi-1]).IsTiny(point_tolerance) )
return false;
if ( N[qi]*N[qi-1] < cos_angle_tolerance )
return false;
if ( !PrincipalCurvaturesAreContinuous(bSmoothCon,kappa1[qi],kappa2[qi],kappa1[qi-1],kappa2[qi-1],curvature_tolerance) )
return false;
}
}
if ( !(P[3]-P[0]).IsTiny(point_tolerance) )
return false;
if ( N[3]*N[0] < cos_angle_tolerance )
return false;
if ( !PrincipalCurvaturesAreContinuous(bSmoothCon,kappa1[3],kappa2[3],kappa1[0],kappa2[0],curvature_tolerance) )
return false;
}
break;
default:
// intentionally ignoring other ON::continuity enum values
break;
}
return true;
}
bool
ON_Surface::IsSingular(int side) const
{
return false;
}
bool ON_Surface::IsSolid() const
{
const bool bIsClosed0 = ( IsClosed(0) || ( IsSingular(1) && IsSingular(3) ) );
const bool bIsClosed1 = ( IsClosed(1) || ( IsSingular(0) && IsSingular(2) ) );
if ( bIsClosed0 && bIsClosed1 )
return true;
const ON_Extrusion* extrusion = ON_Extrusion::Cast(this);
if ( 0 != extrusion && extrusion->IsSolid() )
return true;
return false;
}
bool
ON_Surface::IsAtSingularity(double s, double t,
bool bExact //true by default
) const
{
if (bExact){
if (s == Domain(0)[0]){
if (IsSingular(3))
return true;
}
else if (s == Domain(0)[1]){
if (IsSingular(1))
return true;
}
if (t == Domain(1)[0]){
if (IsSingular(0))
return true;
}
else if (t == Domain(1)[1]){
if (IsSingular(2))
return true;
}
return false;
}
if (IsAtSingularity(s, t, true))
return true;
bool bCheckPartials[2] = {false, false};
int i;
double m[2];
for (i=0; i<2; i++)
m[i] = Domain(i).Mid();
if (s < m[0]){
if (IsSingular(3))
bCheckPartials[1] = true;
}
else {
if (IsSingular(1))
bCheckPartials[1] = true;
}
if (!bCheckPartials[0] && !bCheckPartials[1]){
if (t < m[1]){
if (IsSingular(0))
bCheckPartials[0] = true;
}
else {
if (IsSingular(2))
bCheckPartials[0] = true;
}
}
if (!bCheckPartials[0] && !bCheckPartials[1])
return false;
ON_3dPoint P;
ON_3dVector M[2], S[2];
if (!Ev1Der(s, t, P, S[0], S[1]))
return false;
if (!Ev1Der(m[0], m[1], P, M[0], M[1]))
return false;
for (i=0; i<2; i++){
if (!bCheckPartials[i])
continue;
if (S[i].Length() < 1.0e-6 * M[i].Length())
return true;
}
return false;
}
int ON_Surface::IsAtSeam(double s, double t) const
{
int rc = 0;
int i;
for (i=0; i<2; i++){
if (!IsClosed(i))
continue;
double p = (i) ? t : s;
if (p == Domain(i)[0] || p == Domain(i)[1])
rc += (i+1);
}
return rc;
}
ON_3dPoint
ON_Surface::PointAt( double s, double t ) const
{
ON_3dPoint p(0.0,0.0,0.0);
EvPoint( s, t, p );
return p;
}
ON_3dVector
ON_Surface::NormalAt( double s, double t ) const
{
ON_3dVector N(0.0,0.0,0.0);
EvNormal( s, t, N );
return N;
}
bool ON_Surface::FrameAt( double u, double v, ON_Plane& frame) const
{
bool rc = false;
ON_3dPoint origin;
ON_3dVector udir, vdir, normal;
if( EvNormal( u, v, origin, udir, vdir, normal))
{
if ( udir.Unitize() )
vdir = ON_CrossProduct( normal, udir);
else if ( vdir.Unitize() )
udir = ON_CrossProduct( vdir, normal);
frame.CreateFromFrame( origin, udir, vdir);
rc = frame.IsValid();
}
return rc;
}
bool
ON_Surface::EvPoint( // returns false if unable to evaluate
double s, double t, // evaluation parameters
ON_3dPoint& point,
int side, // optional - determines which side to evaluate from
// 0 = default
// 1 from NE quadrant
// 2 from NW quadrant
// 3 from SW quadrant
// 4 from SE quadrant
int* hint // optional - evaluation hint (int[2]) used to speed
// repeated evaluations
) const
{
bool rc = false;
double ws[128];
double* v;
if ( Dimension() <= 3 ) {
v = &point.x;
point.x = 0.0;
point.y = 0.0;
point.z = 0.0;
}
else if ( Dimension() <= 128 ) {
v = ws;
}
else {
v = (double*)onmalloc(Dimension()*sizeof(*v));
}
rc = Evaluate( s, t, 0, Dimension(), v, side, hint );
if ( Dimension() > 3 ) {
point.x = v[0];
point.y = v[1];
point.z = v[2];
if ( Dimension() > 128 )
onfree(v);
}
return rc;
}
bool
ON_Surface::Ev1Der( // returns false if unable to evaluate
double s, double t, // evaluation parameters
ON_3dPoint& point,
ON_3dVector& ds,
ON_3dVector& dt,
int side, // optional - determines which side to evaluate from
// 0 = default
// 1 from NE quadrant
// 2 from NW quadrant
// 3 from SW quadrant
// 4 from SE quadrant
int* hint // optional - evaluation hint (int[2]) used to speed
// repeated evaluations
) const
{
bool rc = false;
const int dim = Dimension();
double ws[3*32];
double* v;
point.x = 0.0;
point.y = 0.0;
point.z = 0.0;
ds.x = 0.0;
ds.y = 0.0;
ds.z = 0.0;
dt.x = 0.0;
dt.y = 0.0;
dt.z = 0.0;
if ( dim <= 32 ) {
v = ws;
}
else {
v = (double*)onmalloc(3*dim*sizeof(*v));
}
rc = Evaluate( s, t, 1, dim, v, side, hint );
point.x = v[0];
ds.x = v[dim];
dt.x = v[2*dim];
if ( dim > 1 ) {
point.y = v[1];
ds.y = v[dim+1];
dt.y = v[2*dim+1];
if ( dim > 2 ) {
point.z = v[2];
ds.z = v[dim+2];
dt.z = v[2*dim+2];
if ( dim > 32 )
onfree(v);
}
}
return rc;
}
bool
ON_Surface::Ev2Der( // returns false if unable to evaluate
double s, double t, // evaluation parameters
ON_3dPoint& point,
ON_3dVector& ds,
ON_3dVector& dt,
ON_3dVector& dss,
ON_3dVector& dst,
ON_3dVector& dtt,
int side, // optional - determines which side to evaluate from
// 0 = default
// 1 from NE quadrant
// 2 from NW quadrant
// 3 from SW quadrant
// 4 from SE quadrant
int* hint // optional - evaluation hint (int[2]) used to speed
// repeated evaluations
) const
{
bool rc = false;
const int dim = Dimension();
double ws[6*16];
double* v;
point.x = 0.0;
point.y = 0.0;
point.z = 0.0;
ds.x = 0.0;
ds.y = 0.0;
ds.z = 0.0;
dt.x = 0.0;
dt.y = 0.0;
dt.z = 0.0;
dss.x = 0.0;
dss.y = 0.0;
dss.z = 0.0;
dst.x = 0.0;
dst.y = 0.0;
dst.z = 0.0;
dtt.x = 0.0;
dtt.y = 0.0;
dtt.z = 0.0;
if ( dim <= 16 ) {
v = ws;
}
else {
v = (double*)onmalloc(6*dim*sizeof(*v));
}
rc = Evaluate( s, t, 2, dim, v, side, hint );
point.x = v[0];
ds.x = v[dim];
dt.x = v[2*dim];
dss.x = v[3*dim];
dst.x = v[4*dim];
dtt.x = v[5*dim];
if ( dim > 1 ) {
point.y = v[1];
ds.y = v[dim+1];
dt.y = v[2*dim+1];
dss.y = v[3*dim+1];
dst.y = v[4*dim+1];
dtt.y = v[5*dim+1];
if ( dim > 2 ) {
point.z = v[2];
ds.z = v[dim+2];
dt.z = v[2*dim+2];
dss.z = v[3*dim+2];
dst.z = v[4*dim+2];
dtt.z = v[5*dim+2];
if ( dim > 16 )
onfree(v);
}
}
return rc;
}
bool
ON_Surface::EvNormal( // returns false if unable to evaluate
double s, double t, // evaluation parameters (s,t)
ON_3dVector& normal, // unit normal
int side, // optional - determines which side to evaluate from
// 0 = default
// 1 from NE quadrant
// 2 from NW quadrant
// 3 from SW quadrant
// 4 from SE quadrant
int* hint // optional - evaluation hint (int[2]) used to speed
// repeated evaluations
) const
{
ON_3dPoint point;
ON_3dVector ds, dt;
return EvNormal( s, t, point, ds, dt, normal, side, hint );
}
bool
ON_Surface::EvNormal( // returns false if unable to evaluate
double s, double t, // evaluation parameters (s,t)
ON_3dPoint& point, // returns value of surface
ON_3dVector& normal, // unit normal
int side, // optional - determines which side to evaluate from
// 0 = default
// 1 from NE quadrant
// 2 from NW quadrant
// 3 from SW quadrant
// 4 from SE quadrant
int* hint // optional - evaluation hint (int[2]) used to speed
// repeated evaluations
) const
{
ON_3dVector ds, dt;
return EvNormal( s, t, point, ds, dt, normal, side, hint );
}
bool
ON_Surface::EvNormal( // returns false if unable to evaluate
double s, double t, // evaluation parameters (s,t)
ON_3dPoint& point, // returns value of surface
ON_3dVector& ds, // first partial derivatives (Ds)
ON_3dVector& dt, // (Dt)
ON_3dVector& normal, // unit normal
int side, // optional - determines which side to evaluate from
// 0 = default
// 1 from NE quadrant
// 2 from NW quadrant
// 3 from SW quadrant
// 4 from SE quadrant
int* hint // optional - evaluation hint (int[2]) used to speed
// repeated evaluations
) const
{
// simple cross product normal - override to support singular surfaces
bool rc = Ev1Der( s, t, point, ds, dt, side, hint );
if ( rc ) {
const double len_ds = ds.Length();
const double len_dt = dt.Length();
// do not reduce the tolerance used here - there is a retry in the code
// below.
if ( len_ds > ON_SQRT_EPSILON*len_dt && len_dt > ON_SQRT_EPSILON*len_ds )
{
ON_3dVector a = ds/len_ds;
ON_3dVector b = dt/len_dt;
normal = ON_CrossProduct( a, b );
rc = normal.Unitize();
}
else
{
// see if we have a singular point
ON_3dVector v[6];
int normal_side = side;
bool bOnSide = false;
ON_Interval sdom = Domain(0);
ON_Interval tdom = Domain(1);
if (s == sdom.Min()) {
normal_side = (normal_side >= 3) ? 4 : 1;
bOnSide = true;
}
else if (s == sdom.Max()) {
normal_side = (normal_side >= 3) ? 3 : 2;
bOnSide = true;
}
if (t == tdom.Min()) {
normal_side = (normal_side == 2 || normal_side == 3) ? 2 : 1;
bOnSide = true;
}
else if (t == tdom.Max()) {
normal_side = (normal_side == 2 || normal_side == 3) ? 3 : 4;
bOnSide = true;
}
if ( !bOnSide )
{
// 2004 November 11 Dale Lear
// Added a retry again with a more generous tolerance
if ( len_ds > ON_EPSILON*len_dt && len_dt > ON_EPSILON*len_ds )
{
ON_3dVector a = ds/len_ds;
ON_3dVector b = dt/len_dt;
normal = ON_CrossProduct( a, b );
rc = normal.Unitize();
}
else
{
rc = false;
}
}
else {
rc = Evaluate( s, t, 2, 3, &v[0].x, normal_side, hint );
if ( rc ) {
rc = ON_EvNormal( normal_side, v[1], v[2], v[3], v[4], v[5], normal);
}
}
}
}
if ( !rc )
{
normal = ON_3dVector::ZeroVector;
}
return rc;
}
//virtual
ON_Curve* ON_Surface::IsoCurve(
int dir, // 0 first parameter varies and second parameter is constant
// e.g., point on IsoCurve(0,c) at t is srf(t,c) - Horizontal
// 1 first parameter is constant and second parameter varies
// e.g., point on IsoCurve(1,c) at t is srf(c,t) - Vertical
double c // value of constant parameter
) const
{
return nullptr;
}
//virtual
bool ON_Surface::Trim(
int dir,
const ON_Interval& domain
)
{
return false;
}
//virtual
bool ON_Surface::Extend(
int dir,
const ON_Interval& domain
)
{
return false;
}
//virtual
bool ON_Surface::Split(
int dir,
double c,
ON_Surface*& west_or_south_side,
ON_Surface*& east_or_north_side
) const
{
return false;
}
// virtual
int ON_Surface::GetNurbForm(
ON_NurbsSurface& nurbs_surface,
double tolerance
) const
{
return 0;
}
// virtual
int ON_Surface::HasNurbForm( ) const
{
return 0;
}
bool ON_Surface::GetSurfaceParameterFromNurbFormParameter(
double nurbs_s, double nurbs_t,
double* surface_s, double* surface_t
) const
{
// NOTE: ON_SumSurface and ON_RevSurface override this virtual function
*surface_s = nurbs_s;
*surface_t = nurbs_t;
return true;
}
bool ON_Surface::GetNurbFormParameterFromSurfaceParameter(
double surface_s, double surface_t,
double* nurbs_s, double* nurbs_t
) const
{
// NOTE: ON_SumSurface and ON_RevSurface override this virtual function
*nurbs_s = surface_s;
*nurbs_t = surface_t;
return true;
}
ON_NurbsSurface* ON_Surface::NurbsSurface(
ON_NurbsSurface* pNurbsSurface,
double tolerance,
const ON_Interval* s_subdomain,
const ON_Interval* t_subdomain
) const
{
ON_NurbsSurface* nurbs_surface = pNurbsSurface;
if ( !nurbs_surface )
nurbs_surface = new ON_NurbsSurface();
int rc = GetNurbForm( *nurbs_surface, tolerance );
if ( !rc )
{
if (!pNurbsSurface)
delete nurbs_surface;
nurbs_surface = nullptr;
}
return nurbs_surface;
}
ON_SurfaceArray::ON_SurfaceArray( int initial_capacity )
: ON_SimpleArray<ON_Surface*>(initial_capacity)
{}
ON_SurfaceArray::~ON_SurfaceArray()
{
Destroy();
}
void ON_SurfaceArray::Destroy()
{
int i = m_capacity;
while ( i-- > 0 ) {
if ( m_a[i] ) {
delete m_a[i];
m_a[i] = nullptr;
}
}
Empty();
}
bool ON_SurfaceArray::Duplicate( ON_SurfaceArray& dst ) const
{
dst.Destroy();
dst.SetCapacity( Capacity() );
const int count = Count();
int i;
ON_Surface* surface;
for ( i = 0; i < count; i++ )
{
surface = 0;
if ( m_a[i] )
{
surface = m_a[i]->Duplicate();
}
dst.Append(surface);
}
return true;
}
bool ON_SurfaceArray::Write( ON_BinaryArchive& file ) const
{
bool rc = file.BeginWrite3dmChunk( TCODE_ANONYMOUS_CHUNK, 0 );
if (rc) rc = file.Write3dmChunkVersion(1,0);
if (rc )
{
int i;
rc = file.WriteInt( Count() );
for ( i = 0; rc && i < Count(); i++ ) {
if ( m_a[i] )
{
rc = file.WriteInt(1);
if ( rc )
rc = file.WriteObject( *m_a[i] ); // polymorphic surfaces
}
else
{
// nullptr surface
rc = file.WriteInt(0);
}
}
if ( !file.EndWrite3dmChunk() )
rc = false;
}
return rc;
}
bool ON_SurfaceArray::Read( ON_BinaryArchive& file )
{
int major_version = 0;
int minor_version = 0;
ON__UINT32 tcode = 0;
ON__INT64 big_value = 0;
int flag;
Destroy();
bool rc = file.BeginRead3dmBigChunk( &tcode, &big_value );
if (rc)
{
rc = ( tcode == TCODE_ANONYMOUS_CHUNK );
if (rc) rc = file.Read3dmChunkVersion(&major_version,&minor_version);
if (rc && major_version == 1) {
ON_Object* p;
int count = 0;
rc = file.ReadInt( &count );
if (rc && count < 0)
rc = false;
if (rc) {
SetCapacity(count);
SetCount(count);
Zero();
int i;
for ( i = 0; rc && i < count && rc; i++ ) {
flag = 0;
rc = file.ReadInt(&flag);
if (rc && flag==1) {
p = 0;
rc = file.ReadObject( &p ); // polymorphic surfaces
m_a[i] = ON_Surface::Cast(p);
if ( !m_a[i] )
delete p;
}
}
}
}
else {
rc = false;
}
if ( !file.EndRead3dmChunk() )
rc = false;
}
return rc;
}
bool ON_Surface::HasBrepForm() const
{
return true;
}
ON_Brep* ON_Surface::BrepForm( ON_Brep* brep ) const
{
ON_Brep* pBrep = nullptr;
if ( brep )
brep->Destroy();
// 26 August 2008 Dale Lear - fixed bug
// When this function was called on an ON_SurfaceProxy
//ON_Surface* pSurface = Duplicate();
ON_Surface* pSurface = DuplicateSurface();
if ( pSurface )
{
if ( brep )
pBrep = brep;
else
pBrep = new ON_Brep();
if ( !pBrep->Create(pSurface) )
{
if ( pSurface )
{
delete pSurface;
pSurface = nullptr;
}
if ( !brep )
delete pBrep;
pBrep = nullptr;
}
}
return pBrep;
}
void ON_Surface::DestroySurfaceTree()
{
DestroyRuntimeCache(true);
}
ON_SurfaceProperties::ON_SurfaceProperties()
{
memset(this,0,sizeof(*this));
}
void ON_SurfaceProperties::Set( const ON_Surface* surface )
{
m_surface = surface;
if ( 0 == m_surface )
{
m_bIsSet = false;
m_bHasSingularity = false;
m_bIsSingular[0] = false;
m_bIsSingular[1] = false;
m_bIsSingular[2] = false;
m_bIsSingular[3] = false;
m_bHasSeam = false;
m_bIsClosed[0] = false;
m_bIsClosed[1] = false;
m_domain[0].Set(0.0,0.0);
m_domain[1].Set(0.0,0.0);
}
else
{
m_bIsSet = true;
m_bIsSingular[0] = m_surface->IsSingular(0)?true:false;
m_bIsSingular[1] = m_surface->IsSingular(1)?true:false;
m_bIsSingular[2] = m_surface->IsSingular(2)?true:false;
m_bIsSingular[3] = m_surface->IsSingular(3)?true:false;
m_bHasSingularity = (m_bIsSingular[0] || m_bIsSingular[1] || m_bIsSingular[2] || m_bIsSingular[3]);
m_bIsClosed[0] = m_surface->IsClosed(0)?true:false;
m_bIsClosed[1] = m_surface->IsClosed(1)?true:false;
m_bHasSeam = (m_bIsClosed[0] || m_bIsClosed[1]);
m_domain[0] = m_surface->Domain(0);
m_domain[1] = m_surface->Domain(1);
}
}