/********************************************************************/
/* Copyright (c) 2019 System fugen G.K. and Yuzi Mizuno          */
/* All rights reserved.                                             */
/* Granted under the MIT license (see mg/MGCL.h for detail)         */
/********************************************************************/

#include "StdAfx.h"
#include <iomanip>
#include "mg/Box.h"
#include "mg/Curve.h"
#include "mg/Tolerance.h"
#include "mg/CParam_list.h"
#include "mg/SurfCurve.h"
#include "topo/LEPoint.h"
#include "topo/CPointsVec.h"
#include "topo/LCisect_vector.h"
#include "topo/Loop.h"
#include "topo/Edge.h"
#include "topo/Face.h"
#include "Tl2/TL2Polyline.h"
#include "Tl2/TL2Face.h"

#if defined(_DEBUG)
#define new DEBUG_NEW
#undef THIS_FILE
static char THIS_FILE[] = __FILE__;
#endif

typedef std::pair<int, double> SplitParam;//defines split parameter of a loop.
//SplitParam.first=0: second is u-value, 1: second is v-value, -1:undefined.
//SplitParam.second=parameter value.

//Split the face giving bridges over loops. Splitting is done by finding the smallest
//closed areas out of bridges.
void split_by_bridges(
	const MGFace& face,
	std::vector<UniqueTL2PlBridge>& bridges,
	std::vector<UniqueFace>& faces//Result trimmed face(s) will be appended.
){
	double rzero=MGTolerance::rc_zero();
	const MGBox& bxuv=face.box_param();
	double minimum_loop_area=bxuv[0].length();
	minimum_loop_area+=bxuv[1].length().value();
	minimum_loop_area*=rzero;

	mgCPointsVec cpointsVec(face);
		//Found loops that are connected to
		//an inner loop of loop id loop_id will be stored
		//in this cpointsVec[loop_id].

	int nbridges=(int)bridges.size();
	for(int ii=nbridges-1; ii>=0; ii--)
		cpointsVec.build_CPointsVec(*(bridges[ii]));
	int nboundaries=face.number_of_boundaries();
	for(int i=0; i<nboundaries; i++){
		mgCPoints& lveci=cpointsVec[i];
		lveci.sort();
	}

	for(int i=0; i<nboundaries; i++){
		mgCPoints& lveci=cpointsVec[i];
		int nlv=lveci.size();
		assert((nlv%2)==0);//nlv must be even(pair of coming_in and going_out).
		if(!nlv)
			continue;

		for(int j=0; j<nlv; j++){
			std::unique_ptr<MGFace> nface2(new MGFace);
			std::unique_ptr<MGLoop> loop;
			std::deque<MGTrimLoop*> tloops;
			mgCONNECT_POINT& fpoint=lveci[j];
			int extract_inf=cpointsVec.extract_loop(fpoint,loop,tloops);
			if(!extract_inf)
				continue;
		
			std::vector<bool> used_loops(nboundaries,false);
			size_t ntloop=tloops.size();
			for(size_t k=0; k<ntloop; k++){
				tloops[k]->set_used_loop_flag(used_loops);
				tloops[k]->set_null();
			}

			double area=loop->area();
			if(area<=minimum_loop_area && area>= -minimum_loop_area)
				continue;

			nface2->add_boundary(loop.release());
			for(int j=nboundaries-1; j>=0; j--){
				if(!used_loops[j]){
					const UniqueLoop& lpj=face.loop(j);
					MGPosition uvS=lpj->start_point();
					if(nface2->in_range(uvS))
						nface2->add_boundary(face.loop(j)->clone());
				}
			}
			faces.emplace_back(nface2.release());
		}
	}
}

//Extract MGTrimLoop & mgCPoints info, and build this mgCPointsVec data.
void mgCPointsVec::build_CPointsVec(
	mgTL2PlBridge& bridge	//original loops to trim.
							//includes uv is output when return code =0 or 1.
){
	for(int i=0; i<2; i++){
		MGLoop* lp=new MGLoop;
		MGLEPoint leS, leE;
		MGEdge* edg;
		if(i==0){
			leS=bridge.le_start(), leE=bridge.le_end();
			edg=bridge.copy_edge();
		}else{
			leE=bridge.le_start(), leS=bridge.le_end();
			edg=bridge.release_edge();
			edg->negate();
		}
		lp->append(edg);

		int lidS=leS.loop()->get_loop_id_in_face();
		int lidE=leE.loop()->get_loop_id_in_face();
		MGTrimLoop* trloop=new MGTrimLoop(lp,-lidS,leS,-lidE,leE);
		m_trim_loops.emplace_back(trloop);
		push_at(lidS,mgCONNECT_POINT(trloop,mgCONNECT_POINT::going_out));
		push_at(lidE,mgCONNECT_POINT(trloop,mgCONNECT_POINT::coming_in));
	}
}

