/********************************************************************/
/* Copyright (c) 2019 System fugen G.K. and Yuzi Mizuno          */
/* All rights reserved.                                             */
/********************************************************************/
#ifndef _MGFSurface_HH_
#define _MGFSurface_HH_

#include <memory>
#include <vector>
#include <deque>
#include "mg/Default.h"
#include "mg/Unit_vector.h"
#include "mg/Position_list.h"
#include "mg/LBRep.h"
#include "mgGL/VBO.h"

class MGCompositeCurve;
class MGHHisects;
class MGSSisects;
class MGCSisects;
class MGObject;
class MGCurve;
class MGFSurface;
class MGShell;
class MGFace;
class MGDrawParam;
class mgSysGL;
class mgVBO;
class mgTL2Triangles;

/** @addtogroup MGObjectRelated
 *  @{
 */

///MGFSurface is an abstract class to provide the common interfaces to MGFace and MGSurface.
class MG_DLL_DECLR MGFSurface{

friend class MGFace;
friend class MGSurface;
friend class MGCurve;
friend class MGCompositeCurve;
friend class MGTrimmedCurve;
friend class MGStraight;
friend class MGLBRep;
friend class MGSurfCurve;

public:
//////// Constructor ////////

///Null FSurface.
MGFSurface(){;}

///Copy constructor.
MGFSurface(const MGFSurface& fsurf){;};

//////////// Virtual Destructor ////////////
virtual ~MGFSurface(){;};

//////// operator overload////////

///Comparison operator.
bool operator<(const MGFSurface& f2)const;
bool operator>(const MGFSurface& f2)const{return f2<(*this);};

/////////Member function////////

///Return space dimension
virtual int sdimFS() const=0;

///Generate arrow data of the tangent along u and v and the normal
///at the parameter value (u,v) of the FSurface.
///data[0] is the origin of the u-tangent arrow, data[1] is the top of the u-tangent arrow,
///data[2], [3] are two bottoms of u-tangent arrowhead.
///data[0], [4], [5], [6] are the points of v-tangent arrow.
///data[0], [7], [8], [9] are the points of v-tangent arrow.
virtual void arrow(double u,double v, MGPosition data[10])const=0;
virtual void arrow(const MGPosition& uv, MGPosition data[10])const=0;

///Get the box of the object.
const MGBox& get_box() const;

///Return box of the parameter space of the FSurface.
///After trimmed one.
virtual const MGBox box_param2()const=0;

///Get the clone of this MGFSurface.
virtual MGFSurface* clone_fsurface()const=0;

///Get the clone of this as a MGFace.
///If this is MGSurface, it is converted to MGFace.
virtual MGFace* clone_as_face()const=0;

///Compute closest point from a point.
///Returned is the parameter value of the FSurface that is closest to point.
virtual MGPosition closest(
	const MGPosition& point///<the target point.
)const=0;

///Compute closest point from a line to the boundary of the MGFSurface.
///Returned is the parameter value of the FSurface that is closest to sl.
virtual MGPosition closest_on_boundary(
	const MGStraight& sl///<The target line.
)const=0;

///Display direction arrows on the surface.
void display_arrowsFS(
	mgSysGL& sgl
	,int udiv=4
	, int vdiv=4
)const;

///Draw 3D curve in world coordinates.
///The object is converted to curve(s) and is drawn.
void drawWireFS(
	mgVBO& vbo,///<The target graphic object.
	int line_density=1	///<line density to draw a surface in wire mode.
)const;

///Draw 3D curve in world coordinates.
///The object is converted to curve(s) and is drawn.
void drawWireFS_to_highlight(
	mgVBO& vbo,///<The target graphic object.
	int line_density	///<line density to draw a surface in wire mode.
)const;

///Evaluation
///Input parameter value is not checked if it is in_range() or not.
///Even if it is not in_range(), surFSurface evaluation will be executed.
virtual MGVector eval(
	double u,///<FSurface parameter u of value(u,v)
	double v,///<FSurface parameter v of value(u,v)
	int ndu=0,///<Order of derivative along u.
	int ndv=0///<Order of derivative along v.
)const=0;

virtual MGVector eval(
	const MGPosition& uv,///<FSurface parameter value(u,v)
	int ndu=0,///<Order of derivative along u.
	int ndv=0///<Order of derivative along v.
)const=0;

///Evaluate deviations of two faces(this and face2) at npoint
///discrete points.
///(1)Search the common edges which have the distance within tolerance.
///(2)Compute the nearest points from npoint discrete points of this to face2.
///Let uvuvi=uvuvs[i], then
///uvuvi[0], [1] are this face's parameter value(u1,v1), and uvuvi[2], [3] are
///parameter value(u2,v2) of face2 which is the nearest point from the point (u1, v1).
void eval_discrete_deviation(
	const MGFSurface& face2,///< The target 2nd surface.
	std::vector<MGPosition>& uvuvs,///< the nearest points pairs are set.
	int npoint=20,		///<indicates how many discrete points be obtained.
	double tolerance=0.1///<tolerance to get two edge to compute deviation.
)const;

///Obtain all the boundaries(i.e., outer boundary and all the inner boundaries)
std::vector<UniqueCurve> get_all_boundaries(void)const;

///get face pointer if this is MGFace, else null will be returned.
virtual MGFace* get_face_pointer()=0;
virtual const MGFace* get_face_pointer()const=0;

///Get number of inner boundaries as the output of the function.
virtual int get_number_of_boundaries()const=0;

///get surface pointer. Null will never be returned if this is valid MGFSurface.
///That is, if this is MGFace, base surface will be returned.
virtual MGSurface* get_surface_pointer()=0;
virtual const MGSurface* get_surface_pointer()const=0;

///Test if this and 2nd object has common area about their box(),
///taking error into account.
virtual bool has_commonFS(const MGObject& obj2) const=0;

///Test if this FSurface has inner boundary loops or not.
///If this has one, return true.
bool hasInnerBoundaryLoop()const{ return number_of_inner_boundaries()>0;};

///Test if this FSurface has boundary loops or not in the specified box.
///If this has one, return true.
virtual bool hasLoop(const MGBox& uvbox)const{return false;};

///Obtain i-th inner_boundary curves(world coordinates representation)
///of the FSurface. Let the output of inner_boundary(i) be wcurves and
///of inner_boundary_param(i) be pcurves, then wcurves[j] corresponds
///to pcurves[j] one to one. Number of inner_boundary can be obtained
///by the function number_of_inner_boundary().
virtual std::vector<UniqueCurve> inner_boundary(int i)const=0;

///Obtain i-th inner_boundary curves(world coordinates representation)
///of the FSurface. Let the output of inner_boundary(i) be wcurves and
///of inner_boundary_param(i) be pcurves, then wcurves[j] corresponds
///to pcurves[j] one to one. Number of inner_boundary can be obtained
///by the function number_of_inner_boundary().
virtual std::vector<UniqueCurve> inner_boundary_param(int i)const=0;

///Test if parameter value (u,v) is in the range of the FSurface parameter.
virtual bool in_range(double u, double v)const=0;
virtual bool in_range(const MGPosition& uv) const=0;

///Test if (u,v) is inside the face.
///Function's return value is:
///  0:outside the face.
///  1:unknown.
///  2:inside the face, not on a boundary.
///  <0:(u,v) is on an inner boundary, and abs(return code) is the loop id.
///  4:(u,v) is on the outer boundary.
///  >=10: (u,v) is on a perimeter, (10+perimeter number) will be returned.
int in_range_with_on(double u, double v)const{return in_range_with_on(MGPosition(u,v));};
virtual int in_range_with_on(const MGPosition& uv)const=0;

///Intersection.
virtual MGHHisects isectFS(const MGShell& shell2) const;
virtual MGSSisects isectFS(const MGFSurface& fsurf) const;
virtual MGSSisects isectFS(const MGFace& fsurf) const = 0;
virtual MGSSisects isectFS(const MGSurface& fsurf) const = 0;
virtual MGCSisects isectFS(const MGCurve& curv) const = 0;
MGisects isectFS(const MGObject& obj2) const;

///Compute all the intersection points of this face's boundaries(out or inner)
///with face2, and vice versa.
///These intersection points are used to compute surface to surface
///intersection lines.
void intersect12Boundary(
	const MGFSurface& face2,		///<The second surface.
	MGPosition_list& uvuv_list	///<The intersection points will be output,
			///<One member of uvuv_list is (u1,v1,u2,v2), where (u1,v1) is
			///<a parameter of this surface and (u2,v2) is a parameter of 
			///<surf.
)const;

///Compute intersection points of this face's boundary(outer and inners) with
///face2. If intersection points are found and the boundary is a loop,
///the point's edge pointer(of this) will be stored in a member uvuv of uvuvs.
///uvuv[7] is the edge pointer. If the boundary is not a loop(that is, a perimeter of
///Surfaces), uvuv.sdim()==7 and an edge pointer is not returned.
///When uvuv.sdim()==8, the edge pointer of uvuv[7] is accessed through union mgEdgeP.
///uvuvs[i] is i-th intersection points.
virtual int isect_boundary(
	const MGFSurface& face2,///< The target 2nd face.
	MGPosition_list& uvuvs,
	///<id1 and id2 are the ids of uvuv where this face's and f2's parameters
	///<are to be stored in a member of uvuvs.
	///<This face's (u,v) is stored in uvuv(id1) and (id1+1).
	///<f2's (u,v) is stored in uvuv(id2) and (id2+1).
	///<id2=0 if id1=2, and id2=2 if id1=0.
	int id1=0///<id indicator of uvuvs as above.
)const=0;

///Compute intersection points between the boundary of iid-th inner boundary
///of this face and face2 to compute intersections of face with face2.
///Function's return value is the number of ip's obtained before appending
///into uvuv_list, may not be equal to the enlarged size of uvuv_list.
virtual int isect_incurves(
	const MGFSurface& face2,///<2nd surface
	int iid,	///<Inner loop id of this face(from 0)
	MGPosition_list& uvuv_list,	///<intersection points will be appended,
		///<One member in the list is of sdim 8,
		///<(4,5,6) is the direction vector, and (7) is Edge pointer of the point.
	int id1	///<id of uvuv(a member of uvuv_list),
		///<uvuv(id1) for this face parameter uvuv(id2) for face2 parameter,
		///<id2=0 if id1=2, and id2=2 if id1=0.
) const=0;

///Compute intersection points of outer boundary curves of this face 
///with face2 to compute intersections.
///Function's return value is the number of ip's obtained(appended)
///into uvuv_list, may not be equal to the enlarged size of uvuv_list.
virtual int isect_outcurves(
	const MGFSurface& face2,///<2nd surface
	MGPosition_list& uvuv_list,	///<intersection points will be appended,
		///<One member in the list is of sdim 7,
		///<and the last three elements are the ip direction vector.
	int id1			///<id of uvuv(a member of uvuv_list),
		///<uvuv(id1) for this face parameter uvuv(id2) for srf or face2 parameter,
		///<id2=0 if id1=2, and id2=2 if id1=0.
)const=0;

/// "isect_guess" computes one intersection point of surface and a curve,
/// given initail guess parameter values of surface and curve.
///Function's return value is 1(true) if i.p. obtained, and 0(false) if not.
virtual int isect_guess(
	const MGCurve& crv,		///<Curve.
	const MGPosition& uvi,	///<Input initial guess parameter value
						///< of the i.p. of the surface. 
	double ti,			///<Input initial guess parameter value of the line.
	MGPosition& uv,		///< Output parameter value obtained. 
	double& t			///< Output parameter value obtained. 
)const;

/// "isect_guess" computes one intersection point of surface and a curve,
/// given initail guess parameter values of surface and curve.
///Function's return value is 1(true) if i.p. obtained, and 0(false) if not.
virtual int isect_guess(
	const MGStraight& sl,	///<Curve.
	const MGPosition& uvi,	///<Input initial guess parameter value
						///< of the i.p. of the surface. 
	double ti,			///<Input initial guess parameter value of the line.
	MGPosition& uv,		///< Output parameter value obtained. 
	double& t			///< Output parameter value obtained. 
)const{return isect_guess_straight(sl,ti,uvi,t,uv);};

/// "isect_guess" computes one intersection point of surface and a curve,
/// given initail guess parameter values of surface and curve.
///Function's return value is 1(true) if i.p. obtained, and 0(false) if not.
virtual int isect_guess(
	const MGCompositeCurve& crv,	///<Curve
	const MGPosition& uvi,	///<Input initial guess parameter value
						///< of the i.p. of the surface. 
	double ti,			///<Input initial guess parameter value of the line.
	MGPosition& uv,		///< Output parameter value obtained. 
	double& t			///< Output parameter value obtained. 
)const{return isect_guess_composite(crv,uvi,ti,uv,t);};

/// "isect_guess" computes one intersection point of surface and a curve,
/// given initail guess parameter values of surface and curve.
///Function's return value is 1(true) if i.p. obtained, and 0(false) if not.
virtual int isect_guess_composite(
	const MGCompositeCurve& crv,	///<Curve
	const MGPosition& uvi,	///<Input initial guess parameter value
						///< of the i.p. of the surface. 
	double ti,			///<Input initial guess parameter value of the line.
	MGPosition& uv,		///< Output parameter value obtained. 
	double& t			///< Output parameter value obtained. 
)const;

/// "isect_guess_straight" computes one intersection point of surface and
///a straight line, given initail guess parameter values of the surface and 
///the straight line.
///Function's return value is 1(true) if i.p. obtained, and 0(false) if not.
virtual int isect_guess_straight(
	const MGStraight& sl,	///<Straight line.
	double ti,			///<Initial guess parameter value of the straight.
	const MGPosition& uvi,	///<Input initial guess parameter value
						///< of the i.p. of the surface. 
	double& t,			///<Straight parameter obtained.
	MGPosition& uv		///<Surface parameter value obtained(u,v). 
)const;

///Access to i-th element of u knot
virtual double knot_u(int i) const{return 0.0;}

///Access to i-th element of v knot
virtual double knot_v(int i) const{return 0.0;}

///Returns the u knot vector.
virtual const MGKnotVector& knot_vector_u() const=0;
virtual MGKnotVector& knot_vector_u()=0;

///Returns the v knot vector.
virtual const MGKnotVector& knot_vector_v() const=0;
virtual MGKnotVector& knot_vector_v()=0;

///This must be a newed MGFace or MGSurface object.
///If this is a MGFace, returns this pointer.
///If this is a MGSurface, construct a newed MGFace using this newed MGSurface,
///and returns the MGFace*.
virtual MGFace* make_face()=0;

///Negate the FSurface.
virtual void negateFS()=0;

///Compute normal vector(not unit) at uv.
virtual MGVector normal(const MGPosition& uv) const=0;

///Compute normal vector(not unit) at (u,v).
virtual MGVector normal(double u,double v) const=0;

///Get number of inner boundaries as the output of the function.
virtual int number_of_inner_boundaries()const{return 0;};

///Get the object point of this MGFSurface.
virtual const MGObject* object_pointer()const=0;
virtual MGObject* object_pointer()=0;

///Offset.
///distance is plus value if the direction is toward normal vector of the
///FSurface. Minus if against the normal vector.
///エラーコード 0:成功 -1:曲率半径以上のオフセット不可 -3:面生成コンストラクタエラー
virtual int offset_fs(double distance, std::vector<UniqueFSurface>& vecOfsFSurface)const=0;

///Test if a point P is on the FSurface.
///Returned is true if the point P is on the FSurface.
///false(0) if P was not on the FSurface.
virtual bool on(
	const MGPosition& P,///< The target point.
	MGPosition& uv	///<Parameter value of the FSurface is returrned,
					///<Even if P is not on the FSurface, nearest point
					///<parameter value will be returned.
)const=0;

///Test if input (u,v) is parameter value on a perimeter of the base surface.
///If u or v is on a perimeter, they will be updated to the perimeter value.
virtual bool on_a_perimeter(
	double& u,///<Surface parameter u of (u,v)
	double& v,	///<Surface parameter v of (u,v)
	int& perim_num	///<if function returns true,
					///<the perimete number is output.
)const=0;

/// Output virtual function.
virtual std::ostream& outFS(std::ostream& ostrm) const=0;


///Obtain outer_boundary curves(world coordinates representation) of the FSurface.
///Let the output of outer_boundary() be wcurves and of outer_boundary_param()
///be pcurves, then wcurves[i] corresponds to pcurves[i] one to one.
///The output curves can be considered as a continuous counter-clockwise ordered
///boundary of the surface.
virtual std::vector<UniqueCurve> outer_boundary()const=0;

///Obtain boundary curves(parameter space representation) of the FSurface.
///Let the output of boundary() be wcurves and of boundary_parameter()
///be pcurves, then wcurves[i] corresponds to  pcurves[i] one to one.
virtual std::vector<UniqueCurve> outer_boundary_param()const=0;

///Obtain parameter value of the FSurface whose world coordinates are P.
virtual MGPosition param(const MGPosition& P)const=0;

/// Return ending parameter value.
virtual double param_e_u()const=0;
virtual double param_e_v()const=0;

/// Return starting parameter value of the base surface.
virtual double param_s_u()const=0;
virtual double param_s_v()const=0;

///Obtain parameter curves.
///In the case of MGSurface, parameter curve is only one. However, in the case
///of MGFace,  number of parameter curves are more than one.
virtual std::vector<UniqueCurve> parameter_curves(
	int is_u,	///<True(!=0) if x is u-value.(i.e. obtain u=const line)
	double x	///<parameter value. u or v-value accordint to is_u.
)const=0;

///Obtain parameter space error.
double param_error() const;

///Obtain u-direction parameter space error.
double param_error_u() const;

///Obtain v-direction parameter space error.
double param_error_v() const;

///Return parameter value of the middle point of the surface.
///The middle point is defined as the parameter (u,v) where
///u=(param_s_u()+param_e_u())/2, and v likewise.
virtual MGPosition param_mid()const;

/// パラメータ範囲を返す。
///Return parameter range.
virtual MGBox param_range() const=0;

///Return the foot of the perpendicular straight line from P.
///Computation is done from the guess parameter value.
///Function's return value is whether point is obtained(true) or not(false).
virtual bool perp_guess(
	const MGPosition& P,		///<Point
	const MGPosition& uvguess,	///< guess parameter value of the shell
	MGPosition& uv				///< Parameter value will be returned.
)const=0;

///Compute perpendicular points of a curve and the FSurface,
///given guess starting paramter values.
///Function's return value is:
///   perp_guess=true if perpendicular points obtained,
///   perp_guess=false if perpendicular points not obtained,
virtual bool perp_guess(
	const MGCurve& curve,	///<curve.
	const MGPosition& uvguess,	///<Guess parameter value of the FSurface.
	double tguess,			///<Guess parameter value of the curve.
	MGPosition& uv,			///<perpendicular point's parameter values of the shell
	double& t				///<will be output.
)const=0;

///指定点から最も近い、垂線の足とパラメータ値を返す。
///Return the foot of the perpendicular straight line from p that is 
///nearest to point p.
/// Function's return value is whether point is obtained(1) or not(0)
virtual int perp_point(
	const MGPosition& p,///< 指定点(point)
	MGPosition& uv,		///<Parameter value of the surFSurface will be returned.
	const MGPosition* uvguess=0	///< guess parameter value of surFSurface
)const=0;

///Compute perpendicular points on the FSurface from a point P((x,y,z)).

///MGPosition uv in the MGPosition_list is:
///uv(0): u parameter, and uv(1): v parameter of the FSurface.
///Generally number of uv are more than one.
virtual MGPosition_list perps(const MGPosition& P) const=0;

///指定点から最も近い、垂線の足とパラメータ値を返す。

///Return the foot of the perpendicular straight line from p that is 
///nearest to point P.
///Function's return value is whether point is obtained(>0) or not(0)
virtual int perp_one(
	const MGPosition& P, ///< 指定点(point)
	MGPosition& uv 		///<Parameter value of the FSurface will be returned.
)const;

///Obtain the projected curve of a curve onto the FSurface.

///引数vecが与えられないとき、面直投影する。
///投影曲線は面上のパラメータ曲線と3次元曲線としてそれぞれ順番に、vec_crv_uv, vec_crvに格納される。
///uv曲線のトレランスはrc_zero()を、3次元曲線はline_zero()をそれぞれ使用している。
///戻り値：
///		投影曲線の数:		投影曲線が求まった
///		0:			投影曲線が求まらなかった
///		-1:			内部処理エラー
///		-2:			収束処理エラー（収束しなかった）
///The direction of the projection is along the vector vec if the vec is not
///NULL, and normal to the FSurface if the vec is NULL.
///Output of 'project' is two kind of curves:
///one is general world coordinate curves('vec_crv'), and the other is (u,v) curves of
///the parameter space of the FSurface(vec_crv_uv).
///vec_crv_uv.size() is equal to vec_crv.size(). Let the size be n, then
///(vec_crv_uv[i], vec_crv[i]) is one pair for 0<=i<n.
///Function's return value is:
/// >=0: number of curves obtained, <0 : Some error detected.
virtual int project(
	const MGCurve& crv,						///<given curve.
	std::vector<UniqueCurve>& vec_crv_uv,
		///<Projected curve(surface parameter (u,v) representation) will be appended.
	std::vector<UniqueCurve>& vec_crv,
		///<Projected curve(world coordinate(x,y,z) representation) will be appended.
	const MGVector& vec = mgNULL_VEC	///<projection vector,
							///<if vec = NULL then calculate perpendicular project.
)const;

///Obtain the projected curve of a curve onto the FSurface.

///引数vecが与えられないとき、面直投影する。
///投影曲線は3次元曲線としてvec_crvに格納される。
///uv曲線のトレランスはline_zero()を使用している。
///戻り値：
///		投影曲線の数:		投影曲線が求まった
///		0:			投影曲線が求まらなかった
///		-1:			内部処理エラー
///		-2:			収束処理エラー（収束しなかった）
///The direction of the projection is along the vector vec if the vec is not
///NULL, and normal to the FSurface if the vec is NULL.
///Output of 'project' is general world coordinate curves('vec_crv')
///Function's return value is:
/// >=0: number of curves obtained, <0 : Some error detected.
virtual int project(
	const MGCurve& crv,	///<given curve.
	std::vector<UniqueCurve>& vec_crv,
		///<Projected curve(world coordinate(x,y,z) representation) will be appended.
	const MGVector& vec = mgNULL_VEC
		///<if vec = NULL then calculate perpendicular project.
)const;

///Define curve division number when a curve crv be projected onto this MGFSurface.
///The result is used in prj2GetParamRange().
virtual int get_proj_divnum(const MGCurve& crv)const=0;

///Round the input parameter (u,v) of the FSurface to the nearest point of
///the FSurface parameter range.
virtual MGPosition range(const MGPosition& uv)const=0;

///Obtain main parameter lines of the FSurface without boundaries.
///inner_skeleton includes only inner parameter lines without boundaries.
///density indicates how many inner parameter lines are necessary
///for both u and v directions.
std::vector<UniqueCurve> inner_skeleton(int density)const;

///Obtain boundary and main parameter lines of the FSurface.
///skeleton includes boundary() and inner parameter lines.
///density indicates how many inner parameter lines are necessary
///for both u and v directions.
virtual std::vector<UniqueCurve> skeleton(int density=1)const;

///Obtain all the parameter curves at knots of u and v knot vector.
virtual std::vector<UniqueCurve> skeleton_at_knots()const;

///split this fsurface at the parameter param.
virtual void split(
	double param,///<parameter value of this fsurface. if is_u is true, param is u-value,
				///<else v-value.
	bool is_u,	///<indicates if param is u or v of the surface parameter (u,v).
	std::vector<UniqueFSurface>& surfaces///<splitted surfaces will be output.
)const=0;


///split this fsurface with splitters. splitters are 3D (x,y,z) curves that may not
///lie on the surface.
void split(
	const std::vector<const MGCurve*>& splitters,	//splitter world curves.
	const MGVector&  dir,	//splitter projection direction.
							//If dir.is_null(), normal projection will be performed.
	std::vector<std::unique_ptr<MGFace>>& faces//Result splitted face(s) will be appended.
			//If no splitting was performed, no faces will be appended.
)const;

///split this fsurface with splitters. splitters are 2D (u,v) surfaces's parameter curves.
void split(
	const std::vector<UniqueCurve>& splitters,//splitter (u,v) curves.
	std::vector<std::unique_ptr<MGFace>>& faces//Result splitted face(s) will be appended.
			//If no splitting was performed, no faces will be appended.
)const;

///Extract a sub surface with trimmers. trimmers are 3D curves and will be projected
///onto this surface toword the direction dir. If dir is null vector, surface normal
///projection will be performed. Extraction is so performed that the smallest region
///enclosed by trimmers that includes the surface point uv is extracted. 
void extract(
	const std::vector<const MGCurve*>& trimmers,	///<Trimmer curves
	const MGVector&  dir,	///<trimmers projection direction.
	const MGPosition& uv,	///<surface parameter (u,v) that indicates the region to extract.
							///<The smallest region that inclued uv will be extracted.
	std::unique_ptr<MGFace>& eface///<Result extracted face will be output.
)const;
void extract(
	const std::vector<UniqueCurve>& trimmers,	///<Trimmer curves
	const MGVector&  dir,	///<trimmers projection direction.
	const MGPosition& uv,	///<surface parameter (u,v) that indicates the region to extract.
							///<The smallest region that inclued uv will be extracted.
	std::unique_ptr<MGFace>& eface///<Result extracted face will be output.
)const;

///Shade the object in world coordinates.
void shadeFS(
	mgVBO& vbo,
	const MGDrawParam& para,
	MGCL::DRAW_TARGET target= MGCL::SHADING
)const;

///Triangulate this object(MGShell, MGFace, or MGSurface is the target).
virtual void triangulate(
	const MGDrawParam& para,
	MGCL::TL_DATA_KIND dkind,
	std::vector<mgTL2Triangles>& trisVec
)const=0;

///Trim this fsurface with trimmers.

///trimmers are 3D curves and will be projected
///onto this surface tword the direction dir. If dir is null vector, surface normal
///prjection will be performed. Trimming is so performed that the smallest region enclosed
///by trimmers that includes the surface point uv will be removed. 
void trim(
	const std::vector<const MGCurve*>& trimmers,	///<Trimmer curves
	const MGVector&  dir,	///<trimmers projection direction.
	const MGPosition& uv,	///<surface parameter (u,v) that indicates the region to remove,
							///<The smallest region that inclued uv will be removed.
	std::vector<UniqueFace>& faces///<Result trimmed face(s) will be appended,
			///<If no trimming was performed, no faces will be appended.
)const;
void trim(
	const std::vector<UniqueCurve>& trimmers,	///<Trimmer curves
	const MGVector&  dir,	///<trimmers projection direction.
	const MGPosition& uv,	///<surface parameter (u,v) that indicates the region to remove,
							///<The smallest region that inclued uv will be removed.
	std::vector<UniqueFace>& faces///<Result trimmed face(s) will be appended,
			///<If no trimming was performed, no faces will be appended.
)const;

///Compute unit normal vector at uv.
virtual MGUnit_vector unit_normal(const MGPosition& uv) const=0;

///Compute unit normal vector at (u,v).
virtual MGUnit_vector unit_normal(double u,double v) const=0;

protected:

///Obtain coefficient's space dimension.
///This function is used in isect_start etc.
virtual int coef_sdim() const=0;

///isect_area_length() returns initial area length for the intersection
///line.
virtual int isect_area_length() const=0;

///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.
virtual int 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=1.	///<acuracy ratio.
)const;

