// // Copyright (c) 1993-2022 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 . // //////////////////////////////////////////////////////////////// #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 = (jGetNurbForm(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, intersect calculated boxes with // input box srf_ws->BoundingBox(); // calculates srf_ws->m_bbox if ( srf_ws->m_bbox.IsValid() && input_bbox.IsValid() ) srf_ws->m_bbox.Intersection(input_bbox); srf_en->BoundingBox(); // calculates 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 { // 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; abox.Set(arc.plane.origin,true); arc.radius = arc.plane.xaxis.Length(); if ( !arc.plane.xaxis.Unitize() ) break; double dot = arc.plane.xaxis*arc.plane.zaxis; if ( fabs(dot) > 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(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(); }