//Compute all the intersections of face boundaries and sl, which are all parameter
//representation.
MGLCisect_vector isectBoundary(
	double error,//error to get intersection of the loops and sl.
	const MGFace& face,//face to split
	const MGStraight& sl
){
	MGLCisect_vector pvec=face.loop(0)->isectLoop(error,sl);
	int n_inner=face.number_of_inner_boundaries();
	for(int k=0; k<n_inner; k++){
		const UniqueLoop& innerLoopk=face.loop(k);
		pvec.append(innerLoopk->isectLoop(error,sl));
	}
	return pvec;
}

//Compute the intersections of coordinate f (u or v value of
//the surface parameter representation) and loop lp.
void isect1DTl(
	const mgTL2parameter& tlparam,
	int coordinate,	// f coordinate kind, =0: u, =1:v
	double f,
	const MGLoop& lp,
	std::vector<MGLEPoint>& pvec, //intersections be appended.
	std::pair<MGLEPoint,bool>* limit
		//when limit->second==true, isect->first be employed as the isect.
){
	double error = coordinate ? tlparam.get_VError() : tlparam.get_UError(); 
	mgTolSetWCZero wczeroSet(error);//Set&save the error.
	
	double flimit=0.;
	bool employ_greater=true;
	int coordinate2=(coordinate+1)%2;
	if(limit){
		MGLEPoint uvLE=limit->first;
		pvec.push_back(uvLE);
		employ_greater=limit->second;
		MGPosition uv=uvLE.eval();
		flimit=uv[coordinate2];
		if(employ_greater)
			flimit+=error*10.;
		else
			flimit-=error*10.;
	}
	MGComplex::const_iterator ei=lp.pcell_begin(), ee=lp.pcell_end();
	for(; ei!=ee; ei++){
		const mgTL2Polyline* ipoly=TL2Polyline(ei);
		if(limit){
			const MGInterval& frange=ipoly->box()[coordinate2];
			if(employ_greater){
				if(frange<flimit)
					continue;
			}else{
				if(frange>flimit)
					continue;
			}
		}

		int nbd=ipoly->bdim();
		short idLast=short(nbd-1);
		MGCParam_list tlist=ipoly->isect_1D(f, coordinate);
		MGCParam_list::iterator ti=tlist.begin(), te=tlist.end();
		MGComplex::const_iterator ei2=ei;
		double ts=ipoly->param_s();

		for(;ti!=te; ti++){
			double  t=*ti;
			if(limit){
				MGPosition uv=ipoly->eval(t);
				double ft=uv[coordinate2];
				if(employ_greater){
					if(ft<flimit)
						continue;
				}else{
					if(ft>flimit)
						continue;
				}
			}
			//change the parameter value to indicate a vetex point of the polygon.
			short idt=short(t-ts+.5);
			MGLEPoint le(ei,ts+double(idt));
			if(idt==idLast)
				le=le.start_of_next_edge();
			size_t npvec=pvec.size();
			size_t i=0;
			for(; i<npvec; i++)
				if(pvec[i]==le)
					break;
			if(i==npvec)
				pvec.push_back(le);
		}
	}
}

class pvecComp{
	int m_kcod;
public:
	pvecComp(int kcod):m_kcod(kcod){;};
	bool operator()(const MGLEPoint& le1, const MGLEPoint& le2);
};
bool pvecComp::operator()(const MGLEPoint& le1, const MGLEPoint& le2){
	MGPosition uv1=le1.eval(), uv2=le2.eval();
	return uv1[m_kcod]<uv2[m_kcod];
}