///isect_direction_with_direction() is used by isect_start() to define which constant
///parameter line should be use to compute intersection, and what
///incremental value be used for the parameter.
///Function's return value isect_direction() is 1 for u=const parameter
///line, and 0 for v=const parameter line.
int isect_direction_with_direction(
	double u,	///<start parameter u of (u,v) of this surface.
	double v,	///<start parameter v of (u,v) of this surface.
	const MGVector& tangent,///<To indicate which direction isect line
							///<should march toward.
	double& du,				///<Incremental value sign of the parameter kind of
	double& dv		///<isect_direction_with_direction.
)const;

///isect_dt computes incremental values du and dv for the intersection
///computation at parameter position (u,v).
void isect_dt(
	double u,///<Surface parameter u of (u,v).
	double v,///< v of (u,v).
	double& du,///< u incremental value is output.
	double& dv,///< v incremental value is output.
	double acuRatio=1.	///<acuracy ratio.
) const;

///"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).
virtual void 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, ///<Input u incremental value,
	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=1.	///<Accurate ratio.
)const;

///isect_dt_coef provides coef of how fine parameter increment should be,
///given num of intersection points computed so far.
double isect_dt_coef(int n) const;

///isect_div_id_max is maximum id of array of sect_div defined in
///isect_dt_coef. That is, isect_div_id_max+1 is the length of the array
///sect_div.
int isect_div_id_max()const;

