Page MenuHome
Differential D6807 Diff 22906 extern/opennurbs/opennurbs_revsurface.cpp

Changeset View

Standalone View

View Options

extern/opennurbs/opennurbs_revsurface.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_OBJECT_IMPLEMENT( ON_RevSurface, ON_Surface, "A16220D3-163B-11d4-8000-0010830122F0");
void ON_RevSurface::DestroyRuntimeCache( bool bDelete )
{
ON_Surface::DestroyRuntimeCache(bDelete);
if ( 0 != m_curve )
m_curve->DestroyRuntimeCache(bDelete);
// 15 August 2003 Dale Lear
// Added the call to destroy m_bbox.
m_bbox.Destroy();
}
ON_RevSurface* ON_RevSurface::New()
{
return new ON_RevSurface();
}
ON_RevSurface* ON_RevSurface::New( const ON_RevSurface& rev_surface )
{
return new ON_RevSurface(rev_surface);
}
ON_RevSurface::ON_RevSurface() : m_curve(0),
m_axis( ON_3dPoint::Origin, ON_3dVector::ZAxis ),
m_angle( 0.0, 2.0*ON_PI ),
m_t( 0.0, 2.0*ON_PI ),
m_bTransposed(0)
{
ON__SET__THIS__PTR(m_s_ON_RevSurface_ptr);
}
ON_RevSurface::~ON_RevSurface()
{
Destroy();
}
void ON_RevSurface::Destroy()
{
DestroySurfaceTree();
if ( m_curve)
{
delete m_curve;
m_curve = 0;
}
m_axis.Create( ON_3dPoint::Origin, ON_3dVector::ZAxis );
m_angle.Set(0.0,2.0*ON_PI);
m_t = m_angle;
m_bTransposed = false;
m_bbox.Destroy();
}
ON_RevSurface::ON_RevSurface( const ON_RevSurface& src ) : ON_Surface(src)
{
ON__SET__THIS__PTR(m_s_ON_RevSurface_ptr);
m_curve = src.m_curve ? src.m_curve->Duplicate() : nullptr;
m_axis = src.m_axis;
m_angle = src.m_angle;
m_t = src.m_t;
m_bTransposed = src.m_bTransposed;
m_bbox = src.m_bbox;
}
unsigned int ON_RevSurface::SizeOf() const
{
unsigned int sz = ON_Surface::SizeOf();
if ( m_curve )
sz += m_curve->SizeOf();
return sz;
}
ON__UINT32 ON_RevSurface::DataCRC(ON__UINT32 current_remainder) const
{
if ( m_curve )
current_remainder = m_curve->DataCRC(current_remainder);
current_remainder = ON_CRC32(current_remainder,sizeof(m_axis),&m_axis);
current_remainder = ON_CRC32(current_remainder,sizeof(m_angle),&m_angle);
current_remainder = ON_CRC32(current_remainder,sizeof(m_t),&m_t);
current_remainder = ON_CRC32(current_remainder,sizeof(m_bTransposed),&m_bTransposed);
return current_remainder;
}
ON_RevSurface& ON_RevSurface::operator=( const ON_RevSurface& src )
{
if ( this != &src )
{
Destroy();
ON_Surface::operator=(src);
if ( src.m_curve )
m_curve = src.m_curve->Duplicate();
m_axis = src.m_axis;
m_angle = src.m_angle;
m_t = src.m_t;
m_bTransposed = src.m_bTransposed;
m_bbox = src.m_bbox;
}
return *this;
}
bool ON_RevSurface::SetAngleRadians(
double start_angle_radians,
double end_angle_radians
)
{
bool rc = false;
double d = end_angle_radians-start_angle_radians;
if ( d >= 0.0 )
{
if ( d <= ON_ZERO_TOLERANCE || d > 2.0*ON_PI )
{
end_angle_radians = start_angle_radians + 2.0*ON_PI;
}
m_angle.Set( start_angle_radians, end_angle_radians );
rc = true;
DestroySurfaceTree();
}
return rc;
}
bool ON_RevSurface::SetAngleDegrees (
double start_angle_degrees,
double end_angle_degrees
)
{
return SetAngleRadians(
start_angle_degrees*ON_PI/180.0,
end_angle_degrees*ON_PI/180.0
);
}
bool ON_RevSurface::IsValid( ON_TextLog* text_log ) const
{
if ( !m_curve )
{
if ( text_log )
text_log->Print( "ON_RevSurface.m_curve is nullptr.\n");
return false;
}
if ( !m_curve->IsValid(text_log) )
{
if ( text_log )
text_log->Print( "ON_RevSurface.m_curve is not valid.\n");
return false;
}
int dim = m_curve->Dimension();
if ( dim != 3 )
{
if ( text_log )
text_log->Print( "ON_RevSurface.m_curve->Dimension()=%d (should be 3).\n",dim);
return false;
}
if ( !m_axis.IsValid() )
{
if ( text_log )
text_log->Print( "ON_RevSurface.m_axis is not valid.\n");
return false;
}
if ( !m_angle.IsIncreasing() )
{
if ( text_log )
text_log->Print( "ON_RevSurface.m_angle = (%g,%g) (should be an increasing interval)\n",
m_angle[0],m_angle[1]);
return false;
}
double length = m_angle.Length();
if ( length > 2.0*ON_PI + ON_ZERO_TOLERANCE )
{
if ( text_log )
text_log->Print( "ON_RevSurface.m_angle.Length() = %g (should be <= 2*pi radians).\n",length);
return false;
}
if ( m_angle.Length() <= ON_ZERO_TOLERANCE )
{
if ( text_log )
text_log->Print( "ON_RevSurface.m_angle.Length() = %g (should be > ON_ZERO_TOLERANCE).\n",length);
return false;
}
if ( !m_t.IsIncreasing() )
{
if ( text_log )
text_log->Print( "ON_RevSurface.m_t = (%g,%g) (should be an increasing interval)\n",
m_t[0],m_t[1]);
return false;
}
return true;
}
void ON_RevSurface::Dump( ON_TextLog& dump ) const
{
ON_Object::Dump(dump); // print class id
dump.PushIndent();
if ( m_bTransposed )
dump.Print("Paramerization: (curve,angle)\n");
else
dump.Print("Paramerization: (angle,curve)\n");
dump.Print("Axis: ");
dump.Print(m_axis.from);
dump.Print(" to ");
dump.Print(m_axis.to);
dump.Print("\n");
dump.Print("Rotation angle: %g to %g radians.\n",m_angle[0],m_angle[1]);
dump.Print("Angle evaluation parameter interval: [%g,%g].\n",m_t[0],m_t[1]);
if ( m_curve ) {
dump.Print("Revolute: \n");
dump.PushIndent();
m_curve->Dump(dump);
dump.PopIndent();
}
dump.PopIndent();
}
bool ON_RevSurface::Write( ON_BinaryArchive& file ) const
{
bool rc = file.Write3dmChunkVersion(2,0);
if (rc)
{
rc = file.WriteLine( m_axis );
rc = file.WriteInterval( m_angle );
rc = file.WriteInterval( m_t );
rc = file.WriteBoundingBox( m_bbox );
rc = file.WriteInt( m_bTransposed?1:0 );
if ( m_curve )
{
rc = file.WriteChar((char)1);
if (rc) rc = file.WriteObject(*m_curve);
}
else
{
rc = file.WriteChar((char)0);
}
}
return rc;
}
bool ON_RevSurface::Read( ON_BinaryArchive& file )
{
int major_version = 0;
int minor_version = 0;
char bHaveCurve = 0;
bool rc = file.Read3dmChunkVersion(&major_version,&minor_version);
if (rc && major_version == 1)
{
rc = file.ReadLine( m_axis );
rc = file.ReadInterval( m_angle );
rc = file.ReadBoundingBox( m_bbox );
int bTransposedAsInt = m_bTransposed ? 1 : 0;
rc = file.ReadInt( &bTransposedAsInt );
if (rc)
m_bTransposed = bTransposedAsInt ? true : false;
rc = file.ReadChar( &bHaveCurve );
if ( bHaveCurve )
{
ON_Object* obj = 0;
rc = file.ReadObject(&obj);
if ( obj )
{
m_curve = ON_Curve::Cast(obj);
if ( !m_curve )
delete obj;
}
}
m_t[0] = m_angle.Min();
m_t[1] = m_angle.Max();
}
else if (rc && major_version == 2)
{
rc = file.ReadLine( m_axis );
rc = file.ReadInterval( m_angle );
rc = file.ReadInterval( m_t );
rc = file.ReadBoundingBox( m_bbox );
int bTransposedAsInt = m_bTransposed ? 1 : 0;
rc = file.ReadInt( &bTransposedAsInt );
if (rc)
m_bTransposed = bTransposedAsInt ? true : false;
rc = file.ReadChar( &bHaveCurve );
if ( bHaveCurve )
{
ON_Object* obj = 0;
rc = file.ReadObject(&obj);
if ( obj )
{
m_curve = ON_Curve::Cast(obj);
if ( !m_curve )
delete obj;
}
}
}
return rc;
}
int ON_RevSurface::Dimension() const
{
return 3;
}
bool ON_RevSurface::Transform( const ON_Xform& xform )
{
DestroyRuntimeCache();
TransformUserData(xform);
bool rc = (m_curve) ? m_curve->Transform(xform) : false;
ON_3dVector X, Y, Z;
Z = m_axis.Tangent();
X.PerpendicularTo( Z );
X.Unitize();
Y = ON_CrossProduct( Z, X );
if ( !m_axis.Transform(xform) )
rc = false;
ON_3dVector transZ = m_axis.Tangent();
if ( transZ.Length() == 0.0 )
{
// transformation collapsed axis
m_axis.to = m_axis.from + Z;
}
else
{
// see if axis needs to be reversed.
// (Happens with transformations that
// have negative determinant - like mirroring.)
ON_3dVector transX = xform*X;
ON_3dVector transY = xform*Y;
ON_3dVector transXxY = ON_CrossProduct( transX, transY );
double d = transXxY*transZ;
if ( d < 0.0 )
m_axis.to = m_axis.from - m_axis.Direction();
}
m_bbox.Destroy();
m_bbox = BoundingBox();
return rc;
}
bool ON_RevSurface::Evaluate( // returns false if unable to evaluate
double s, // angle
double t, // curve_parameter
// evaluation parameters
int der_count, // number of derivatives (>=0)
int v_stride, // array stride (>=Dimension())
double* v, // array of length stride*(ndir+1)*(ndir+2)/2
int side, // optional - determines which quadrant 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 ds = 1.0;
double x,y,z;
int i, j, k, src_i, dst_i;
ON_3dPoint pt;
if ( m_bTransposed )
{
x = s; s = t; t = x;
if ( side == 2 ) side = 4; else if ( side == 4 ) side = 2;
}
if ( m_t != m_angle )
{
if ( m_t.m_t[1] != m_t.m_t[0] )
{
ds = (m_angle.m_t[1] - m_angle.m_t[0])/(m_t.m_t[1] - m_t.m_t[0]);
x = m_t.NormalizedParameterAt(s);
y = m_angle.ParameterAt(x);
s = y;
}
}
double a = cos(s);
double b = sin(s);
const double ca[4] = {a, -b, -a, b}; // cosine derivatives
const double sa[4] = {b, a, -b, -a}; // sine derivatives
const int curve_dim = m_curve ? m_curve->Dimension() : 0;
if ( curve_dim == 2 || curve_dim == 3 )
{
int curve_side = 0;
switch(side)
{
case 1:
case 2:
curve_side = 1;
break;
case 3:
case 4:
curve_side = -1;
break;
}
rc = m_curve->Evaluate( t, der_count, v_stride, v, curve_side, hint )?true:false;
if ( rc )
{
ON_3dVector zaxis = m_axis.Tangent();
ON_3dVector xaxis, yaxis;
xaxis.PerpendicularTo(zaxis);
xaxis.Unitize();
yaxis = ON_CrossProduct(zaxis,xaxis);
// move curve derivatives to pure t partial spaces in v[]
if ( curve_dim == 2 )
{
for ( i = der_count; i >= 1; i-- )
{
// move curve derivative to proper spots
src_i = v_stride*i;
dst_i = v_stride*((i+1)*(i+2)/2 - 1);
v[dst_i++] = v[src_i++];
v[dst_i++] = 0.0;
v[dst_i] = v[src_i];
}
}
else
{
for ( i = der_count; i >= 1; i-- )
{
// move curve derivative to proper spots
src_i = v_stride*i;
dst_i = v_stride*((i+1)*(i+2)/2 - 1);
v[dst_i++] = v[src_i++];
v[dst_i++] = v[src_i++];
v[dst_i] = v[src_i];
}
}
// convert location coordinates to local frame with origin at m_axis.from
{
pt = ON_3dPoint(v)-m_axis.from;
v[0] = pt*xaxis;
v[1] = pt*yaxis;
v[2] = pt*zaxis;
}
// convert curve derivative coordinates to local frame
for ( i = 1; i <= der_count; i++ )
{
dst_i = v_stride*((i+1)*(i+2)/2 - 1);
pt = ON_3dPoint(v+dst_i); // pt = curve derivative in world coords
v[dst_i++] = pt*xaxis;
v[dst_i++] = pt*yaxis;
v[dst_i] = pt*zaxis;
}
for ( i = der_count; i >= 0; i-- )
{
// i = total order of derivative
double f = 1.0; // f = chain rule scale factor
for ( j = i; j >= 0; j-- )
{
// j = number of partials w.r.t curve parameter
// i-j = number of partials w.r.t angular parameter
dst_i = v_stride*(i*(i+1)/2 + j); //
src_i = v_stride*((j+1)*(j+2)/2 - 1); // curve derivative
k=(i-j)%4;
a = f*ca[k];
b = f*sa[k];
f *= ds;
// calculate derivative in local frame
x = a*v[src_i] - b*v[src_i+1];
y = b*v[src_i] + a*v[src_i+1];
z = (j<i) ? 0.0 : v[src_i+2];
// store answer in world coordinates
pt = x*xaxis + y*yaxis + z*zaxis;
v[dst_i++] = pt.x;
v[dst_i++] = pt.y;
v[dst_i] = pt.z;
}
}
// translate location
v[0] += m_axis.from.x;
v[1] += m_axis.from.y;
v[2] += m_axis.from.z;
if ( m_bTransposed )
{
for ( i = 1; i <= der_count; i++ )
{
for ( j = 0, k = i; j < k; j++, k-- )
{
dst_i = i*(i+1)/2;
src_i = dst_i + k;
dst_i += j;
src_i *= v_stride;
dst_i *= v_stride;
x = v[src_i]; v[src_i++] = v[dst_i]; v[dst_i++] = x;
x = v[src_i]; v[src_i++] = v[dst_i]; v[dst_i++] = x;
x = v[src_i]; v[src_i] = v[dst_i]; v[dst_i] = x;
}
}
}
}
}
return rc;
}
class ON_RevolutionTensor : public ON_TensorProduct
{
public:
ON_3dPoint O;
ON_3dVector X;
ON_3dVector Y;
ON_3dVector Z;
int DimensionA() const;
int DimensionB() const;
int DimensionC() const;
bool Evaluate( double, // a
const double*, // A
double, // b
const double*, // B
double* // C
);
};
int ON_RevolutionTensor::DimensionA() const
{
return 2;
}
int ON_RevolutionTensor::DimensionB() const
{
return 3;
}
int ON_RevolutionTensor::DimensionC() const
{
return 3;
}
bool ON_RevolutionTensor::Evaluate( double a, const double* ArcPoint, double b, const double* ShapePoint, double* SrfPoint )
{
double x, y, z, c, s, rx, ry, A[2], B[3];
if (a != 1.0) {
A[0] = a*ArcPoint[0];
A[1] = a*ArcPoint[1];
ArcPoint = A;
}
if (b != 1.0) {
B[0] = b*ShapePoint[0];
B[1] = b*ShapePoint[1];
B[2] = b*ShapePoint[2];
ShapePoint = B;
}
x = (ShapePoint[0] - O.x)*X.x + (ShapePoint[1] - O.y)*X.y + (ShapePoint[2] - O.z)*X.z;
y = (ShapePoint[0] - O.x)*Y.x + (ShapePoint[1] - O.y)*Y.y + (ShapePoint[2] - O.z)*Y.z;
z = (ShapePoint[0] - O.x)*Z.x + (ShapePoint[1] - O.y)*Z.y + (ShapePoint[2] - O.z)*Z.z;
c = ArcPoint[0];
s = ArcPoint[1];
rx = c*x - s*y;
ry = s*x + c*y;
SrfPoint[0] = O.x + rx*X.x + ry*Y.x + z*Z.x;
SrfPoint[1] = O.y + rx*X.y + ry*Y.y + z*Z.y;
SrfPoint[2] = O.z + rx*X.z + ry*Y.z + z*Z.z;
return true;
}
int ON_RevSurface::GetNurbForm(class ON_NurbsSurface& srf , double tolerance ) const
{
int rc = 0;
if ( 0 != m_curve )
{
ON_NurbsCurve a, c;
ON_Arc arc;
arc.plane.CreateFromNormal( ON_3dPoint::Origin, ON_3dVector::ZAxis );
arc.radius = 1.0;
arc.SetAngleRadians(m_angle[1]-m_angle[0]);
if ( arc.GetNurbForm(a) )
{
if ( m_t.IsIncreasing() )
a.SetDomain( m_t[0], m_t[1] );
rc = m_curve->GetNurbForm(c,tolerance);
if (rc)
{
if ( 2 == c.m_dim )
{
// Increasing the dimension of a 2d curve to 3d fixes
// was added to make the Scale1D operation in
// bug # 103845 work.
ON_WARNING("ON_RevSurface.m_curve is 2-dimensional.");
c.ChangeDimension(3);
}
if ( 3 != c.m_dim )
{
ON_ERROR("ON_RevSurface.m_curve is not valid.");
return 0;
}
if ( m_angle[0] != 0.0 )
{
c.Rotate( m_angle[0], m_axis.Direction(), m_axis.from );
}
ON_RevolutionTensor rho;
rho.O = m_axis.from;
rho.Z = m_axis.Direction();
rho.Z.Unitize();
rho.X.PerpendicularTo(rho.Z);
rho.X.Unitize();
rho.Y = ON_CrossProduct(rho.Z,rho.X);
rho.Y.Unitize();
if ( !srf.TensorProduct( a, c, rho ) )
{
// Testing for false here prevents crashes
// when and was added as part of investigating
// bug # 103845. A change a few lines up
// made it so that particular bug no longer
// fails to create a nurbs surface.
return 0;
}
// make singular points "spot on"
ON_3dPoint C0 = c.PointAtStart();
ON_3dPoint C1 = c.PointAtEnd();
ON_3dPoint A0, A1;
ON_4dPoint CV;
double t0 = ON_UNSET_VALUE;
double t1 = ON_UNSET_VALUE;
if (m_axis.ClosestPointTo(C0,&t0) && ON_IsValid(t0) )
{
A0 = m_axis.PointAt(t0);
if ( C0.DistanceTo(A0) <= ON_ZERO_TOLERANCE )
{
// SouthPole
int j = 0;
for ( int i = 0; i < srf.m_cv_count[0]; i++ )
{
CV.w = srf.Weight(i,j);
CV.x = CV.w*A0.x;
CV.y = CV.w*A0.y;
CV.z = CV.w*A0.z;
srf.SetCV(i,j,CV);
}
}
}
if (m_axis.ClosestPointTo(C1,&t1) && ON_IsValid(t1) )
{
A1 = m_axis.PointAt(t1);
if ( C1.DistanceTo(A1) <= ON_ZERO_TOLERANCE )
{
// NorthPole
int j = srf.m_cv_count[1]-1;
for ( int i = 0; i < srf.m_cv_count[0]; i++ )
{
CV.w = srf.Weight(i,j);
CV.x = CV.w*A1.x;
CV.y = CV.w*A1.y;
CV.z = CV.w*A1.z;
srf.SetCV(i,j,CV);
}
}
}
if ( m_bTransposed )
srf.Transpose();
}
}
}
return (rc > 0) ? 2 : 0;
}
int ON_RevSurface::HasNurbForm() const
{
if (!IsValid())
return 0;
return 2;
}
bool ON_RevSurface::GetSurfaceParameterFromNurbFormParameter(
double nurbs_s, double nurbs_t,
double* surface_s, double* surface_t
) const
{
// NOTE: overrides ON_Surface virtual function
bool rc = (0 != m_curve);
if ( m_bTransposed )
{
double* pTemp = surface_s; surface_s = surface_t; surface_t = pTemp;
double temp = nurbs_s; nurbs_s = nurbs_t; nurbs_t = temp;
}
*surface_s = nurbs_s;
*surface_t = nurbs_t;
ON_Circle circle( ON_xy_plane, 1.0 );
ON_Arc arc( circle, m_angle );
ON_ArcCurve arc_curve(arc, m_t[0], m_t[1]);
if ( !arc_curve.GetCurveParameterFromNurbFormParameter( nurbs_s, surface_s ) )
rc = false;
if ( m_curve )
{
if (!m_curve->GetCurveParameterFromNurbFormParameter( nurbs_t, surface_t ))
rc = false;
}
return rc;
}
bool ON_RevSurface::GetNurbFormParameterFromSurfaceParameter(
double surface_s, double surface_t,
double* nurbs_s, double* nurbs_t
) const
{
// NOTE: overrides ON_Surface virtual function
bool rc = (0 != m_curve);
if ( m_bTransposed )
{
double temp = surface_s; surface_s = surface_t; surface_t = temp;
double* pTemp = nurbs_s; nurbs_s = nurbs_t; nurbs_t = pTemp;
}
*nurbs_s = surface_s;
*nurbs_t = surface_t;
ON_Circle circle( ON_xy_plane, 1.0 );
ON_Arc arc( circle, m_angle );
ON_ArcCurve arc_curve(arc, m_t[0], m_t[1]);
if ( !arc_curve.GetNurbFormParameterFromCurveParameter( surface_s, nurbs_s ) )
rc = false;
if ( m_curve )
{
if (!m_curve->GetNurbFormParameterFromCurveParameter( surface_t, nurbs_t ))
rc = false;
}
return rc;
}
ON_Arc ON_RevSurface::IsoArc(
double curve_parameter
) const
{
for (;;)
{
if (nullptr == m_curve)
break;
// 8 December 2003 Chuck - fix iso extraction bug
// when m_angle[0] != 0.
ON_Circle circle;
ON_3dPoint P = m_curve->PointAt(curve_parameter);
if (false == P.IsValid())
break;
circle.plane.origin = m_axis.ClosestPointTo(P);
circle.plane.zaxis = m_axis.Tangent();
circle.plane.xaxis = P - circle.plane.origin;
circle.radius = circle.plane.xaxis.Length();
if (!circle.plane.xaxis.Unitize())
{
// 8 December 2003 Dale Lear - get valid zero radius
// arc/circle when revolute hits the axis.
// First: try middle of revolute for x-axis
P = m_curve->PointAt(m_curve->Domain().ParameterAt(0.5));
ON_3dPoint Q = m_axis.ClosestPointTo(P);
circle.plane.xaxis = P - Q;
if (!circle.plane.xaxis.Unitize())
{
// Then: just use a vector perp to zaxis
circle.plane.xaxis.PerpendicularTo(circle.plane.zaxis);
}
}
circle.plane.yaxis = ON_CrossProduct(circle.plane.zaxis, circle.plane.xaxis);
circle.plane.yaxis.Unitize();
circle.plane.UpdateEquation();
ON_Arc arc(circle, m_angle);
return arc;
}
ON_Arc arc;
arc.plane = ON_Plane::NanPlane;
arc.radius = ON_DBL_QNAN;
return arc;
}
ON_Curve* ON_RevSurface::IsoCurve( int dir, double c ) const
{
if ( dir < 0 || dir > 1 || !m_curve )
return nullptr;
ON_Curve* crv = 0;
if ( m_bTransposed )
dir = 1-dir;
if ( dir == 0 )
{
// 8 December 2003 Chuck - fix iso extraction bug
// when m_angle[0] != 0.
ON_Circle circle;
ON_3dPoint P = m_curve->PointAt(c);
circle.plane.origin = m_axis.ClosestPointTo(P);
circle.plane.zaxis = m_axis.Tangent();
circle.plane.xaxis = P - circle.plane.origin;
circle.radius = circle.plane.xaxis.Length();
if ( !circle.plane.xaxis.Unitize() )
{
// 8 December 2003 Dale Lear - get valid zero radius
// arc/circle when revolute hits the axis.
// First: try middle of revolute for x-axis
P = m_curve->PointAt(m_curve->Domain().ParameterAt(0.5));
ON_3dPoint Q = m_axis.ClosestPointTo(P);
circle.plane.xaxis = P-Q;
if ( !circle.plane.xaxis.Unitize() )
{
// Then: just use a vector perp to zaxis
circle.plane.xaxis.PerpendicularTo(circle.plane.zaxis);
}
}
circle.plane.yaxis = ON_CrossProduct( circle.plane.zaxis, circle.plane.xaxis );
circle.plane.yaxis.Unitize();
circle.plane.UpdateEquation();
ON_Arc arc( circle, m_angle );
crv = new ON_ArcCurve(arc, m_t[0], m_t[1]);
}
else if ( dir == 1 && m_curve )
{
crv = m_curve->DuplicateCurve();
if ( crv )
{
double a = c;
if ( m_t != m_angle )
{
double t = m_t.NormalizedParameterAt(c);
a = m_angle.ParameterAt(t);
}
if ( a != 0.0 )
{
crv->Rotate( a, m_axis.Direction(), m_axis.from );
}
}
}
return crv;
}
bool ON_RevSurface::Trim( int dir, const ON_Interval& domain )
{
bool rc = false;
if ( dir != 0 && dir != 1 )
return false;
if ( !domain.IsIncreasing() )
return false;
if ( m_bTransposed )
dir = 1-dir;
if ( dir == 0 )
{
ON_Interval dom;
dom.Intersection(domain,m_t);
if ( !dom.IsIncreasing() || !m_t.IsIncreasing() || !m_angle.IsIncreasing() )
return false;
double t0 = m_t.NormalizedParameterAt(dom[0]);
double t1 = m_t.NormalizedParameterAt(dom[1]);
ON_Interval a;
a[0] = m_angle.ParameterAt(t0);
a[1] = m_angle.ParameterAt(t1);
double d = a.Length();
if ( fabs(d) > ON_ZERO_TOLERANCE && fabs(d) <= 2.0*ON_PI+ON_ZERO_TOLERANCE )
{
m_angle = a;
m_t = domain;
rc = true;
}
}
else if ( dir == 1 && m_curve )
{
rc = m_curve->Trim( domain );
}
if ( rc )
{
// update bounding box
ON_BoundingBox bbox0 = m_bbox;
m_bbox.Destroy();
BoundingBox();
if ( m_bbox.IsValid() && bbox0.IsValid() )
m_bbox.Intersection(bbox0);
}
return rc;
}
bool ON_RevSurface::Extend(
int dir,
const ON_Interval& domain
)
{
if ( dir != 0 && dir != 1 )
return false;
if (IsClosed(dir))
return false;
bool do_it = false;
ON_Interval dom = Domain(dir);
if (domain[0] < dom[0])
{
dom[0] = domain[0];
do_it = true;
}
if (domain[1] > dom[1])
{
dom[1] = domain[1];
do_it = true;
}
if (!do_it)
return false;
if ( m_bTransposed )
dir = 1-dir;
bool rc = false;
if ( dir == 0 )
{
double t0 = m_t.NormalizedParameterAt(dom[0]);
double t1 = m_t.NormalizedParameterAt(dom[1]);
ON_Interval a;
a[0] = m_angle.ParameterAt(t0);
a[1] = m_angle.ParameterAt(t1);
if (a.Length() > 2.0*ON_PI+ON_ZERO_TOLERANCE) a[1] = a[0]+2.0*ON_PI;
m_angle = a;
m_t = dom;
rc = true;
}
else if ( dir == 1 && m_curve )
{
rc = m_curve->Extend(dom);
}
if ( rc )
{
DestroySurfaceTree();
// update bounding box
//ON_BoundingBox bbox0 = m_bbox;
m_bbox.Destroy();
BoundingBox();
}
return rc;
}
bool ON_RevSurface::Split(
int dir,
double c,
ON_Surface*& west_or_south_side,
ON_Surface*& east_or_north_side
) const
{
bool rc = false;
ON_RevSurface* srf_ws=ON_RevSurface::Cast(west_or_south_side);
ON_RevSurface* srf_en=ON_RevSurface::Cast(east_or_north_side);
if ( srf_ws && srf_ws == srf_en )
return false;
if ( west_or_south_side && !srf_ws )
return false;
if ( east_or_north_side && !srf_en )
return false;
if ( dir != 0 && dir != 1 )
return false;
if ( m_bTransposed )
dir = 1-dir;
ON_Curve* left_side = 0;
ON_Curve* right_side = 0;
ON_Interval left_angle, right_angle;
ON_Interval left_t, right_t;
left_angle = m_angle;
right_angle = m_angle;
left_t = m_t;
right_t = m_t;
if ( dir == 0 )
{
double t = m_t.NormalizedParameterAt(c);
if ( m_t.Includes(c,true) && 0.0 < t && t < 1.0 )
{
double a = m_angle.ParameterAt( t );
if ( m_angle.Includes(a,true) )
{
rc = true;
left_angle[1] = a;
right_angle[0] = a;
left_t[1] = c;
right_t[0] = c;
if ( srf_ws == this )
left_side = m_curve;
else
left_side = m_curve->Duplicate();
if ( srf_en == this )
right_side = m_curve;
else
right_side = m_curve->Duplicate();
}
}
}
else if ( dir == 1 && m_curve )
{
rc = m_curve->Split( c, left_side, right_side );
if ( rc )
{
if ( this == srf_ws )
{
delete m_curve;
srf_ws->m_curve = left_side;
}
else if ( this == srf_en )
{
delete m_curve;
srf_en->m_curve = right_side;
}
}
}
if ( rc )
{
// save input bounding box
const ON_BoundingBox input_bbox(m_bbox);
if ( !srf_ws )
west_or_south_side = srf_ws = new ON_RevSurface();
else if ( srf_ws != this && srf_ws->m_curve )
{
delete srf_ws->m_curve;
srf_ws->m_curve = 0;
}
if ( !srf_en )
east_or_north_side = srf_en = new ON_RevSurface();
if ( srf_en != this && srf_en->m_curve )
{
delete srf_en->m_curve;
srf_en->m_curve = 0;
}
srf_ws->m_axis = m_axis;
srf_ws->m_angle = left_angle;
srf_ws->m_t = left_t;
srf_ws->m_bTransposed = m_bTransposed;
srf_ws->m_curve = left_side;
srf_ws->m_bbox.Destroy();
srf_en->m_axis = m_axis;
srf_en->m_angle = right_angle;
srf_en->m_t = right_t;
srf_en->m_bTransposed = m_bTransposed;
srf_en->m_curve = right_side;
srf_en->m_bbox.Destroy();
// update bounding boxes
// if input box was valid, intesect calculated boxes with
// input box
srf_ws->BoundingBox(); // calcluates srf_ws->m_bbox
if ( srf_ws->m_bbox.IsValid() && input_bbox.IsValid() )
srf_ws->m_bbox.Intersection(input_bbox);
srf_en->BoundingBox(); // calcluates srf_en->m_bbox
if ( srf_en->m_bbox.IsValid() && input_bbox.IsValid() )
srf_en->m_bbox.Intersection(input_bbox);
}
return rc;
}
bool ON_RevSurface::Transpose()
{
m_bTransposed = m_bTransposed ? false : true;
return true;
}
bool ON_RevSurface::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
{
bool rc = true;
if ( m_curve )
{
double curve_t = m_bTransposed ? s : t;
rc = m_curve->IsContinuous( desired_continuity, curve_t, hint,
point_tolerance, d1_tolerance, d2_tolerance,
cos_angle_tolerance, curvature_tolerance );
}
return rc;
}
bool ON_RevSurface::Reverse( int dir )
{
bool rc = false;
if ( m_bTransposed )
dir = dir ? 0 : 1;
if ( dir == 0 )
{
m_axis.Reverse();
double a0 = m_angle[0];
double a1 = m_angle[1];
m_angle.Set( 2.0*ON_PI - a1, 2.0*ON_PI - a0 );
m_t.Reverse();
rc = true;
}
else if ( dir == 1 && m_curve )
{
rc = m_curve->Reverse();
}
return rc;
}
bool ON_RevSurface::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
bool rc = false;
if ( (m_bTransposed ? 1 : 0) == dir )
{
// angle direction
ON_Circle circle(ON_xy_plane,1.0);
ON_Arc arc( circle, m_angle );
ON_ArcCurve arc_curve(arc, m_t[0], m_t[1]);
rc = arc_curve.GetNextDiscontinuity(
c,
t0,t1,t,
(hint? &hint[dir] : 0),
dtype,cos_angle_tolerance,
curvature_tolerance);
}
else
{
rc = m_curve->GetNextDiscontinuity(
c,
t0,t1,t,
(hint? &hint[dir] : 0),
dtype,cos_angle_tolerance,
curvature_tolerance);
}
return rc;
}
bool ON_RevSurface::IsSingular( int side ) const // 0 = south, 1 = east, 2 = north, 3 = west
{
bool rc = false;
ON_3dPoint P, Q;
//double d, tol;
if ( side < 0 || side > 3 )
return false;
if ( m_bTransposed )
{
switch(side)
{
case 0: side = 3; break;
case 1: side = 2; break;
case 2: side = 1; break;
case 3: side = 0; break;
}
}
if ( 0 == side || 2 == side )
{
P = (0 == side)
? m_curve->PointAtStart()
: m_curve->PointAtEnd();
Q = m_axis.ClosestPointTo(P);
// 26 Feb 2003 Dale Lear
//
// I changed the code to handle cases when the
// "Q" has some coordinates that are small and
// other coordinates that are very large. The
// fabs(Q.*)*ON_SQRT_EPSILON term is required
// in cases where the evaluations in PointAtStart()
// and/or ClosestPointTo() have more numerical
// error than ON_ZERO_TOLERANCE. The numerical
// tolerance term has to be calculated on a
// coordinate-by-coordinate basis. See RR 9683.
//
// old test:
//d = P.DistanceTo(Q);
//if ( d <= ON_ZERO_TOLERANCE )
// rc = true;
// 12 July 2012 Dale Lear
// Use ON_PointsAreCoincident() for all IsSingular queries
//
//d = fabs(P.x - Q.x);
//tol = ON_ZERO_TOLERANCE + fabs(Q.x)*ON_SQRT_EPSILON;
//if ( d <= tol )
//{
// d = fabs(P.y - Q.y);
// tol = ON_ZERO_TOLERANCE + fabs(Q.y)*ON_SQRT_EPSILON;
// if ( d <= tol )
// {
// d = fabs(P.z - Q.z);
// tol = ON_ZERO_TOLERANCE + fabs(Q.z)*ON_SQRT_EPSILON;
// if ( d <= tol )
// rc = true;
// }
//}
rc = ON_PointsAreCoincident(3,0,&P.x,&Q.x);
}
return rc;
}
bool ON_RevSurface::IsPeriodic( int dir ) const // dir 0 = "s", 1 = "t"
{
bool rc = false;
if ( m_bTransposed )
dir = dir ? 0 : 1;
if ( dir == 0 )
{
if (m_angle.Length() >= 2.0*ON_PI-ON_ZERO_TOLERANCE )
rc = true;
}
else if ( dir == 1 && m_curve )
{
rc = m_curve->IsPeriodic();
}
return rc;
}
bool ON_RevSurface::IsPlanar(
ON_Plane* plane,
double tolerance
) const
{
bool rc = false;
if( IsValid()){
ON_Plane AxisNormal( m_curve->PointAtStart(), m_axis.Tangent() );
rc = m_curve->IsInPlane( AxisNormal, tolerance);
if(rc && plane)
*plane = AxisNormal;
}
return rc;
}
bool ON_RevSurface::IsClosed( int dir ) const // dir 0 = "s", 1 = "t"
{
bool rc = false;
if ( m_bTransposed )
dir = dir ? 0 : 1;
if ( dir == 0 )
{
if (m_angle.Length() >= 2.0*ON_PI-ON_ZERO_TOLERANCE )
rc = true;
}
else if ( dir == 1 && m_curve )
{
rc = m_curve->IsClosed();
}
return rc;
}
bool ON_RevSurface::GetParameterTolerance( // returns tminus < tplus: parameters tminus <= s <= tplus
int dir, // 0 gets first parameter, 1 gets second parameter
double t, // t = parameter in domain
double* tminus, // tminus
double* tplus // tplus
) const
{
bool rc = false;
if ( m_bTransposed )
dir = dir ? 0 : 1;
if ( dir == 0 )
{
if ( m_t.IsIncreasing() )
rc = ON_GetParameterTolerance( m_t[0], m_t[1], t, tminus, tplus );
}
else if ( dir == 1 && m_curve )
{
rc = m_curve->GetParameterTolerance( t, tminus, tplus );
}
return rc;
}
bool ON_RevSurface::SetDomain(
int dir, // 0 sets first parameter's domain, 1 gets second parameter's domain
double t0,
double t1
)
{
bool rc = false;
if ( m_bTransposed )
dir = 1-dir;
if ( dir == 0 )
{
if ( t0 < t1 )
{
m_t.Set(t0,t1);
DestroyRuntimeCache();
rc = true;
}
}
else if ( dir == 1 && m_curve )
{
rc = m_curve->SetDomain( t0, t1 ) ? true : false;
DestroyRuntimeCache();
}
return rc;
}
ON_Interval ON_RevSurface::Domain( int dir ) const
{
ON_Interval d;
if ( m_bTransposed )
dir = 1-dir;
if ( dir == 0 )
{
d = m_t;
}
else if ( dir == 1 && m_curve )
{
d = m_curve->Domain();
}
return d;
}
bool ON_RevSurface::GetSurfaceSize(
double* width,
double* height
) const
{
bool rc = false;
if ( m_bTransposed )
{
double* ptr = width;
width = height;
height = ptr;
}
if ( m_curve )
{
rc = true;
ON_Interval cdom = m_curve->Domain();
int i, hint = 0;
int imax = 64;
double d = 1.0/((double)imax);
ON_3dPoint pt0 = ON_3dPoint::UnsetPoint;
ON_3dPoint pt;
double length_estimate = 0.0;
if ( width != nullptr || height != nullptr )
{
double radius_estimate = 0.0;
double r;
for ( i = 0; i <= imax; i++ )
{
if ( m_curve->EvPoint( cdom.ParameterAt(i*d), pt, 0, &hint ) )
{
r = m_axis.DistanceTo(pt);
if ( r > radius_estimate )
radius_estimate = r;
if ( pt0 != ON_3dPoint::UnsetPoint )
length_estimate += pt0.DistanceTo(pt);
pt0 = pt;
}
}
if ( width != nullptr )
*width = m_angle.Length()*radius_estimate;
}
if ( height != nullptr )
{
*height = length_estimate;
}
}
else
{
if ( width )
*width = 0.0;
if ( height )
*height = 0.0;
}
return rc;
}
int ON_RevSurface::SpanCount( int dir ) const
{
int span_count = 0;
if ( m_bTransposed )
dir = 1-dir;
if ( dir==0 && m_t.IsIncreasing() )
{
double a = (0.5 + ON_SQRT_EPSILON)*ON_PI;
double da = fabs(m_angle.Length());
if (da <= a)
span_count = 1;
else if (da <= 2.0*a)
span_count = 2;
else
span_count = 4;
}
else if ( dir == 1 && m_curve )
span_count = m_curve->SpanCount();
return span_count;
}
bool ON_RevSurface::GetSpanVector( int dir, double* s ) const
{
bool rc = false;
if ( m_bTransposed )
dir = 1-dir;
if ( dir==0 && m_t.IsIncreasing() ) {
int span_count = SpanCount(m_bTransposed?1-dir:dir);
if ( span_count > 0 )
{
double d = 1.0/span_count;
s[0] = m_t[0];
for ( int i = 1; i < span_count; i++ )
s[i] = m_t.ParameterAt( i*d );
s[span_count] = m_t[1];
rc = true;
}
}
else if ( dir==1 && m_curve ) {
rc = m_curve->GetSpanVector(s);
}
return rc;
}
int ON_RevSurface::Degree( int dir ) const
{
int span_degree = 0;
if ( m_bTransposed )
dir = 1-dir;
if ( dir == 0 )
span_degree = 2;
else if ( dir == 1 && m_curve )
span_degree = m_curve->Degree();
return span_degree;
}
void ON_RevSurface::ClearBoundingBox()
{
m_bbox.Destroy();
}
bool ON_RevSurface::GetBBox( // returns true if successful
double* boxmin, // boxmin[dim]
double* boxmax, // boxmax[dim]
bool bGrowBox // true means grow box
) const
{
bool rc = m_bbox.IsValid();
if ( !rc )
{
ON_BoundingBox bbox, cbox, abox;
rc = m_curve->GetBoundingBox( cbox );
if (rc)
{
//Dec 16 2010 - Chuck - if the angle range does not include 0, m_curve is not part of the surface.
//bbox = cbox;
ON_3dPointArray corners;
cbox.GetCorners(corners);
ON_3dPoint P;
ON_Arc arc;
arc.plane.zaxis = m_axis.Tangent();
arc.SetAngleRadians(m_angle[1] - m_angle[0]);
int i;
double t;
for ( i = 0; i < corners.Count(); i++ )
{
P = corners[i];
abox.Set(P,false);
while( m_axis.ClosestPointTo(P,&t) ) // not a loop - used for flow control
{
abox.Set(arc.plane.origin,true);
// If we cannot construct a valid arc, then P and the point on the axis
// are added to the bounding box. One case where this happens is when
// P is on the axis of revolution. See bug 84354.
arc.plane.origin = m_axis.PointAt(t);
arc.plane.xaxis = P-arc.plane.origin;
arc.radius = arc.plane.xaxis.Length();
if ( !arc.plane.xaxis.Unitize() )
break;
if ( fabs(arc.plane.xaxis*arc.plane.zaxis) > 0.0001 )
break;
arc.plane.yaxis = ON_CrossProduct(arc.plane.zaxis,arc.plane.xaxis);
if ( !arc.plane.yaxis.Unitize() )
break;
arc.plane.UpdateEquation();
arc.plane.Rotate( m_angle[0], arc.plane.zaxis );
if ( !arc.IsValid() )
break;
abox = arc.BoundingBox();
break;
}
bbox.Union(abox);
}
if ( bbox.IsValid() )
{
ON_RevSurface* ptr = const_cast<ON_RevSurface*>(this);
ptr->m_bbox = bbox;
rc = true;
}
}
}
if ( rc )
{
if ( boxmin )
{
if (bGrowBox){
if (m_bbox.m_min.x < boxmin[0]) boxmin[0] = m_bbox.m_min.x;
if (m_bbox.m_min.y < boxmin[1]) boxmin[1] = m_bbox.m_min.y;
if (m_bbox.m_min.z < boxmin[2]) boxmin[2] = m_bbox.m_min.z;
}
else {
boxmin[0] = m_bbox.m_min.x;
boxmin[1] = m_bbox.m_min.y;
boxmin[2] = m_bbox.m_min.z;
}
}
if ( boxmax )
{
if (bGrowBox){
if (m_bbox.m_max.x > boxmax[0]) boxmax[0] = m_bbox.m_max.x;
if (m_bbox.m_max.y > boxmax[1]) boxmax[1] = m_bbox.m_max.y;
if (m_bbox.m_max.z > boxmax[2]) boxmax[2] = m_bbox.m_max.z;
}
else {
boxmax[0] = m_bbox.m_max.x;
boxmax[1] = m_bbox.m_max.y;
boxmax[2] = m_bbox.m_max.z;
}
}
}
return rc;
}
bool ON_RevSurface::IsSpherical(ON_Sphere* sphere, double tolerance ) const
{
bool rc = false;
if ( m_curve )
{
ON_Plane plane;
ON_Arc arc;
ON_3dPoint P = m_curve->PointAt( m_curve->Domain().Mid() );
plane.origin = m_axis.from;
plane.yaxis = m_axis.Tangent();
plane.zaxis = ON_CrossProduct( P-plane.origin, plane.yaxis );
plane.zaxis.Unitize();
plane.xaxis = ON_CrossProduct( plane.yaxis, plane.zaxis );
plane.UpdateEquation();
if ( plane.IsValid() )
{
if ( m_curve->IsArc( &plane, &arc, tolerance ) )
{
P = m_axis.ClosestPointTo( arc.Center() );
if ( P.DistanceTo(arc.Center()) <= tolerance )
{
rc = true;
if ( sphere )
{
sphere->plane.origin = arc.Center();
sphere->plane.zaxis = m_axis.Tangent();
sphere->plane.yaxis = arc.plane.zaxis;
sphere->plane.xaxis = ON_CrossProduct( sphere->plane.zaxis, sphere->plane.yaxis );
sphere->plane.UpdateEquation();
sphere->radius = arc.radius;
}
}
}
}
}
return rc;
}
bool ON_Surface::IsSphere( ON_Sphere* sphere, double tolerance ) const
{
if ( !ON_IsValid(tolerance) || tolerance <= 0.0 )
tolerance = ON_ZERO_TOLERANCE;
const ON_RevSurface* rs = ON_RevSurface::Cast(this);
if (rs)
{
return rs->IsSpherical(sphere,tolerance) ? true : false;
}
ON_Curve* crv = IsoCurve(0,Domain(1).Mid());
if ( !crv )
return false;
ON_Arc arc0;
int bIsArc0 = crv->IsArc(0,&arc0,tolerance > ON_ZERO_TOLERANCE ? tolerance : 0.0);
delete crv;
crv = 0;
if ( !bIsArc0 )
return false;
crv = IsoCurve(1,Domain(0).Mid());
if ( !crv )
return false;
ON_Arc arc1;
int bIsArc1 = crv->IsArc(0,&arc1,tolerance > ON_ZERO_TOLERANCE ? tolerance : 0.0);
delete crv;
crv = 0;
if ( !bIsArc1 )
return false;
// Determine if one of these arcs is a a portion of a
// great circle.
ON_Sphere sph0;
sph0.plane = arc0.plane;
sph0.radius = arc0.radius;
bool bTestSphere0 = sph0.IsValid();
ON_Sphere sph1;
sph1.plane = arc1.plane;
sph1.radius = arc1.radius;
bool bTestSphere1 = sph1.IsValid();
if ( !bTestSphere0 && !bTestSphere1 )
return false;
double sph0tol = 0.0;
double sph1tol = 0.0;
double tol = 0.5*ON_SQRT_EPSILON*(arc0.radius+arc1.radius);
ON_3dPoint P, S;
double a, d;
for ( a = 0.0; a < 1.0; a += 0.25 )
{
P = arc0.PointAt(a*2.0*ON_PI);
if ( bTestSphere0 )
{
S = sph0.ClosestPointTo(P);
d = S.DistanceTo(P);
if ( d > tol )
{
bTestSphere0 = false;
if ( !bTestSphere1 )
return false;
}
else if ( d > sph0tol )
sph0tol = d;
}
if ( bTestSphere1 )
{
S = sph1.ClosestPointTo(P);
d = S.DistanceTo(P);
if ( d > tol )
{
bTestSphere1 = false;
if ( !bTestSphere0 )
return false;
}
else if ( d > sph1tol )
sph1tol = d;
}
P = arc1.PointAt(a*2.0*ON_PI);
if ( bTestSphere0 )
{
S = sph0.ClosestPointTo(P);
d = S.DistanceTo(P);
if ( d > tol )
{
bTestSphere0 = false;
if ( !bTestSphere1 )
return false;
}
else if ( d > sph0tol )
sph0tol = d;
}
if ( bTestSphere1 )
{
S = sph1.ClosestPointTo(P);
d = S.DistanceTo(P);
if ( d > tol )
{
bTestSphere1 = false;
if ( !bTestSphere0 )
return false;
}
else if ( d > sph1tol )
sph1tol = d;
}
}
// If the arc's are both great circles, then
// both will be true unless we have a bug or
// numerical issues.
if (!bTestSphere0 && !bTestSphere1)
return false;
if ( tol < tolerance )
tol = tolerance;
double u, v;
int sc0 = SpanCount(0);
int sc1 = SpanCount(1);
double* s = (double*)onmalloc( (sc0+sc1+2)*sizeof(s[0]) );
double* t = s + (sc0+1);
GetSpanVector(0,s);
GetSpanVector(1,t);
for ( int i = 0; i < sc0; i++ )
{
for ( int ii = i?1:0; ii <= 4; ii++ )
{
u = 0.25*((4-ii)*s[i] + ii*s[i+1]);
for ( int j = 0; j <= sc1; j++ )
{
for ( int jj = j?1:0; jj <= 4; jj++ )
{
v = 0.25*((4-jj)*t[j] + jj*t[j+1]);
P = PointAt(u,v);
if ( bTestSphere0 )
{
S = sph0.ClosestPointTo(P);
d = S.DistanceTo(P);
if ( d > tol )
{
bTestSphere0 = false;
if ( !bTestSphere1 )
{
onfree(s);
return false;
}
}
else if ( d > sph0tol )
sph0tol = d;
}
if ( bTestSphere1 )
{
S = sph1.ClosestPointTo(P);
d = S.DistanceTo(P);
if ( d > tol )
{
bTestSphere1 = false;
if ( !bTestSphere0 )
{
onfree(s);
return false;
}
}
else if ( d > sph1tol )
sph1tol = d;
}
}
}
}
}
onfree(s);
bool rc = (bTestSphere0 || bTestSphere1);
if ( rc && sphere )
{
if (!bTestSphere0)
*sphere = sph1;
else if (!bTestSphere1)
*sphere = sph0;
else if (sph0tol <= sph1tol)
*sphere = sph0;
else
*sphere = sph1;
}
return rc;
}
bool ON_Surface::IsCylinder( ON_Cylinder* cylinder, double tolerance ) const
{
if ( !ON_IsValid(tolerance) || tolerance <= 0.0 )
tolerance = ON_ZERO_TOLERANCE;
const ON_RevSurface* rs = ON_RevSurface::Cast(this);
bool rc = rs && rs->IsCylindrical(cylinder,tolerance);
if ( !rc && !rs )
{
ON_Curve* crv = IsoCurve(0,Domain(1).Mid());
if ( !crv )
return false;
ON_Arc arc;
ON_Line line;
int bIsLine = 0;
int bIsArc = crv->IsArc(0,&arc,tolerance > ON_ZERO_TOLERANCE ? tolerance : 0.0);
if ( !bIsArc )
{
bIsLine = crv->IsLinear(tolerance > ON_ZERO_TOLERANCE ? tolerance : 0.0);
if ( bIsLine )
{
line.from = crv->PointAtStart();
line.to = crv->PointAtEnd();
}
}
delete crv;
crv = 0;
if ( !bIsArc && !bIsLine )
return false;
crv = IsoCurve(1,Domain(0).Mid());
if ( !crv )
return false;
if ( !bIsArc )
bIsArc = crv->IsArc(0,&arc,tolerance > ON_ZERO_TOLERANCE ? tolerance : 0.0);
else if ( !bIsLine )
{
bIsLine = crv->IsLinear(tolerance > ON_ZERO_TOLERANCE ? tolerance : 0.0);
if ( bIsLine )
{
line.from = crv->PointAtStart();
line.to = crv->PointAtEnd();
}
}
delete crv;
crv = 0;
if ( !bIsArc || !bIsLine )
return false;
double tol = 0.5*ON_SQRT_EPSILON*(arc.radius);
if ( tol < tolerance )
tol = tolerance;
double r = arc.plane.origin.DistanceTo(arc.plane.ClosestPointTo(line.from));
if ( fabs(arc.radius - r) > tol )
return false;
r = arc.plane.origin.DistanceTo(arc.plane.ClosestPointTo(line.to));
if ( fabs(arc.radius - r) > tol )
return false;
ON_3dPoint P;
double u, v;
int sc0 = SpanCount(0);
int sc1 = SpanCount(1);
double* s = (double*)onmalloc( (sc0+sc1+2)*sizeof(s[0]) );
double* t = s + (sc0+1);
GetSpanVector(0,s);
GetSpanVector(1,t);
for ( int i = 0; i < sc0; i++ )
{
for ( int ii = i?1:0; ii <= 4; ii++ )
{
u = 0.25*((4-ii)*s[i] + ii*s[i+1]);
for ( int j = 0; j < sc1; j++ )
{
for ( int jj = j?1:0; jj <= 4; jj++ )
{
v = 0.25*((4-jj)*t[j] + jj*t[j+1]);
P = PointAt(u,v);
r = arc.plane.origin.DistanceTo(arc.plane.ClosestPointTo(P));
if ( fabs(arc.radius - r) > tol )
{
onfree(s);
return false;
}
}
}
}
}
onfree(s);
rc = true;
if ( cylinder )
{
cylinder->Create(arc);
rc = cylinder->IsValid();
}
}
return rc;
}
bool ON_Surface::IsCone( ON_Cone* cone, double tolerance ) const
{
if ( !ON_IsValid(tolerance) || tolerance <= 0.0 )
tolerance = ON_ZERO_TOLERANCE;
const ON_RevSurface* rs = ON_RevSurface::Cast(this);
if (rs)
{
return rs->IsConical(cone,tolerance) ? true : false;
}
ON_Curve* crv = IsoCurve(0,Domain(1).Mid());
if ( !crv )
return false;
ON_Arc arc;
ON_Line line;
int bIsLine = 0;
int bIsArc = crv->IsArc(0,&arc,tolerance > ON_ZERO_TOLERANCE ? tolerance : 0.0);
if ( !bIsArc )
{
bIsLine = crv->IsLinear(tolerance > ON_ZERO_TOLERANCE ? tolerance : 0.0);
if ( bIsLine )
{
line.from = crv->PointAtStart();
line.to = crv->PointAtEnd();
}
}
delete crv;
crv = 0;
if ( !bIsArc && !bIsLine )
return false;
crv = IsoCurve(1,Domain(0).Mid());
if ( !crv )
return false;
if ( !bIsArc )
bIsArc = crv->IsArc(0,&arc,tolerance > ON_ZERO_TOLERANCE ? tolerance : 0.0);
else if ( !bIsLine )
{
bIsLine = crv->IsLinear(tolerance > ON_ZERO_TOLERANCE ? tolerance : 0.0);
if ( bIsLine )
{
line.from = crv->PointAtStart();
line.to = crv->PointAtEnd();
}
}
delete crv;
crv = 0;
if ( !bIsArc || !bIsLine )
return false;
double r0 = arc.plane.origin.DistanceTo(arc.plane.ClosestPointTo(line.from));
double r1 = arc.plane.origin.DistanceTo(arc.plane.ClosestPointTo(line.to));
if ( fabs(r0-r1) <= ON_ZERO_TOLERANCE )
return false;
double h0 = arc.plane.plane_equation.ValueAt(line.from);
double h1 = arc.plane.plane_equation.ValueAt(line.to);
if ( fabs(h0-h1) <= ON_ZERO_TOLERANCE )
return false;
double tol = 0.5*ON_SQRT_EPSILON*(r0+r1);
ON_Cone cn;
cn.height = (h0*r1 - r0*h1)/(r1-r0);
if ( !ON_IsValid(cn.height) || fabs(cn.height) <= ON_ZERO_TOLERANCE )
return false;
cn.plane = arc.plane;
cn.plane.origin = cn.plane.origin + cn.height*cn.plane.zaxis;
cn.plane.UpdateEquation();
if ( r0 >= r1 )
{
cn.radius = r0;
cn.height = h0-cn.height;
}
else
{
cn.radius = r1;
cn.height = h1-cn.height;
}
if ( !cn.IsValid() )
return false;
tol *= fabs(cn.height);
// if (r - h*cn.radius/cn.height > tolerance) return false;
double h = cn.plane.plane_equation.ValueAt(line.from);
if ( fabs(r0*cn.height - h*cn.radius) > tol )
return false;
h = cn.plane.plane_equation.ValueAt(line.to);
if ( fabs(r1*cn.height - h*cn.radius) > tol )
return false;
ON_3dPoint P;
double u, v, r;
int sc0 = SpanCount(0);
int sc1 = SpanCount(1);
double* s = (double*)onmalloc( (sc0+sc1+2)*sizeof(s[0]) );
double* t = s + (sc0+1);
GetSpanVector(0,s);
GetSpanVector(1,t);
tol = 0.5*ON_SQRT_EPSILON*(r0+r1);
if ( tol < tolerance )
tol = tolerance;
tol *= fabs(cn.height);
for ( int i = 0; i < sc0; i++ )
{
for ( int ii = i?1:0; ii <= 4; ii++ )
{
u = 0.25*((4-ii)*s[i] + ii*s[i+1]);
for ( int j = 0; j < sc1; j++ )
{
for ( int jj = j?1:0; jj <= 4; jj++ )
{
v = 0.25*((4-jj)*t[j] + jj*t[j+1]);
P = PointAt(u,v);
h = cn.plane.plane_equation.ValueAt(P);
r = cn.plane.origin.DistanceTo(cn.plane.ClosestPointTo(P));
// if (r - h*cn.radius/cn.height > tolerance) return false;
if ( fabs(r*cn.height - h*cn.radius) > tol )
{
onfree(s);
return false;
}
}
}
}
}
onfree(s);
if ( cone )
{
*cone = cn;
}
return true;
}
bool ON_Surface::IsTorus( ON_Torus* torus, double tolerance ) const
{
if ( !ON_IsValid(tolerance) || tolerance <= 0.0 )
tolerance = ON_ZERO_TOLERANCE;
ON_Curve* crv = IsoCurve(0,Domain(1).Mid());
if ( !crv )
return false;
ON_Arc arc0;
int bIsArc0 = crv->IsArc(0,&arc0,tolerance > ON_ZERO_TOLERANCE ? tolerance : 0.0);
delete crv;
crv = 0;
if ( !bIsArc0 )
return false;
crv = IsoCurve(1,Domain(0).Mid());
if ( !crv )
return false;
ON_Arc arc1;
int bIsArc1 = crv->IsArc(0,&arc1,tolerance > ON_ZERO_TOLERANCE ? tolerance : 0.0);
delete crv;
crv = 0;
if ( !bIsArc1 )
return false;
double tol = 0.5*ON_SQRT_EPSILON*(arc0.radius+arc1.radius);
ON_Torus tr0;
tr0.plane = arc0.plane;
double h = tr0.plane.plane_equation.ValueAt(arc1.plane.origin);
tr0.plane.origin = tr0.plane.origin + h*tr0.plane.zaxis;
tr0.plane.UpdateEquation();
tr0.major_radius = tr0.plane.origin.DistanceTo(arc1.plane.origin);
tr0.minor_radius = arc1.radius;
ON_Torus tr1;
tr1.plane = arc1.plane;
h = tr1.plane.plane_equation.ValueAt(arc0.plane.origin);
tr1.plane.origin = tr1.plane.origin + h*tr1.plane.zaxis;
tr1.plane.UpdateEquation();
tr1.major_radius = tr1.plane.origin.DistanceTo(arc0.plane.origin);
tr1.minor_radius = arc0.radius;
bool bTestTorus0 = tr0.IsValid()?true:false;
bool bTestTorus1 = tr1.IsValid()?true:false;
if ( !bTestTorus0 && !bTestTorus1 )
return false;
double tr0tol = 0.0;
double tr1tol = 0.0;
ON_3dPoint P, T;
double a, d;
for ( a = 0.0; a < 1.0; a += 0.25 )
{
P = arc0.PointAt(a*2.0*ON_PI);
if ( bTestTorus0 )
{
T = tr0.ClosestPointTo(P);
d = T.DistanceTo(P);
if ( d > tol )
{
bTestTorus0 = false;
if ( !bTestTorus1 )
return false;
}
else if ( d > tr0tol )
tr0tol = d;
}
if ( bTestTorus1 )
{
T = tr1.ClosestPointTo(P);
d = T.DistanceTo(P);
if ( d > tol )
{
bTestTorus1 = false;
if ( !bTestTorus0 )
return false;
}
else if ( d > tr1tol )
tr1tol = d;
}
P = arc1.PointAt(a*2.0*ON_PI);
if ( bTestTorus0 )
{
T = tr0.ClosestPointTo(P);
d = T.DistanceTo(P);
if ( d > tol )
{
bTestTorus0 = false;
if ( !bTestTorus1 )
return false;
}
else if ( d > tr0tol )
tr0tol = d;
}
if ( bTestTorus1 )
{
T = tr1.ClosestPointTo(P);
d = T.DistanceTo(P);
if ( d > tol )
{
bTestTorus1 = false;
if ( !bTestTorus0 )
return false;
}
else if ( d > tr1tol )
tr1tol = d;
}
}
// If the arc's planes are perpendicular, then
// both will be true unless we have a bug or
// numerical issues.
if (!bTestTorus0 && !bTestTorus1)
return false;
if ( tol < tolerance )
tol = tolerance;
double u, v;
int sc0 = SpanCount(0);
int sc1 = SpanCount(1);
double* s = (double*)onmalloc( (sc0+sc1+2)*sizeof(s[0]) );
double* t = s + (sc0+1);
GetSpanVector(0,s);
GetSpanVector(1,t);
for ( int i = 0; i < sc0; i++ )
{
for ( int ii = i?1:0; ii <= 4; ii++ )
{
u = 0.25*((4-ii)*s[i] + ii*s[i+1]);
for ( int j = 0; j < sc1; j++ )
{
for ( int jj = j?1:0; jj <= 4; jj++ )
{
v = 0.25*((4-jj)*t[j] + jj*t[j+1]);
P = PointAt(u,v);
if ( bTestTorus0 )
{
T = tr0.ClosestPointTo(P);
d = T.DistanceTo(P);
if ( d > tol )
{
bTestTorus0 = false;
if ( !bTestTorus1 )
{
onfree(s);
return false;
}
}
else if ( d > tr0tol )
tr0tol = d;
}
if ( bTestTorus1 )
{
T = tr1.ClosestPointTo(P);
d = T.DistanceTo(P);
if ( d > tol )
{
bTestTorus1 = false;
if ( !bTestTorus0 )
{
onfree(s);
return false;
}
}
else if ( d > tr1tol )
tr1tol = d;
}
}
}
}
}
onfree(s);
bool rc = (bTestTorus0 || bTestTorus1);
if ( rc && torus )
{
if (!bTestTorus0)
*torus = tr1;
else if (!bTestTorus1)
*torus = tr0;
else if (tr0tol <= tr1tol)
*torus = tr0;
else
*torus = tr1;
}
return rc;
}
static
bool ON__IsCylConeHelper(
const ON_Line& axis,
const ON_Curve* curve,
double tolerance,
ON_Plane& plane,
ON_Line& line,
double r[2],
double& h
)
{
if ( !axis.IsValid() )
return false;
if ( !curve )
return false;
line.from = curve->PointAtStart();
line.to = curve->PointAtEnd();
if ( !line.IsValid() )
return false;
if ( line.Length() <= ON_ZERO_TOLERANCE )
return false;
plane.zaxis = axis.Tangent();
h = plane.zaxis*line.Direction();
if ( !ON_IsValid(h) )
return false;
if ( h < 0.0 )
{
plane.zaxis = -plane.zaxis;
h = -h;
}
if ( h <= ON_ZERO_TOLERANCE )
return false;
double a0 = ON_UNSET_VALUE;
if ( !axis.ClosestPointTo(line.from,&a0) )
return false;
double a1 = ON_UNSET_VALUE;
if ( !axis.ClosestPointTo(line.to,&a1) )
return false;
if ( !ON_IsValid(a0) || !ON_IsValid(a1) )
return false;
ON_3dPoint A0 = axis.PointAt(a0);
ON_3dPoint A1 = axis.PointAt(a1);
plane.origin = A0;
ON_3dVector x0 = line.from - A0;
ON_3dVector x1 = line.to - A1;
r[0] = x0.Length();
r[1] = x1.Length();
if ( x0*x0 < 0.0 && r[0] > ON_ZERO_TOLERANCE && r[1] > ON_ZERO_TOLERANCE )
return false;
plane.xaxis = ( r[0] >= r[1] ) ? x0 : x1;
if (fabs(plane.xaxis.Length()) <= ON_ZERO_TOLERANCE)
return false;
if ( !plane.xaxis.Unitize() )
return false;
plane.yaxis = ON_CrossProduct(plane.zaxis,plane.xaxis);
if ( !plane.yaxis.Unitize() )
return false;
plane.UpdateEquation();
if ( !plane.IsValid() )
return false;
x0 = line.Tangent();
if ( fabs(x0*plane.yaxis) > ON_ZERO_TOLERANCE )
return false;
return (curve->IsLinear( tolerance )?true:false);
}
bool ON_RevSurface::IsCylindrical(
ON_Cylinder* cylinder,
double tolerance
) const
{
ON_Cylinder c;
ON_Line line;
double r[2] = {0.0,0.0};
double h = 0.0;
if ( !ON_IsValid(tolerance) || tolerance <= 0.0 )
tolerance = ON_ZERO_TOLERANCE;
if ( !ON__IsCylConeHelper(m_axis,m_curve,tolerance,c.circle.plane,line,r,h) )
return false;
if ( fabs(r[0]-r[1]) > tolerance )
return false;
if ( fabs(line.Tangent()*c.circle.plane.xaxis) > ON_ZERO_TOLERANCE )
return false;
c.circle.radius = (r[0]==r[1]) ? r[0] : (0.5*(r[0]+r[1]));
c.height[0] = 0.0;
c.height[1] = h;
if ( cylinder )
*cylinder = c;
return c.IsValid();
}
bool ON_RevSurface::IsConical(
ON_Cone* cone,
double tolerance
) const
{
ON_Cone c;
ON_Line line;
double r[2] = {0.0,0.0};
double h = 0.0;
if ( !ON_IsValid(tolerance) || tolerance <= 0.0 )
tolerance = ON_ZERO_TOLERANCE;
if ( !ON__IsCylConeHelper(m_axis,m_curve,tolerance,c.plane,line,r,h) )
return false;
double dr = r[0]-r[1];
if ( fabs(dr) <= ON_ZERO_TOLERANCE )
return false;
if ( 0.0 == r[0] )
{
c.radius = r[1];
c.height = h;
}
else if ( 0.0 == r[1] )
{
c.plane.origin += h*c.plane.zaxis;
c.plane.UpdateEquation();
c.radius = r[0];
c.height = -h;
}
else if ( dr > 0.0 )
{
h *= (r[0]/dr);
c.plane.origin += h*c.plane.zaxis;
c.plane.UpdateEquation();
c.radius = r[0];
c.height = -h;
}
else
{
double d = h*r[0]/dr;
c.plane.origin += d*c.plane.zaxis;
c.plane.UpdateEquation();
c.radius = r[1];
c.height = h-d;
}
if ( cone )
*cone = c;
return c.IsValid();
}