//Compute minimum sine angle of crv and the edge at pv[] in their world representation.
//When sin angle is small enough crv is parallel to loop(boundary of the face) at
//the start or end point.
//Function's return value is min of sin angle at both end.
double angle_edge_crv(
	const MGSurface& srf,
	const MGCurve& crv,//(u,v) parameter line representation
	std::vector<MGLEPoint>::const_iterator pv
		//a point on the loop of the face at start and end of crv.
){
	MGSurfCurve crvOnsrf(srf,crv);//Build world curve rep.
	double t[2]={crvOnsrf.param_s(), crvOnsrf.param_e()};
	double angmin=3., ang1;
	for(int i=0; i<2; i++){//test for crv's start and end point.
		MGVector crvDir=crvOnsrf.eval(t[i],1);//crv's direction
		if(i)
			crvDir*=-1.;//When end point of crvOnsrf(to direct to inner of the loop)

		const MGLEPoint& pvi=pv[i];
		const MGEdge* edgi=pvi.edge();
		const mgTL2Polyline* edgeiCrv=TL2Polyline(edgi);
		if(pvi.is_Edge_start_point()){
			MGVector v0=edgeiCrv->direWorld_start();//edge start point direction.
			ang1= parallelismTl(crvDir,v0);

			const MGEdge* epre=edgi->pre_edge();
			const mgTL2Polyline* eCrvPre=TL2Polyline(epre);
			MGVector Vpre=eCrvPre->direWorld_end();//previous edge end point direction.
			Vpre*=-1.;//To direct away from the vertex.
			double ang2=parallelismTl(crvDir,Vpre);
			if(ang2<ang1)
				ang1=ang2;
		}else{
			MGVector v0=edgeiCrv->direWorld_with_non_normalized(pvi.param());
			ang1= parallelismTl(crvDir,v0);
		}
		if(ang1<angmin)
			angmin=ang1;
	}
	return angmin;
}

//Get_split_network.
//Function's return value is true, if network is obtained, false if not.
//False means the parameter value t_to_divide was not adequate for the split.
bool get_split_bridge(
	int tryNum,//input how many try number this is, from 1.
	SplitParam para[2],
	const MGFace& face,//Target face to split.
	const mgTL2parameter& tlparam,
	std::vector<UniqueTL2PlBridge>& bridges,
	std::pair<MGLEPoint,bool>* limit
){
	const MGSurface& sufaceOrg=tlparam.get_surface();
	double serror=sufaceOrg.parameter_error();
	bool obtained=false;
	int nloop=face.number_of_loops();
	double zeroAngle=LOOSE_ZERO_ANGLE/double(tryNum);

	for(int numtry=0; numtry<2; numtry++){

	if(para[numtry].first==-1)
		break;
	int kcod0= para[numtry].first;
	double& t_to_divide=para[numtry].second;

	int kcod1= (kcod0+1)%2;
	std::vector<MGLEPoint> pvec;
	for(int k=0; k<nloop; k++)
		isect1DTl(tlparam,kcod0,t_to_divide,*(face.loop(k)),pvec,limit);
	int nisecm1=(int)pvec.size()-1;
	if(nisecm1>1)
		std::sort(pvec.begin(),pvec.end(), pvecComp(kcod1));

	auto pveci=pvec.begin();
	for(int i=0;i<nisecm1; i++, pveci++){
		MGLEPoint& pvi=pvec[i]; MGLEPoint& pvip1=pvec[i+1];
		MGPosition uv0=pvi.eval(), uv1=pvip1.eval();
		MGPosition uvmid=(uv0+uv1)*.5;
		int in=face.in_range_with_on(uvmid);
		if(in!=2)
			continue;

		MGStraight sl(uv1,uv0);
		if(angle_edge_crv(sufaceOrg,sl,pveci)<=zeroAngle)//When nearly parallel
			break;

		const mgTL2Polyline* poly=TL2Polyline(pvi.edge());
		if(int(pvi.param()-poly->param_s()+.5)==poly->bdim()-2)
			if(pvi.start_of_next_edge()==pvip1)//If pvi and pvipi are the same.
			continue;

		MGLCisect_vector tempI=isectBoundary(serror,face,sl);
		if(tempI.entries()>2)
			continue;

		mgTL2Polyline* uvpolyline=new mgTL2Polyline(pvip1,pvi);
		mgTL2PlBridge* bridge=new mgTL2PlBridge(uvpolyline,pvi,pvip1);
		bridges.emplace_back(bridge);
		obtained=true;
	}
	if(obtained)
		break;
	else
		bridges.clear();

	}
	return obtained;
}