///"isect1_incr_pline" is a dedicated function of isect_start_incr, will get
/// shortest parameter line necessary to compute intersection.
virtual MGCurve* 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,	///<Incremental u parameter length.
	double dv,	///<Incremental v parameter length.
	double& u,				///<next u value will be output
	double& v,				///<next v value will be output
	int incr=0			///<Incremental valuse of B-coef's id.
)const=0;

///isect_start compute one intersection line of two surfaces, this and sf2,
/// given starting intersetion point uvuv((u1,v1) of this and (u2,v2) of sf2)
/// and direction(stan). isect_start halts the computation when intersection
/// reached to a boundary of this or sf2, or reached to one of the points
/// in uvuv_list.
///The function's return value is:
/// =0: Intersection was not obtained.
/// !=0: Intersection was obtained as follows:
///    =1: End point is a point on a perimeter of one of the surfaces.
///    =3: End point is one of boundary points in uvuv_list.
///    =4: End point is the starting point.
///    =7: isect_start halted the computation since intersection was lost
///     during the computation.
int isect_start(
	const MGPosition& uvuv_startIn, ///<Starting point of the intersection line.
	MGPosition_list& uvuv_list,	///<isect_start will halt when ip reached one of 
		///<the point in uvuv_list. isect_start does not change uvuv_list(actually
		///<uvuv_list is const.) uvuv's space dimension is at least 4,
		///<and the first 2 is (u,v) of this and the next 2 is (u,v) of sf2. 
	const MGFSurface& sf2,	///<2nd surface.
	MGSSisect& ssi,			///<Surface-surface intersection line will be output.
	MGPosition_list::iterator& uvuv_id,
			///<When the end point of ip was one of the points of uvuv_list,
			///<uvuv_list's iterator of the point will be returned, that is, 
			///<when the function's return value was 3 or 5.
			///<When was not a point of uvuv_list, end() of uvuv_list will be
			///<returned.
	int& m1	///<id that indicates which surface was used as the main surface.
			///<m1=0: this surface, m1=2: sf2.
)const;

