/********************************************************************/ /* Copyright (c) 2017 System fugen G.K. and Yuzi Mizuno */ /* All rights reserved. */ /********************************************************************/ #include "MGCLStdAfx.h" #include "mg/Position_list.h" #include "mg/CCisect_list.h" #include "mg/CSisect_list.h" #include "mg/LBRep.h" #include "mg/RLBRep.h" #include "mg/SurfCurve.h" #include "mg/BSumCurve.h" #include "mg/Transf.h" #include "mg/Plane.h" #include "mg/SPhere.h" #include "mg/Cylinder.h" #include "mg/SBRep.h" #include "mg/RSBRep.h" #include "mg/BSumSurf.h" #include "mg/Tolerance.h" #if defined(_DEBUG) #define new DEBUG_NEW #undef THIS_FILE static char THIS_FILE[] = __FILE__; #endif using namespace std; // MGCylinder.cpp // Implementation of class MGCylinder //MGCylinder is a Cylinder in 3D space. //Cylinder is expressed by an ellipse and a straight line. //Cylinder function f(u,v) = m_ellipse(u) + m_axis(v), //where u and v are two parameter of surface representation. //Here, m_axis is a straight line that passes through the origin. //m_axis'es m_root_point is always set to the origin. //m_ellise is the ellipse when v=0(when v=0, m_axis(v)=the origin); // MGCylinderクラスは3次元空間における円筒面を表すクラスである。 // MGCylinderクラスでは以下のようなパラメータ表現を使用します。 // f(u,v) = m_ellipse(u) + m_axis(v); /////////////Constructor コンストラクタ//////////// //Void constructor 初期化なしで柱面を生成する。 MGCylinder::MGCylinder(void):MGSurface(){;} // Construct a Cylinder by changing this space dimension or // ordering the coordinates. MGCylinder::MGCylinder( int dim, // New space dimension. const MGCylinder& cyl, // Original Cylinder. int start1, // Destination order of new Surface. int start2) // Source order of original Surface. :MGSurface(cyl),m_ellipse(dim,cyl.m_ellipse,start1,start2), m_axis(dim,cyl.m_axis,start1,start2){ update_mark(); m_ortho=m_ellipse.normal().parallel(m_axis.direction()); } // Construct a cylinder of whole circle whose bottom center is bottom //and top center is bottom+axis. MGCylinder::MGCylinder( const MGPosition& bottom, //Location on axis to position the cylinder, //defines zero v. const MGVector axis,//The axis vector for the cylinder. double radius, //The radius of the cylinder. bool outgoing //Indicates if the surface normal is going to //outside of the cylinder(true) or inside(false). ):MGSurface(),m_ellipse(bottom,radius,axis), m_axis(mgORIGIN+axis,mgORIGIN), m_ortho(true){ if(!outgoing) m_ellipse.negate(); } //Cylinder from an ellipse and the axis straight line. //When axis'es start point is not ellipse's center, axis'es start point //will be set to the center. MGCylinder::MGCylinder( const MGEllipse& ellipse, //ellispe of the cylinder const MGStraight& axis //axis of the cylinder. //axis's root point will be neglected, always be set as //the origin. ):MGSurface(),m_ellipse(ellipse),m_axis(axis){ m_axis.m_root_point=mgORIGIN; m_ortho=m_ellipse.normal().parallel(axis.direction()); } //////////Operator overload 演算子の多重定義///////////// //Assignment. //When the leaf object of this and srf2 are not equal, this assignment //does nothing. MGCylinder& MGCylinder::operator=(const MGCylinder& gel2){ if(this==&gel2) return *this; MGSurface::operator=(gel2); m_axis=gel2.m_axis; m_ellipse=gel2.m_ellipse; m_ortho=gel2.m_ortho; return *this; } MGCylinder& MGCylinder::operator=(const MGGel& gel2){ const MGCylinder* gel2_is_this=dynamic_cast<const MGCylinder*>(&gel2); if(gel2_is_this) operator=(*gel2_is_this); return *this; } //Translation of the Cylinder MGCylinder MGCylinder::operator+ (const MGVector& vec) const{ return MGCylinder(m_ellipse+vec,m_axis); } MGCylinder operator+ (const MGVector& v, const MGCylinder& cyl){ return cyl+v; } //Translation of the Cylinder MGCylinder& MGCylinder::operator+= (const MGVector& vec){ m_ellipse+=vec; if(m_box) (*m_box)+=vec; return *this; } //Translation of the Cylinder MGCylinder MGCylinder::operator- (const MGVector& vec) const{ return MGCylinder(m_ellipse-vec,m_axis); } //Translation of the Cylinder MGCylinder& MGCylinder::operator-= (const MGVector& vec){ m_ellipse-=vec; if(m_box) (*m_box)-=vec; return *this; } //柱面のスケーリングを行い,柱面を作成する。 //Scaling of the Cylinder by a double. MGCylinder MGCylinder::operator* (double s) const{ return MGCylinder(m_ellipse*s, m_axis*s); } //Scaling of the Cylinder by a double. MGCylinder operator* (double scale, const MGCylinder& cyl){ return cyl*scale; } //Scaling of the Cylinder by a double. MGCylinder& MGCylinder::operator*= (double s){ m_ellipse*=s; m_axis*=s; update_mark(); return *this; } //Transformation of the Cylinder by a matrix. MGCylinder MGCylinder::operator* (const MGMatrix& mat) const{ return MGCylinder(m_ellipse*mat, m_axis*mat); } //Transformation of the Cylinder by a matrix. MGCylinder& MGCylinder::operator*= (const MGMatrix& mat){ m_ellipse*=mat; m_axis*=mat; m_ortho=m_ellipse.normal().parallel(m_axis.direction()); update_mark(); return *this; } //Transformation of the Cylinder by a MGTransf. MGCylinder MGCylinder::operator* (const MGTransf& tr) const{ return MGCylinder(m_ellipse*tr, m_axis*tr.affine()); } //Transformation of the Cylinder by a MGTransf. MGCylinder& MGCylinder::operator*= (const MGTransf& tr){ m_ellipse*=tr; m_axis*=tr.affine(); m_ortho=m_ellipse.normal().parallel(m_axis.direction()); update_mark(); return *this; } //Comparison between Cylinder and a surface. bool MGCylinder::operator==(const MGCylinder& srf2)const{ if(m_ellipse!=srf2.m_ellipse) return false; if(m_axis!=srf2.m_axis) return false; return true; } bool MGCylinder::operator<(const MGCylinder& gel2)const{ if(m_ellipse==gel2.m_ellipse) return m_axis<gel2.m_axis; return m_ellipse<gel2.m_ellipse; } bool MGCylinder::operator==(const MGGel& gel2)const{ const MGCylinder* gel2_is_this=dynamic_cast<const MGCylinder*>(&gel2); if(gel2_is_this) return operator==(*gel2_is_this); return false; } bool MGCylinder::operator<(const MGGel& gel2)const{ const MGCylinder* gel2_is_this=dynamic_cast<const MGCylinder*>(&gel2); if(gel2_is_this) return operator<(*gel2_is_this); return false; } /////////Debug function デバッグ関数/////////// // Output virtual function. //Output to ostream メンバデータを標準出力に出力する。 std::ostream& MGCylinder::out(std::ostream &outpt) const{ outpt<<"MGCylinder::"<<this<<",m_ortho="<<m_ortho<<std::endl<<" ,m_ellipes="<<m_ellipse <<std::endl<<" ,m_axis="<<m_axis; return outpt; } ////////////Member function メンバ関数/////////// //Compute the axis point of the parameter v. MGPosition MGCylinder::axis_point(double v)const{ return m_ellipse.center()+m_axis.eval(v); } // 入力のパラメータ範囲の曲線部分を囲むボックスを返す。 //Box that includes limitted Cylinder by box. MGBox MGCylinder::box_limitted( const MGBox& uvrange // Parameter Range of the surface. ) const{ MGEllipse el(m_ellipse); el.limit(uvrange[0]); MGStraight axis(m_axis); axis.limit(uvrange[1]); MGBox bx=el.box()+axis.start_point(); bx|=(el+=axis.end_point()).box(); return bx; } //Obtain ceter coordinate of the geometry. MGPosition MGCylinder::center() const{ return m_ellipse.center()+m_axis.center(); } //Changing this object's space dimension. MGCylinder& MGCylinder::change_dimension( int sdim, // new space dimension int start1, // Destination order of new object. int start2) // Source order of this object. { m_ellipse.change_dimension(sdim,start1,start2); m_axis.change_dimension(sdim,start1,start2); update_mark(); return *this; } //Compute the closest point parameter value (u,v) of this surface //from a point. MGPosition MGCylinder::closest(const MGPosition& point) const{ MGPosition uv; perp_point(point,uv); return uv; } //Return minimum box that includes whole of the surface. //Returned is a newed object pointer. MGBox* MGCylinder::compute_box() const{ MGBox* bx=new MGBox(m_ellipse.box()+m_axis.start_point()); *bx|=(m_ellipse+m_axis.end_point()).box(); return bx; } //Construct new surface object by copying to newed area. //User must delete this copied object by "delete". MGCylinder* MGCylinder::clone() const{ return new MGCylinder(*this); } //Construct new surface object by changing //the original object's space dimension. //User must delete this copied object by "delete". MGCylinder* MGCylinder::copy_change_dimension( int sdim, // new space dimension int start1, // Destination order of new line. int start2) // Source order of this line. const{ return new MGCylinder(sdim,*this,start1,start2); } //Change parameter range, be able to change the direction by providing //t1 greater than t2. MGCylinder& MGCylinder::change_range( int is_u, //if true, (t1,t2) are u-value. if not, v. double t1, //Parameter value for the start of original. double t2 //Parameter value for the end of original. ){ if(is_u) m_ellipse.change_range(t1,t2); else m_axis.change_range(t1,t2); return *this; } // 与えられた点との距離を返却する。 //Return the distace between Cylinder and the point. double MGCylinder::distance(const MGPosition& point) const{ return (eval(closest(point))-point).len(); } //Evaluate surface data. MGVector MGCylinder::eval( double u, double v // Parameter value of the surface. , int ndu // Order of derivative along u. , int ndv // Order of derivative along v. )const{ if(ndu==0){ if(ndv==0) return m_ellipse.eval(u)+m_axis.eval(v); else return m_axis.eval(v,ndv); }else{ if(ndv==0) return m_ellipse.eval(u,ndu); else return mgZERO_VEC; } } //Evaluate right continuous surface data. //Evaluate all positional data, 1st and 2nd derivatives. void MGCylinder::eval_all( double u, double v, // Parameter value of the surface. MGPosition& f, // Positional data. MGVector& fu, // df(u,v)/du MGVector& fv, // df/dv MGVector& fuv, // d**2f/(du*dv) MGVector& fuu, // d**2f/(du**2) MGVector& fvv // d**2f/(dv**2) )const{ f=m_ellipse.eval(u)+m_axis.eval(v); fu=m_ellipse.eval(u,1); fv=m_axis.eval(v,1); fuu=m_ellipse.eval(u,2); fuv=fvv=mgZERO_VEC; } // Exchange parameter u and v. //This is not allowed. MGSurface& MGCylinder::exchange_uv(){ assert(false);return *this;} //Modify the original Surface by extrapolating the specified perimeter. //The extrapolation is C2 continuous if the order >=4. //The extrapolation is done so that extrapolating length is "length" //at the position of the parameter value "param" of the perimeter. MGCylinder& MGCylinder::extend( int perimeter, //perimeter number of the Surface. // =0:v=min, =1:u=max, =2:v=max, =3:u=min. double param, // parameter value of above perimeter. double length, //chord length to extend at the parameter param of the perimeter. double dk){ //Coefficient of how curvature should vary at // extrapolation start point. When dk=0, curvature keeps same, i.e. // dK/dS=0. When dk=1, curvature becomes zero at length extrapolated point, // i.e. dK/dS=-K/length at extrapolation start point. // (S=parameter of arc length, K=Curvature at start point) // That is, when dk reaches to 1 from 0, curve changes to flat. assert(0<=perimeter && perimeter<4); bool is_start=true;//starting perimeter MGCurve* to_extend; int ndu=0,ndv=0; if(perimeter==1 || perimeter==3){ // Extrapolate to u-direction to_extend=&m_ellipse; if(perimeter==1) is_start=false;//ending perimeter ndu=1; }else{ // Extrapolate to v-direction to_extend=&m_axis; if(perimeter==2) is_start=false;//ending perimeter ndv=1; } MGPosition uv=perimeter_uv(perimeter,param);//Surface parameter value of param. double vlen=eval(uv,ndu,ndv).len(); double slen=length/vlen; to_extend->extend(slen,is_start); return *this; } //Test if the surface is flat or not within the parameter value rectangle of uvbox. //Function's return value is: // true: if the surface is flat, // false: if the surface is not falt. //When this is not falt, the direction that indicates which direction the surface //should be divided will be output. //***** the flatness is tested only approximately. This is for exclusive use of //planar(). bool MGCylinder::flat( const MGBox& uvbox, double tol, //Tolerance allowed to regard flat //(Allowed distance from a Cylinder). int& direction, // 1: u-direction is more non flat. // 0: v-direction is more non flat. MGPosition& P, //Position of the flat plane will be output. MGUnit_vector& N//Normal of the flat plane will be output. )const{ direction=1; N=normal(uvbox.mid()); const MGInterval& urange=uvbox[0]; const MGInterval& vrange=uvbox[1]; MGEllipse el(mgORIGIN_2D, MGVector(m_ellipse.major_len(),0.), MGVector(0.,m_ellipse.minor_len()), urange); const MGBox& elbx=el.box(); P=elbx.mid(); P+=m_ellipse.center(); P+=m_axis.eval(vrange.mid_point()); double a=elbx[0].length().value(), b=elbx[1].length().value(); if(a<b) a=b; double uspan=el.m_prange[1]-el.m_prange[0]; if(uspan<mgPAI){ MGMatrix mat;mat.set_x_axis(el.end_point()-el.start_point()); el*=mat; const MGBox& elbx2=el.box(); double a2=elbx2[0].length().value(), b2=elbx2[1].length().value(); if(a2<b2) a2=b2; if(a>a2) a=a2; } return a<=tol; } //This is the same as flat except that this does not have the arguments P, N. //***** the flatness is tested only approximately. This is for exclusive use of //tessellation. bool MGCylinder::flat_tess( double u0,double u1, //u range from u0 to u1. double v0,double v1, //v range from v0 to v1. double tol, //Tolerance allowed to regart flat //(Allowed distance from a Cylinder). bool& direction,// 1: u-direction is more non flat. // 0: v-direction is more non flat. double max_edge_len )const{ direction=true; double um=(u0+u1)*0.5; int i; //id of Pn[]. MGVector Pn[3], Nn[3]; Pn[0]=eval(u0, v0); Nn[0]=normal(u0,v0); Pn[1]=eval(um, v0); Nn[1]=normal(um,v0); Pn[2]=eval(u1, v0); Nn[2]=normal(u1,v0); MGVector P=Pn[0]; MGVector VN=Nn[0]; for(i=1; i<3; i++){P+=Pn[i]; VN+=Nn[i];} P/=3.; MGUnit_vector N(VN); double x, d=P%N; double dist[3]; bool is_flat=true; for(i=0; i<3; i++){ x=dist[i]=d-Pn[i]%N; if(x<0.) x=-x; if(x>tol) is_flat=false; } if(is_flat){ MGVector P02(Pn[2],Pn[0]); double lenu=P02%P02; MGVector axisV=m_axis.eval(v0)-m_axis.eval(v1); double lenv=axisV%axisV; direction=lenu>=lenv; if(max_edge_len<=0.) return true; double melen2=max_edge_len; melen2*=melen2; if(direction){ if(lenu>melen2) return false; }else{ if(lenv>melen2) return false; } return true; } return false; } bool MGCylinder::in_range(double u, double v) const{ return m_ellipse.in_range(u) && m_axis.in_range(v); } // Surface と Curve の交点を求める。 //Compute intesection of Cylinder and Curve. MGCSisect_list MGCylinder::isect(const MGCurve& curve) const{ return curve.isect(*this); } // Surface と Curve の交点を求める。 //Compute intesection of Cylinder and Curve. MGCSisect_list MGCylinder::isect(const MGRLBRep& curve) const{ return curve.isect(*this); } // Surface と Curve の交点を求める。 //Compute intesection of Cylinder and Curve. MGCSisect_list MGCylinder::isect(const MGEllipse& curve) const{ return curve.isect(*this); } // Surface と Curve の交点を求める。 //Compute intesection of Cylinder and Curve. MGCSisect_list MGCylinder::isect(const MGLBRep& curve) const{ return curve.isect(*this); } // Surface と Curve の交点を求める。 //Compute intesection of Cylinder and Curve. MGCSisect_list MGCylinder::isect(const MGSurfCurve& curve) const{ return curve.isect(*this); } // Surface と Curve の交点を求める。 //Compute intesection of Cylinder and Curve. MGCSisect_list MGCylinder::isect(const MGBSumCurve& curve) const{ return curve.isect(*this); } // Surface と Surface の交線を求める。 //Compute intesection of Sphere and Surface. MGSSisect_list MGCylinder::isect(const MGSurface& srf2) const{ MGSSisect_list list=srf2.isect(*this); list.replace12(); return list; } MGSSisect_list MGCylinder::isect(const MGPlane& srf2) const{ return intersectPl(srf2); } MGSSisect_list MGCylinder::isect(const MGSphere& srf2) const{ return intersect(srf2); } MGSSisect_list MGCylinder::isect(const MGCylinder& srf2) const{ return intersect(srf2); } MGSSisect_list MGCylinder::isect(const MGSBRep& srf2) const{ return intersect(srf2); } MGSSisect_list MGCylinder::isect(const MGRSBRep& srf2) const{ return intersect(srf2); } MGSSisect_list MGCylinder::isect(const MGBSumSurf& srf2) const{ return intersect(srf2); } //isect_direction() is used by isect_startPt() to define which constant //parameter line should be used to compute intersection, and what //incremental value be used for the parameter. //Function's return value is direction to get next intersection(with dt). //When =1: u=const direction, =0: v=const, =-1: cannot get intersection. int MGCylinder::isect_direction( const MGFSurface& sf2, //Second surface for the intersection. int m1, //id of uvuvS that indicates this surface's parameter //position in uvuvS. (uvuvS(m1), uvuvS(m1+1))=(u,v) of this surface. MGPosition& uvuvS,//start parameter (u,v) pair of this surface and sf2. double& du, //Incremental value of the parameter kind of kdt will be output. double& dv, //Right dt will be output according to the function's output =0,1. double acuRatio //acuracy ratio. )const{ const MGPlane* pl=dynamic_cast<const MGPlane*>(&sf2); if(!pl) return MGFSurface::isect_direction(sf2,m1,uvuvS,du,dv,acuRatio); double du2,dv2; isect_dt(uvuvS[m1],uvuvS[m1+1],du2,dv2,acuRatio); du=du2; dv=dv2; const MGUnit_vector& pln=pl->normal(); MGUnit_vector axisv=m_axis.direction(); double angl=pln%axisv; if(fabs(angl)<.01) return 0; return 1; } #define MGCylinder_isect_div_num1 12. #define MGCylinder_isect_div_num2 5. //isect_dt computes incremental values of u and v direction for the intersection //computation at parameter position (u,v). void MGCylinder::isect_dt( double u, double v, double& du, double& dv, double acuRatio //acuracy ratio. )const{ double alfa=isect_dt_coef(0); du=m_ellipse.param_span()/MGCylinder_isect_div_num1*alfa*acuRatio; dv=m_axis.param_span()/MGCylinder_isect_div_num2*alfa*acuRatio; } //"isect1_incr_pline" is a dedicated function of isect_start_incr, will get // shortest parameter line necessary to compute intersection. MGCurve* MGCylinder::isect_incr_pline( const MGPosition& uv, //last intersection point. int kdt, //Input if u=const v-parameter line or not. // true:u=const, false:v=const. double du, double dv,//Incremental parameter length. double& u, //next u value will be output double& v, //next v value will be output int incr //Incremental valuse of B-coef's id. ) const{ //Compute necessary sub-interval length of parameter line of this surface. if(kdt){ u=m_ellipse.range(uv[0]+du); v=uv[1]; //Compute u=const v-parameter line of this surface in pline. return parameter_curve(kdt,u); }else{ v=m_axis.range(uv[1]+dv); u=uv[0]; //Compute v=const u-parameter line in pline. return parameter_curve(kdt,v); } } //"isect_inner_dt" is a dedicated function of isect_startPt, // comutes adequate incremental parameter value(du,dv) and parameter line kind //kdt(u=const or v=const). void MGCylinder::isect_inner_dt( int n, //num of i.p. obtained so far(not include uvnow). const MGPosition& uvnow,//intersection point obtained last(of this). double& du, double& dv, //incremental length from previous to uvnow is input. //New du or dv will be output according to kdt's return value. int& kdt, //Parameter kind used so far is input, will be output as: //=1:parameter line kind(u=const), =0: v=const, //=-1:should halt computation since incremental value is zero. double acuRatio //Accurate ratio. ) const{ double uerr=param_error_u(); double verr=param_error_v(); double abdu=fabs(du), abdv=fabs(dv); double ratio=acuRatio; if(abdu<uerr){ if(abdv<verr){kdt=-1; return;} else kdt=0; }else if(abdv<verr) kdt=1; else{ MGVector dfdu=eval(uvnow,1,0), dfdv=eval(uvnow,0,1); MGVector dfdt=dfdu*du+dfdv*dv; double fuu=dfdu%dfdu, fuv=dfdu%dfdv; double dffu=fuu*du+fuv*dv; double cos2=dffu*dffu/(dfdt%dfdt)/(dfdu%dfdu); double rcz=MGTolerance::rc_zero(); rcz*=5.; rcz*=rcz; if(cos2<=rcz) kdt=0; else kdt=1; } //Define new dt, kdt. double dt,dtold; if(kdt){ dtold=du; dt=m_ellipse.param_span()/MGCylinder_isect_div_num1; }else{ dtold=dv; dt=m_axis.param_span()/MGCylinder_isect_div_num2; } if(dtold<0.) dt=-dt; dt*=isect_dt_coef(n)*ratio; if(n){ //When this is not the 1st call of isect_inner_dt, //dt must not exceed twice or half of the old dt. double dtr1=dt*dt, dtr2=dtold*dtold; if(dtr1 > 2.*dtr2) dt=dtold*2.; else if(dtr1 < .5*dtr2) dt=dtold*.5; } if(kdt) du=dt; else dv=dt; return; } typedef MGPosition_list::iterator plitr; bool uvcompare(plitr i1,plitr i2){return (*i1)[0]<(*i2)[0];}; //Compute the intersection line of this and the plane pl. MGSSisect_list MGCylinder::intersectPl(const MGPlane& pl)const{ MGSSisect_list lst(this,&pl); MGPosition_list uvuv_list2; double error=MGTolerance::set_wc_zero(MGTolerance::line_zero()*.5); //Save the error. isect_boundary(pl,uvuv_list2); MGTolerance::set_wc_zero(error);//Restore the error. int numi=uvuv_list2.size(); if(numi){ int j; std::vector<plitr> iss(numi); plitr i=uvuv_list2.begin(), ie=uvuv_list2.end(); for(j=0; j<numi; j++,i++) iss[j]=i; std::sort(iss.begin(), iss.end(), uvcompare); MGPosition_list uvuv_list; for(j=0; j<numi; j++,i++) uvuv_list.append(*(iss[j])); //std::cout<<uvuv_list; double uspan=m_ellipse.param_span(); double verr=m_axis.param_error(); while(uvuv_list.entries()>1){ plitr pi=uvuv_list.begin(); MGPosition& uvuv1=*(pi++); MGPosition& uvuv2=*pi; if(MGREqual_base(uvuv1[0], uvuv2[0], uspan)){ MGCurve* iline= new MGStraight(pl.eval(uvuv2[2],uvuv2[3]),pl.eval(uvuv1[2],uvuv1[3])); double t1=iline->param_s(), t2=iline->param_e(); MGCurve* param1= new MGStraight(MGPosition(uvuv2[0],uvuv2[1]), MGPosition(uvuv1[0],uvuv1[1])); param1->change_range(t1,t2); MGCurve* param2= new MGStraight(MGPosition(uvuv2[2],uvuv2[3]), MGPosition(uvuv1[2],uvuv1[3])); param2->change_range(t1,t2); lst.append(iline,param1,param2); uvuv_list.pop_front(); uvuv_list.pop_front(); }else if(fabs(uvuv1[1]-uvuv2[1])<=verr){ MGCurve* iline= new MGEllipse(m_ellipse+m_axis.eval(uvuv1[1])); double t1=iline->param_s(), t2=iline->param_e(); MGCurve* param1= new MGStraight(MGPosition(uvuv2[0],uvuv2[1]), MGPosition(uvuv1[0],uvuv1[1])); param1->change_range(t1,t2); MGPosition elcenter=iline->eval((t1+t2)*.5); MGPosition pluvc; pl.on(elcenter,pluvc); MGCurve* param2= new MGEllipse( MGPosition(uvuv1[2],uvuv1[3]), //Starting point pluvc, MGPosition(uvuv2[2],uvuv2[3]));//Ending point param2->change_range(t1,t2); lst.append(iline,param1,param2); uvuv_list.pop_front(); uvuv_list.pop_front(); }else{ lst.append(isect_with_surf(uvuv_list,pl)); } } } return lst; } //Intersection of Surface and a straight line. MGCSisect_list MGCylinder::isectSl( const MGStraight& sl, const MGBox& uvbox //indicates if this surface is restrictied to the parameter //range of uvbox. If uvbox.is_null(), no restriction. )const{ assert(m_ortho);//currently, m_axis must be orthogonal to m_ellipse. //std::cout<<uvbox<<*this<<endl; const MGVector& axisdir=m_axis.direction(); MGMatrix mat;mat.set_axis(axisdir,2); MGEllipse el2(2,m_ellipse*mat);//std::cout<<el2; MGStraight sl21(2,sl*mat);//std::cout<<sl21<<endl; MGCCisect_list isects=el2.isect(sl21);//std::cout<<isects<<endl; MGCSisect_list list(&sl,this); if(!isects.entries()) return list; const MGVector& sldir=sl.direction(); MGVector slbyaxis=sldir*axisdir; MGMatrix mat2; mat2.set_axis(slbyaxis,2); MGStraight sl22(2,sl*mat2);//std::cout<<sl22; MGCCisect_list::CCiterator i=isects.begin(), ie=isects.end(); for(; i!=ie; i++){ double u=(*i).param1(), t1=(*i).param2(); MGStraight* slcyl=static_cast<MGStraight*>(parameter_curve(true,u)); MGStraight slcyl2(2,(*slcyl)*mat2);//std::cout<<slcyl2; delete slcyl; MGCCisect visect; slcyl2.relation(sl22,visect);//std::cout<<visect; double v=visect.param1(), t2=visect.param2(); double t=(t1+t2)*.5;//To increase the accuracy. if(!uvbox.is_null()){ if(!uvbox[0].includes(u)) continue; if(!uvbox[1].includes(v)) continue; } list.append(sl.eval(t),t,MGPosition(u,v)); } return list; } // 柱面を反転する。ノーマルを逆方向にする。 //Negate the normal of the Cylinder. void MGCylinder::negate( int is_u)// Negate along u-direction if is_u is ture, // else along v-direction. { if(is_u) m_ellipse.negate(); else m_axis.negate(); } //Obtain parameter value if this surface is negated by "negate()". //Negate along u-direction if is_u is ture, // else along v-direction. MGPosition MGCylinder::negate_param(const MGPosition& uv, int is_u)const{ MGPosition uvnew(uv); if(is_u) uvnew(0)=m_ellipse.negate_param(uv[0]); else uvnew(1)=m_axis.negate_param(uv[1]); return uvnew; } //C1連続曲面の一定オフセット関数 //オフセット方向は、ノーマル方向を正とする。トレランスはline_zero()を使用している。 //戻り値は、オフセット面のオートポインターが返却される。 //Surface offset. positive offset value is offset normal direction. //the radius of curvature must be larger than offset value. //line_zero() is used. std::unique_ptr<MGSurface> MGCylinder::offset_c1( double ofs_value, //オフセット量 int& error//エラーコード 0:成功 -1:面におれがある* // -2:曲率半径以上のオフセット不可 -3:面生成コンストラクタエラー )const{ MGVector S=m_ellipse.eval(m_ellipse.param_s(),1); const MGVector& M=m_ellipse.major_axis(); const MGVector& A=m_axis.direction(); if((S*A)%M <0.) ofs_value*=-1.; double r1=m_ellipse.m_r; double r2=r1+ofs_value; if(r2<=0.){ error=-2; return std::unique_ptr<MGSurface>(); } MGEllipse el2(m_ellipse); double scale=r2/r1; el2.m_m*=scale; el2.m_n*=scale; el2.m_r*=scale; el2.update_mark(); std::unique_ptr<MGSurface> surf(new MGCylinder(el2,m_axis)); return surf; } // 指定点が面上にあるか調べる。(面上ならばtrue) //Test if a point is on the Cylinder. If on the Cylinder, return true. bool MGCylinder::on( const MGPosition& point, //A point to test 指定点 MGPosition& puv //Parameter value of the Cylinder will be //returned. )const{ puv=closest(point); return ((point-eval(puv)).len()<=MGTolerance::wc_zero()); } //Test if input (u,v) is parameter value on a perimeter of the surface. //If u or v is on a perimeter, (u,v) will be updated to the perimeter value. bool MGCylinder::on_a_perimeter( double& u, double& v, //Surface parameter (u,v) int& perim_num //if function returns true, //the perimete rnumber is output. )const{ if(!in_range(u,v)) return false; MGSTRAIGHT_TYPE slt=m_axis.straight_type(); if(slt!=MGSTRAIGHT_UNLIMIT){ double vspan=m_axis.param_span(); double v0=m_axis.param_s(); if(MGREqual_base(v,v0,vspan)){ perim_num=0; v=v0; return true; } if(slt!=MGSTRAIGHT_HALF_LIMIT){ double v1=m_axis.param_e(); if(MGREqual_base(v,v1,vspan)){ perim_num=2; v=v1; return true; } } } double uspan=m_ellipse.param_span(); double u0=m_ellipse.param_s(); if(MGREqual_base(u,u0,uspan)){ perim_num=3; u=u0; return true; } double u1=m_ellipse.param_e(); if(MGREqual_base(u,u1,uspan)){ perim_num=1; u=u1; return true; } return false; } //Obtain parameter space error. double MGCylinder::param_error() const{ double uerror=m_ellipse.param_error(); double verror=m_axis.param_error(); return sqrt(uerror*uerror+verror*verror); } // パラメータ範囲を返す。 //Return parameter range of the Cylinder(Infinite box). MGBox MGCylinder::param_range() const{ return MGBox(m_ellipse.param_range(), m_axis.param_range()); } // Compute parameter curve. //Returned is newed area pointer, and must be freed by delete. MGCurve* MGCylinder::parameter_curve( int is_u //Indicates x is u-value if is_u is true. , double x //Parameter value. //The value is u or v according to is_u. )const{ if(is_u){ MGStraight* sl=new MGStraight(m_axis); sl->m_root_point=m_ellipse.eval(x); return sl; }else{ MGEllipse* el=new MGEllipse(m_ellipse); el->m_center+=m_axis.eval(x); return el; } } //Compute part of the surface limitted by the parameter range bx. //bx(0) is the parameter (us,ue) and bx(1) is (vs,ve). //That is u range is from us to ue , and so on. MGCylinder* MGCylinder::part( const MGBox& uvbx, int multiple //Indicates if start and end knot multiplicities //are necessary. =0:unnecessary, !=0:necessary. )const{ MGCylinder* cyl=new MGCylinder(*this); cyl->m_ellipse.limit(uvbx[0]); cyl->m_axis.limit(uvbx[1]); return cyl; } // i must be < perimeter_num(). //When perimeter_num()==0, this function is undefined. MGCurve* MGCylinder::perimeter_curve(int i) const{ if(i==0){ return new MGEllipse(m_ellipse+m_axis.start_point()); }else if(i==2){ return new MGEllipse(m_ellipse+m_axis.end_point()); }else if(i==1){ return new MGStraight(m_axis+m_ellipse.end_point()); }else{ return new MGStraight(m_axis+m_ellipse.start_point()); } } // 与えられた点にもっとも近い面上の垂直のパラメータ値を返却する。 //Return the nearest perpendicular point of the Cylinder from P. // Function's return value is whether point is obtained(1) or not(0) int MGCylinder::perp_point( const MGPosition& P,// 与えられた点 MGPosition& uv, //Parameter value of the Cylinder will be output const MGPosition* uvguess // guess )const{ assert(m_ortho); double v; int obtained; int range=vrange(P,v); MGEllipse el(m_ellipse); if(range==1){ el+=m_axis.end_point(); obtained=0; }else if(range==-1){ el+=m_axis.start_point(); obtained=0; }else{ el+=m_axis.eval(v); obtained=1; } double u=el.closest(P); uv=MGPosition(u,v); return obtained; } //Return all(actually one) foots of perpendicular straight lines from P. MGPosition_list MGCylinder::perps( const MGPosition& P // Point of a space(指定点) ) const{ assert(m_ortho); MGPosition_list list; double v; int vinrange=vrange(P,v); if(vinrange==0){ MGPosition p2d= P*MGTransf(m_ellipse.major_axis(), m_ellipse.minor_axis(),axis_point(v)); double theta[4]; // Intersection points are maixmum 4. int nump=m_ellipse.perp2d(p2d.ref(0), p2d.ref(1), theta); for(int i=0; i<nump; i++){ double u=m_ellipse.radian_to_gp(theta[i]); if(m_ellipse.in_range(u)) list.append(MGPosition(u,v)); } return list; }else return list; } // 入力パラメータをパラメータ範囲でまるめて返却する。 //Round the input uv into parameter range of the Cylinder, //return the same value as input. MGPosition MGCylinder::range(const MGPosition& uv) const{ return MGPosition(m_ellipse.range(uv[0]), m_axis.range(uv[1])); } //Return the space dimension. int MGCylinder::sdim() const{return 3;} //Obtain the v parameter value of the neareast point from the point P to the axis. //Function's return value is if the parameter value v is in the range of this //cylinder: -1: below the range, 0:in the range, 1:above the range. int MGCylinder::vrange(const MGPosition& P, double& v)const{ MGStraight sl(m_axis); sl+=m_ellipse.center(); v=sl.perp_param(P); MGInterval vrng=sl.param_range(); int obtained; if(vrng<v){ v=vrng[1].value(); obtained=1; }else if(vrng>v){ v=vrng[0].value(); obtained=-1; }else{ obtained=0; } return obtained; } //メンバデータを読み込む関数 void MGCylinder::ReadMembers(MGIfstream& buf){ MGSurface::ReadMembers(buf); m_ellipse.ReadMembers(buf); m_axis.ReadMembers(buf); } //メンバデータを書き込む関数 void MGCylinder::WriteMembers(MGOfstream& buf) const{ MGSurface::WriteMembers(buf); m_ellipse.WriteMembers(buf); m_axis.WriteMembers(buf); }