double normalized_dif_angle(double angle,double base_angle){
	double dif1=angle-base_angle;
	double dif2=dif1;
	if(dif1<0.)
		dif2+=mgDBLPAI;
	else
		dif1-=mgDBLPAI;
	if(fabs(dif1)>fabs(dif2))
		dif1=dif2;
	return dif1;
}

//Find concave vertex of the outer loop lp.
//Function's return value is true if concave vertex is found.
bool find_concave_at_vertex(
	const MGLoop& lp,//Target loop.
	MGLEPoint& le_vertex,//MGLEPoint of the concave vertex be returned, which is always
					//the starting vertex of an edge.
	bool& is_u,	//split parameter kind, =true if u, false if v.
		//Original kind is input, and is modified fi necessary.
	SplitParam& para,
	bool& employ_greater
){
	int nedges=lp.number_of_edges();
	const MGEdge* pre=lp.last_edge();
	const mgTL2Polyline* crvPre=TL2Polyline(pre);
	const MGSurface& srf=crvPre->TL2param().get_surface();

	double angmin=mgDBLPAI;
	MGVector Nsave;
	MGComplex::const_iterator i=lp.pcell_begin(), ie=lp.pcell_end();
	for(; i!=ie; i++){
		const MGVector& Vpre=crvPre->direWorld_end();//

		const MGEdge* aft=edge_from_iterator(i);
		const mgTL2Polyline* crvAft=TL2Polyline(aft);
		const MGVector& Vaft=crvAft->direWorld_start();//

		const MGPosition& uv_vertex=crvAft->uv_start();
		MGVector N=srf.unit_normal(uv_vertex);
		double ang=Vpre.angle2pai(Vaft,N);
		if(ang>=mgPAI)
			ang-=mgDBLPAI;
		if(ang<angmin){
			angmin=ang;
			le_vertex=MGLEPoint(i,crvAft->param_s());
			Nsave=N;
		}
		pre=aft;
		crvPre=crvAft;
	}

	if(angmin<LOOSE_ZERO_ANGLE2){
		MGPosition uv=le_vertex.eval();
		const MGEdge* aft=le_vertex.edge();
		const mgTL2Polyline* crvAft=TL2Polyline(aft);
		const MGVector& Vaft=crvAft->direWorld_start();//

		const MGEdge* pre=aft->pre_edge();
		const mgTL2Polyline* crvPre=TL2Polyline(pre);
		const MGVector& Vpre=crvPre->direWorld_end();//

		MGVector Vu=srf.eval(uv,1),Vv=srf.eval(uv,0,1);//
		double angupre=Vu.sangle(Vpre), angvpre=Vv.sangle(Vpre);
		double anguaft=Vu.sangle(Vaft), angvaft=Vv.sangle(Vaft);

		//Define is_u. We assume Vu is not parallel Vv(not singular).
		int isoIsParallelToEdge=0;
		if(angupre<=LOOSE_ZERO_ANGLE){
			isoIsParallelToEdge=1;
			if(angvaft>LOOSE_ZERO_ANGLE)
				is_u=true;//case 1
		}else{
			if(angvpre<=LOOSE_ZERO_ANGLE){
				isoIsParallelToEdge=2;
				if(anguaft>LOOSE_ZERO_ANGLE)
					is_u=false;//case 2
			}else if(anguaft<=LOOSE_ZERO_ANGLE){
				isoIsParallelToEdge=3;
				is_u=true;//case 3
			}else if(angvaft<=LOOSE_ZERO_ANGLE){
				isoIsParallelToEdge=4;
				is_u=false;//case 4
			}/*else{//case 5(Neither of Vu or Vv is parallel to Vaft or Vpre).
				if((Vv*Vpre)%(Vv*Vaft)>0.)
					is_u=true;
				else if((Vu*Vpre)%(Vu*Vaft)>0.)
					is_u=false;
			}*/
		}

		double sangVpreVaft=Vpre.sangle(Vaft);
		const MGVector& Viso=is_u ? Vv:Vu;
		if(-LOOSE_ZERO_ANGLE2<=sangVpreVaft && sangVpreVaft<=LOOSE_ZERO_ANGLE2){
			//When angle of concave is small.
			MGVector Visoaft=Viso*Vaft;
			employ_greater=Visoaft%Nsave<0.;
		}else{
			if(isoIsParallelToEdge){
				const MGVector& Vedge=isoIsParallelToEdge<=2 ? Vaft:Vpre;
				if(isoIsParallelToEdge<=2){
					employ_greater=Viso%Vaft<0.;
				}else{
					employ_greater=Viso%Vpre>0.;
				}
			}else{
				double ang=Vaft.angle2pai(Viso,Nsave);
				employ_greater=ang<=mgPAI;
			}
		}
		para.first= is_u ? 0:1;
		para.second= uv[para.first];
		return true;
	}
	return false;
}