/// "isect_start_boundary" is a dedicated function of isect_start.
///"isect_start_boundary" computes one intersection point of two surfaces,
/// this surface's parameter line at uv+dt(according to kdt) and sb2,
/// given previous intersetion point(uv) and incremental value dt.
/// "isect_start_boundary" is used only when no intersection found at dt.
///Function's return value is 0: ip not found.
///                           2: ip found as intersection line end point.
int isect_start_boundary(
	const MGFSurface& sf2,	///<2nd surface b-rep.
	const MGPosition& uvuv_pre,	///<Starting parameter values of ip. 
	int kdt,		///<kdt=true: u=const parameter line,
					///<    else: v=const parameter line of this surface.
	double du,///<Incremental u parameter length.
	double dv,///<Incremental v parameter length.
	int lid1,		///<id of parameter of this surface in uvuv_pre or uvuv_now.
	MGPosition& uvuv_now///<New parameter values of ip will be output.
)const;

///Compute the maximum difference between the intersection line and the two surfaces,
///this and sr2. The difference evaluation is done at the data point tau.
double isect_start_dif(
	const MGNDDArray& tau,	///<data points
	const MGLBRep& line,	///<the intersection line of this and sf2.
	const MGFSurface& sf2	///<second surface.
)const;

///isect_start_incr compute one intersection point of two surfaces,
/// this surface's parameter line at uv1+dt(according to kdt) and sf2,
/// given previous intersetion point(uv1,uv2) and incremental value dt.
///Here uv1 is uvuv_pre(lid1, lid1+1), and uv2 is other two values of uvuvpre.
/// isect_start_incr is a dedicated function of isect_start
/// and isect_start_boundary.
///Function's return value is true: if ip found,
///                           false: if ip not found.
int isect_start_incr(
	const MGFSurface& sf2,	///<2nd surface b-rep.
	const MGPosition& uvuv_pre,	///<Starting parameter values of ip. 
	int kdt,		///<kdt=true: u=const parameter line,
					///<    else: v=const parameter line of this surface.
	double du,///<Incremental u parameter length.
	double dv,///<Incremental v parameter length.
	int lid1,	///<id of parameter of this surface in uvuv_pre or uvuv_now.
	MGPosition& uvuv_now///<New parameter values of ip will be output.
)const;