//Find edges that constitute concavity from the outer loop lp.
//Let numEdge be function's return value.
//Then numEdge continuous edges constitute the concavity if numEdge>0.
int find_concave_over_edges(
	const MGLoop& lp,//Target face to split.
	const MGEdge*& edgeStart,//Starting edge that constitues the concavity will be returned.
		//numEdge continuous edges constitute the concavity.
	double& difAngle
){
	double wzero2=MGTolerance::wc_zero_sqr();
	int numEdge=0;
	double angle3=difAngle=0.;

	int nedges=lp.number_of_edges();
	MGComplex::const_iterator i=lp.pcell_begin(), ie=lp.pcell_end();
	for(int j=0; j<nedges; j++,i++){
		const MGEdge* ei=edge_from_iterator(i);
		const mgTL2Polyline* crvi=TL2Polyline(ei);
		const MGVector& V0=crvi->direWorld_start();//
		if(V0%V0<=wzero2){
			numEdge=0;
			difAngle=0.;
			continue;
		}

		double angle0=crvi->angle2Uline(0.);
		if(numEdge){
			double angle2NewEdge=normalized_dif_angle(angle0,angle3);
			if(angle2NewEdge>=LOOSE_ZERO_ANGLE){
				if(difAngle < LOOSE_ZERO_ANGLE*2.-mgHALFPAI){
					return numEdge;
				}else{
					numEdge=0;
					difAngle=0.;
				}
			}else
				difAngle+=angle2NewEdge;
		}
		double angle1=crvi->angle2Uline(0.3333);
		difAngle+=normalized_dif_angle(angle1,angle0);
		double angle2=crvi->angle2Uline(0.6666);
		difAngle+=normalized_dif_angle(angle2,angle1);
		angle3=crvi->angle2Uline(1.);
		difAngle+=normalized_dif_angle(angle3,angle2);
		if(difAngle>=LOOSE_ZERO_ANGLE){
			numEdge=0;
			difAngle=0.;
		}else{
			if(!numEdge)
				edgeStart=ei;
			numEdge++;
		}
	}
	if(difAngle >= LOOSE_ZERO_ANGLE-mgHALFPAI)
		numEdge=0;
	return numEdge;
}

//Get the split u or v value and the parameter limit at uv of a concave edge.
void compute_split_limit(
	const MGPosition& uv, //surface parameter value at a concave edge.
	const MGSurface& srf,
	const MGVector& Ve,	//direction of the edge at uv in world representation.
	bool is_u,	//split parameter kind, =true if u, false if v.
		//Original kind is input, and may be modified.
	SplitParam& para,
	bool& employ_greater
){
	MGVector Vu=srf.eval(uv,1);
	MGVector Vv=srf.eval(uv,0,1);

	double uangl=Ve.sangle(Vu);
	double vangl=Ve.sangle(Vv);
	if(uangl>vangl*2.)
		is_u=false;
	if(vangl>uangl*2.)
		is_u=true;
	const MGVector& Viso=is_u ? Vv:Vu;
	MGVector N=Vu*Vv;
	MGVector NVe=N*Ve;
	double direction=NVe%Viso;
	employ_greater=direction>0.;

	para.first= is_u ? 0:1;
	para.second= uv[para.first];
}

//Find the 1st vertex of an edge that exeeds the half of concave angle difAngle.
//The search is started from the 1st vertex of concaveEdge and is repeated
//to the numE edges.
int find_uv_to_split(
	int numE,	//Number of edges to search from the edge concaveEdge.
	const MGEdge* concaveEdge,//1st edge that constitutes the concavity. numE edges constitute
				//the total concavity that have the concavity angle difAngle.
	double difAngle,//Total concavity angle.
	const mgTL2Polyline*& crv,//uv at which vertex split should be done will be returned.
	MGLEPoint& uv_le//MGLEPoint data of uv.
){
	MGComplex::const_iterator i=concaveEdge->edge_iterator();
	double halfA=difAngle*.5;
	double  angleSummed=0., angle0=0.;
	for(int j=0; j<numE; j++,i++){
		crv=TL2Polyline(i);
		int nv=crv->number_of_points();
		int k;
		for(k=0; k<nv; k++){
			double anglek=crv->angle2Uline_at_i(k);
			if(j || k){
				angleSummed+=normalized_dif_angle(anglek,angle0);
				if(angleSummed<=halfA)//When the concavity angle exceeds the half.
					break;
			}
			angle0=anglek;
		}
		if(angleSummed<=halfA){//When the concavity angle exceeds the half.
			int nv3=nv/3;
			if(k<=nv3 && j){
				k=0;
			}else if((nv-nv3)<=k && j<(numE-1)){
				k=0;
				crv=TL2Polyline(++i);
			}
			double t=crv->param_s()+double(k);
			uv_le=MGLEPoint(i,t);
			return k;
		}
	}
	return -1;//This should not happen.
}

//Find at what parameter the loop be subdivided for the tessellation.
//Candidate params are stored in para.
void find_where(
	int tryNum,//input how many try number this is, from 1.
	int loopID,	//=0 if lp is outer, >=1 if inner.
	const MGLoop& lp,//Input loop.
	bool& is_u,	//Valid only when lp is outer, and
		//indicates if split is done at u paramete(is_u=true) or not.
	SplitParam para[2],//[0] is 1st candidate, [1] is 2nd.
	double* t_to_avoid//parameter to avoid. if input, employ other value.
){
	para[0].first=para[1].first=-1;//set undefined.
	const MGBox& uvrange=lp.box();
	double uspan=uvrange[0].length();
	double vspan=uvrange[1].length();
	if(loopID&&tryNum<=1){//When inner loop.
		if(uspan>=2.*vspan)
			is_u=true;
		else if(vspan>=2.*uspan)
			is_u=false;
	}
	int id=is_u ? 0:1;
	const MGInterval& trange=uvrange[id];
	double tspan=trange.length();
	double terror=tspan*NEAR_PARAM;
	double deviMaxFromMid=tspan*MAX_DEVIATION_FROM_MIDDLE;//maximum deviation.
	double mid=trange.mid_point();//mid point of trange.
	double lengthMinFromMid, lengthMin2FromMid;
	double& t0=para[0].second;//target parameter value stored(1st candidate).
	double& t1=para[1].second;//target parameter value stored(2nd candidate).
	MGComplex::const_iterator i=lp.pcell_begin(), ie=lp.pcell_end();
	for(; i!=ie; i++){
		const mgTL2Polyline* pli=TL2Polyline(i);
		MGPosition uvi=pli->uv(0);
		double ti=uvi[id];
		double leni=(ti<=mid ? mid-ti:ti-mid);
		if(leni>deviMaxFromMid)//If ti is far away from mid.
			continue;
		if(t_to_avoid && *t_to_avoid-terror<=ti && ti<=*t_to_avoid+terror)
			continue;

		if(para[0].first==-1){
			para[0].first=id, t0=ti; lengthMinFromMid=leni;
		}else if(para[1].first==-1){
			if(leni<lengthMinFromMid){
				para[1]=para[0];lengthMin2FromMid=lengthMinFromMid;
				para[0].first=id, t0=ti; lengthMinFromMid=leni;
			}else{
				para[1].first=id, t1=ti; lengthMin2FromMid=leni;
			}
		}else{
			if(leni<lengthMinFromMid){
				para[1]=para[0]; lengthMin2FromMid=lengthMinFromMid;
				para[0].first=id, t0=ti; lengthMinFromMid=leni;
			}else if(leni<lengthMin2FromMid){
				para[1].first=id, t1=ti; lengthMin2FromMid=lengthMinFromMid;
			}
		}
	}
	if(para[0].first==-1){
		if(tryNum==1){
			para[0].first=id, t0=mid;
			para[1].first=id, t1=mid+deviMaxFromMid*.5;
		}else{
			para[0].first=id, t0=mid-deviMaxFromMid*.5;
			para[1].first=id, t1=mid+deviMaxFromMid*.8;
		}
	}
}