///isect_startPt compute an array of parameter value pairs of this surf and sf2
///for one intersection line of the two surfaces,
/// given starting intersetion point uvuv((u1,v1) of this and (u2,v2) of sf2)
/// and direction(stan).
///isect_startPt is a dedicated function for isect_start.
///isect_startPt halts the computation when intersection
///reached to a boundary of this or sf2, or reached to one of the points
///in uvuv_list.
///The function's return value is:
///  =0: Intersection was not obtained.
/// !=0: Intersection was obtained as follows:
///  =1: End point is a point on a perimeter of one of the surfaces.
///  =3: End point is one of boundary points in uvuv_list.
///  =4: End point is the starting point.
///  =7: isect_start halted the computation since intersection was lost
///     during the computation.
int isect_startPt(
	const MGPosition& uvuv_startIn, ///<Starting point of the intersection line.
	MGPosition_list& uvuv_list,	///<isect_startPt will halt when ip reached one of 
		///<the point in uvuv_list. isect_startPt does not change uvuv_list(actually
		///<uvuv_list is const.) uvuv's space dimension is at least 4,
		///<and the first 2 is (u,v) of this and the next 2 is (u,v) of sf2. 
	const MGFSurface& sf2,	///<2nd surface.
	double acuRatio,///<Accurate ratio, should be decreased by multiplyng .2
		///<(or a number less than 1.).
	MGBPointSeq& point,		///<Surface-surface intersection parameter values
		///<will be returned as:point(.,0) and point(.,1) for(u,v) of this surface
		///<point(.,2) and point(.,3) for(u,v) of this surface.
		///<point has the dimension of(.,7).
	MGPosition_list::iterator& uvuv_id,
			///<When the end point of ip was one of the points of uvuv_list,
			///<uvuv_list's iterator of the point will be returned, that is, 
			///<when the function's return value was 3 or 5,
			///<When was not a point of uvuv_list, end() of uvuv_list will be
			///<returned.
	int& m1	///<id that indicates which surface was used as the main surface,
			///<m1=0: this surface, m1=2: sf2.
)const;

///Compute an intersection line of this surface(a surface that is converted to
///1D SBRep surface(sf1d)) and a plane pl.
///sf1D is so converted that pl be x=0. plane. This surface is the original
/// surface and pl is the original plane.
///The function's return value is:
/// =0: Intersection was not obtained.
/// !=0: Intersection was obtained as follows:
///    =1: End point is a point on a perimeter of one of the surfaces.
///    =3: End point is one of boundary points in uvuv_list.
///    =4: End point is the starting point.
///    =7: isect_start halted the computation since intersection was lost
///     during the computation.
int isect_startPlane(
	const MGPosition& uvuvS,///<Starting parameter value of the intersection.
	MGPosition_list& uvuv_list,
		///<isect_startPlane will halt when ip reached one of the point in
		///<uvuv_list. isect_startPlane does not change uvuv_list(actually
		///<uvuv_list is const.) uvuv's space dimension is at least 4,
		///<and the first 2 is (u,v) of this and the next 2 is (u,v) of pl. 
		///<When uvuv's space dimension is more than 4, it indicates that the uvuv
		///<is used to input approcimate tangent of the intersection.
	const MGPlane& pl,	///<Plane expression(2nd surface for the intersection).
	MGSSisect& ssi,		///<Surface-surface intersection line will be output.
	MGPosition_list::iterator& uvuv_id
			///<When the end point of ip was one of the points of uvuv_list,
			///<uvuv_list's iterator of the point will be returned, that is, 
			///<when the function's return value was 3 or 5.
			///<When was not a point of uvuv_list, end() of uvuv_list will be
			///<returned.
)const;