//Split face.
//Function's return value is true, if split was executed.
//false, if not.
bool splitTl(
	const MGFace& face,//Target face to split.
	const mgTL2parameter& tlparam,
	bool is_u,	//indicates if split be done in u or v direction.
			//If true split is done around the middle of u parameter.
			//This is valid only when face has no inner loops.
	std::vector<UniqueFace>& faces//All of the splitted faces will be output.
){
	//Define at what parameter face be split.

	int nloop=face.number_of_loops();
	int loopID=(nloop>1) ? 1:0;
	const MGLoop& lp=*(face.loop(loopID));

	SplitParam para[2];//[0] is 1st, [1] is 2nd candidate to split.
	para[0].first=is_u?0:1;
	double& t_to_divide=para[0].second;
	para[1].first=-1;

	std::pair<MGLEPoint,bool> limitData;
	MGLEPoint& uv_le=limitData.first;
	bool& employ_greater=limitData.second;

	std::vector<UniqueTL2PlBridge> bridges;
	bool obtained=false;
	if(loopID==0){
		if(find_concave_at_vertex(lp,uv_le,is_u,para[0],employ_greater)){
			//If concave found at a vertex.
			obtained=get_split_bridge(1,para,face,tlparam,bridges,&limitData);
		}
		if(!obtained){
			const MGEdge* concaveEdge;
			double difAngle;
			int numE=find_concave_over_edges(lp,concaveEdge,difAngle);
			if(numE){//If concave found over edges.
				const mgTL2Polyline* crv;
				int idV=find_uv_to_split(numE,concaveEdge,difAngle,crv,uv_le); assert(idV>=0);
				MGPosition uv=crv->uv(idV);
				MGVector Ve=crv->direWorld_with_non_normalized(uv_le.param());
				const MGSurface& srf=tlparam.get_surface();
				compute_split_limit(uv,srf,Ve,is_u,para[0],limitData.second);
				obtained=get_split_bridge(1,para,face,tlparam,bridges,&limitData);
			}
		}
	}

	if(!obtained){
		int tryNum=1;
		find_where(tryNum,loopID,lp,is_u,para,nullptr);
		if(!get_split_bridge(tryNum++,para,face,tlparam,bridges,nullptr)){
			//1st try of find_where.
			double tOld=t_to_divide;
			find_where(tryNum,loopID,lp,is_u,para,&tOld);
			if(!get_split_bridge(tryNum++,para,face,tlparam,bridges,nullptr)){
				//2nd try of find_where.
				is_u=!is_u;
				tryNum=1;
				find_where(tryNum,loopID,lp,is_u,para,nullptr);
				if(!get_split_bridge(tryNum++,para,face,tlparam,bridges,nullptr)){
					//3rd try of find_where.
					tOld=t_to_divide;
					find_where(tryNum,loopID,lp,is_u,para,&tOld);
					get_split_bridge(tryNum++,para,face,tlparam,bridges,nullptr);
					//final try of find_where.
				}
			}
		}
	}

	bool splitted=false;
	if(bridges.size()){
		split_by_bridges(face,bridges,faces);
		splitted=(faces.size()>1);
	}
	return splitted;
}

mgTL2PlBridge::mgTL2PlBridge(
	mgTL2Polyline* line,//newed object, whose ownership is transfered to this.
	const MGLEPoint& le_start,//line's start point's MGLEPoint.
	const MGLEPoint& le_end//line's end point's MGLEPoint.
):m_polyLine(line), m_le_start(le_start), m_le_end(le_end){;}

//Make a cop of m_polyLine and make MGEdge.
//Function's return value is newed MGEdge pointer.
MGEdge* mgTL2PlBridge::copy_edge(){
	mgTL2Polyline* line=new mgTL2Polyline(*m_polyLine);
	return new MGEdge(line);
}

//Release m_polyLine object and make MGEdge.
//Function's return value is newed MGEdge pointer.
//After release_edge, m_polyLine is null.
MGEdge* mgTL2PlBridge::release_edge(){
	mgTL2Polyline* line=m_polyLine.release();
	return new MGEdge(line);
}

std::ostream& operator<< (std::ostream& out, const mgTL2PlBridge& bridge){
	out<<"mgTL2PlBridge:"<<*(bridge.m_polyLine);
	out<<"start="<<bridge.m_le_start<<", end="<<bridge.m_le_end
		<<std::endl<<std::endl;
	return out;
}