///isect_startPlanePt compute an array of parameter value pairs of this surf and pl2
///for one intersection line of a surface and a plane,
/// given starting intersetion point uvuv((u1,v1) of this and (u2,v2) of pl2)
/// and direction in uvuv_startI(4-6) (optionally).
///isect_startPlanePt is a dedicated function for isect_startPlane.
///isect_startPlanePt halts the computation when intersection
///reached to a boundary of this, or reached to one of the points
///in uvuv_list.
///The function's return value is:
///  =0: Intersection was not obtained.
/// !=0: Intersection was obtained as follows:
///  =1: End point is a point on a perimeter of this surfaces.
///  =3: End point is one of boundary points in uvuv_list.
///  =4: End point is the starting point.
///  =7: isect_startPlanePt halted the computation since intersection was lost
///     during the computation.
int isect_startPlanePt(
	const MGPosition& uvuv_startIn, ///<Starting point of the intersection line.
	MGPosition_list& uvuv_list,	///<isect_startPlanePt will halt when ip reached one of 
		///<the point in uvuv_list. isect_startPlanePt does not change uvuv_list(actually
		///<uvuv_list is const.) uvuv's space dimension is at least 4,
		///<and the first 2 is (u,v) of this and the next 2 is (u,v) of pl2. 
	const MGPlane& pl2,	///<2nd surface(MGPlane).
	double acuRatio,///<Accurate ratio, should be decreased by multiplyng .2
		///<(or a number less than 1.).
	MGBPointSeq& point,		///<Surface-surface intersection parameter values
		///<will be returned as:point(.,0) and point(.,1) for(u,v) of this surface
		///<point(.,2) and point(.,3) for(u,v) of pl2.
		///<point will have the dimension of(.,7).
	MGPosition_list::iterator& uvuv_id
			///<When the end point of ip was one of the points of uvuv_list,
			///<uvuv_list's iterator of the point will be returned, that is, 
			///<when the function's return value was 3.
			///<When was not a point of uvuv_list, end() of uvuv_list will be
			///<returned.
)const;

///Compute the intersection lines of this surface and srf2(both are not planes).
MGSSisects isect_with_surf(
	MGPosition_list& uvuv_list,
	///<Let a member of uvuv_list be uvuv. Then, uvuv's space dimension is
	///<at least 4, and the first 2 is (u,v) of this and the next 2 is (u,v) of srf2. 
	///<When uvuv's space dimension is more than 4, it indicates that the uvuv
	///<is used to input approximate tangent of the intersection.
	const MGFSurface& srf2	///<2nd surface for the intersection.
)const;

///Compute the intersection lines of this(is not a plane) surface and a plane pl.
///sd1D that is converted to surf1D() about pl is necessary to input.
MGSSisects isect_with_plane(
	MGPosition_list& uvuv_list,
	///<Let a member of uvuv_list be uvuv. Then, uvuv's space dimension is
	///<at least 4, and the first 2 is (u,v) of this and the next 2 is (u,v) of pl. 
	///<When uvuv's space dimension is more than 4, it indicates that the uvuv
	///<is used to input approximate tangent of the intersection.
	const MGPlane& pl,	///<Plane expression(2nd surface for the intersection).
	const MGFSurface& fsrf2	///<Original FSurface before casting to plane.
)const;

///Obtain the projected curve of a curve onto the surface.

///The direction of the projection is along the vector vec if the vec is not
///NULL, and normal to the surface if the vec is NULL.
///Output of 'project' is two kind of curves:
///one is general world coordinate curves('vec_crv'), and the other is (u,v) curves of
///the parameter space of the surfaces(vec_crv_uv).
///vec_crv_uv.size() is equal to vec_crv.size(). Let the size be n, then
/// (vec_crv_uv[i], vec_crv[i]) is one pair for 0<=i<n.
///戻り値：
///		投影曲線の数:		投影曲線が求まった
///		0:			投影曲線が求まらなかった
///		-1:			内部処理エラー
///		-2:			収束処理エラー（収束しなかった）
int projectbyApproximateAsLBRep(
	const MGCurve& crv,
	std::vector<UniqueCurve>& vec_crv_uv,
		//Projected curve(surface parameter (u,v) representation) will be appended.
	std::vector<UniqueCurve>& vec_crv,
		//Projected curve(world coordinate(x,y,z) representation) will be appended.
	const MGVector& vec
)const;

///Obtain the projected curve of a curve onto the surface.

///The direction of the projection is along the vector vec if the vec is not
///NULL, and normal to the surface if the vec is NULL.
///Output of 'project' is two kind of curves:
///one is general world coordinate curves('vec_crv'), and the other is (u,v) curves of
///the parameter space of the surfaces(vec_crv_uv).
///vec_crv_uv.size() is equal to vec_crv.size(). Let the size be n, then
/// (vec_crv_uv[i], vec_crv[i]) is one pair for 0<=i<n.
///戻り値：
///		投影曲線の数:		投影曲線が求まった
///		0:			投影曲線が求まらなかった
///		-1:			内部処理エラー
///		-2:			収束処理エラー（収束しなかった）
int projectbyRemovKnots(
	const MGCurve& crv,
	std::vector<UniqueCurve>& vec_crv_uv,
		//Projected curve(surface parameter (u,v) representation) will be appended.
	std::vector<UniqueCurve>& vec_crv,
		//Projected curve(world coordinate(x,y,z) representation) will be appended.
	const MGVector& vec
)const;

///カーブを折れで分割して行い、後で接続する
int projVector(
	const MGCurve& crv,
	std::vector<UniqueCurve>& vec_crv_uv,
	std::vector<UniqueCurve>& vec_crv,
	const MGVector& vec
) const;

///カーブを折れで分割して行い、後で接続する
int projNormal(
	const MGCurve& crv,
	std::vector<UniqueCurve>& vec_crv_uv,
	std::vector<UniqueCurve>& vec_crv
)const;

private:

///Test in_range if this is a face, if not, do nothing.
///This is to accelerate the test of in_range in isect_guess().
///See isect_guess(MGCurve).
virtual bool in_range_face(const MGPosition& uv)const{return true;};

///uvuvE_is_a_midpoint() is used by isect_with_surf() to determine if
///uvuvE is a mid point at the intersectio and should be neglected.
bool uvuvE_is_a_midpoint(
	const MGFSurface& f2,
	int m1,		///<id of uvuvE that indicates this surface's parameter
		///<position in uvuvE. (uvuvE(m1), uvuvE(m1+1))=(u,v) of this surface.
	const MGSSisect& ssi, ///<ssi obtained so far.
	const MGPosition& uvuvE///<start parameter (u,v) pair of this surface and sf2.
)const;

///Obtain the projected curve of a curve onto the surface.
///projNormalProc does not divide the curve at the c0 continuity, which is treated by
///project(). See projec().
///The direction of the projection is normal to the surface.
///Output of 'projNormalProc' is curves of space dimension 5, which are (x,y,z,u,v).
///Here (x,y,z) is the world coordinates, and (u,v) is parameter of this surface.
void projNormalProc(
	const MGCurve& crv,	///<The target curve to project.
	std::vector<UniqueLBRep>& crv_xyzuv_vector///<Projected curve of(x,y,z,u,v) will be appended.						
)const;

///ベクトル投影は、crvをvec方向にスイープした面と元の面との交線で求める
int projVectorProc(
	const MGCurve& crv,
	std::vector<UniqueLBRep>& vec_crv,	///<５次元曲線が返る
	const MGVector& vec
) const;

};

///@cond
typedef std::pair<const MGBox*,int> mgBCPair;
//Private class to sort boxes by their distances from a position.
class MG_DLL_DECLR MGBox2PositionCompare{
public:
	MGBox2PositionCompare(const MGPosition& P):m_position(P){;};
	bool operator()(const MGBox* b1, const MGBox* b2)const;
	bool operator()(const mgBCPair& b1, const mgBCPair& b2)const;
	const MGPosition& m_position;
};
///@endcond

/** @} */ // end of MGObjectRelated group
#endif
