/********************************************************************/
/* Copyright (c) 2017 System fugen G.K. and Yuzi Mizuno          */
/* All rights reserved.                                             */
/********************************************************************/
#include "MGCLStdAfx.h"
#include "mg/Box.h"
#include "mg/BPointSeq.h"
#include "mg/Position.h"
#include "mg/Curve.h"
#include "mg/LBRep.h"
#include "mg/RLBRep.h"
#include "mg/Tolerance.h"
using namespace std;

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

// Implementation of knot rebuild.

//曲線列のノットベクトルを再構築する
//入力された複数曲線列を指定オーダーで再構築する。トレランスはline_zero()を使用している。
//オーダーが指定されていないとき曲線列のうちで最も大きいオーダーを使用する。このとき、
//Ellipse, Straightのオーダーは4として考える。
//パラメータ範囲は1次微分値の大きさが１になるようにしたときの長さの平均を使用している。
//戻り値は再構築後の曲線列が返却される。
MGPvector<MGLBRep> rebuild_knot(
	const std::vector<const MGCurve*>& crvl,	//入力曲線列
	int ordr ,			//指定オーダー
	MGLBRep**		tp	//接続面	input and output.
		//if tp[i] for crvl[i] was not null, converted new tp will be output.
){
	//使用するオーダーと空間次元を求める
	int nCrv = (int)crvl.size(), maxOrder = 0, maxSdim = 0, i = 0;
	int idMaxOrder=0;
	MGPvector<MGLBRep> rtnBrepList(nCrv);
	double allParam = 0.0;
	for(i = 0; i < nCrv; i++){
		int ord = crvl[i]->order(), sdim = crvl[i]->sdim();
		//if(sdim < 1){
		//	rtnBrepList.clear();
		//	return rtnBrepList;	//空間次元0はエラーとする
		//}
		if(ord > maxOrder){
			maxOrder = ord; idMaxOrder=i;
		}
		if(sdim > maxSdim)
			maxSdim = sdim;
	}
	if(!ordr)
		ordr = maxOrder;		//オーダーが指定されていないときの処理
	if(ordr <= 2)
		ordr = 4;	//Ellipse, Straightのときオーダー4にする

	//LBRepでないものは、LBRepに変換してリストに入れる
	//全てノットベクトルが等しいLBRepのときリビルドしないでそのままの曲線を返却する
	bool fsame = true;
	for(i = 0; i < nCrv; i++){
		const MGLBRep *ptempBRep = dynamic_cast<const MGLBRep*>(crvl[i]);
		if(ptempBRep){
			rtnBrepList.reset(i, new MGLBRep(*ptempBRep));
			if(i && fsame){
				// KnotVectorが等しいかどうかのチェック
				const MGKnotVector& kv0 = rtnBrepList[i-1]->knot_vector();
				const MGKnotVector& kv1 = rtnBrepList[i]->knot_vector();
				if (kv0.param_s() != kv1.param_s() ||
					kv0.param_e() != kv1.param_e() ||
					kv0 != kv1) fsame = false;  
//				if(rtnBrepList[i]->knot_vector() != rtnBrepList[i - 1]->knot_vector())fsame = false;
			}
		}else{
			rtnBrepList.reset(i, new MGLBRep(*crvl[i], ordr));
			fsame = false;
		}
		//1次微分値の大きさが１になるようなパラメータ範囲を求める
		double t0=rtnBrepList[i]->param_s(), t1=rtnBrepList[i]->param_e();
		double t2=(t0+t1)*.5;
		double mag=rtnBrepList[i]->eval(t0,1).len();
		mag+=rtnBrepList[i]->eval(t1,1).len();
		mag+=rtnBrepList[i]->eval(t2,1).len();
		mag /= 3.;
		allParam += rtnBrepList[i]->param_span() * mag;
	}
	if(fsame && tp==0)
		return rtnBrepList;

	double avgSpan = allParam / nCrv;

	//パラメータ範囲を合わせたノットベクトルを作成する。
	for(i = 0; i< nCrv; i++){
		rtnBrepList[i]->change_range(0.0, avgSpan);
		if(tp){
			if(tp[i])
				tp[i]->change_range(0.0, avgSpan);
		}
	}

	//全てのノットベクトルを足し合わせる。
	MGNDDArray tau_temp;
	tau_temp.update_from_knot(rtnBrepList[idMaxOrder]->knot_vector());
	if(tau_temp.length()<ordr)
		tau_temp.change_number(ordr);
	MGKnotVector mixKnotVector(tau_temp,ordr);
	for(i = 1; i < nCrv; i++){
		MGKnotVector& ti=rtnBrepList[i]->knot_vector();
		int orderti=ti.order();
		int orderi= ordr<=orderti ? ordr: orderti;
		mixKnotVector.mix_knot_vector(MGKnotVector(tau_temp.update_from_knot(ti),orderi));
		//if(tp){
		//	if(tp[i]) mixKnotVector.mix_knot_vector(MGKnotVector(tp[i]->knot_vector(), ordr));
		//}
	}
	//std::cout<<mixKnotVector<<endl;///////////
	//ノットベクトルを十分に増やす。offset_div_num()を使用する。
	MGKnotVector tempKnotVector = mixKnotVector;
	int bd=tempKnotVector.bdim();
	for(i=ordr-1; i<bd; i++){
		int maxNumDiv = 0, j = 0;
		double	spara = tempKnotVector(i),
				epara = tempKnotVector(i + 1);
		if(MGRZero((epara - spara) / avgSpan))continue;	//マルチノットのときの処理
		for(; j < nCrv; j++){	//各スパンの最大分割数を求める
			int ndiv = rtnBrepList[j]->offset_div_num(MGInterval(spara, epara));
			if(ndiv>maxNumDiv)
				maxNumDiv = ndiv;
		}
		double span = (epara - spara) / maxNumDiv, tParam = tempKnotVector(i);
		for(j = 0; j < maxNumDiv - 1; j++){
			tParam += span;
			mixKnotVector.add_data(tParam, ordr - 1);	//ノットベクトルを増やす
		}
	}
	//std::cout<<mixKnotVector<<endl;///////////

	//増やしたノットで曲線を再作成する
	MGNDDArray tau;
	tau.update_from_knot(mixKnotVector);
	int n = tau.length();
	double tplen=0.;
	for(i = 0; i < nCrv; i++){
		double tpleni=0.;
		MGBPointSeq bp(n, maxSdim);
		for(int j = 0; j < n; j++){
			bp.store_at(j, rtnBrepList[i]->eval(tau(j)));
		}
		rtnBrepList.reset(i, new MGLBRep(tau, bp, mixKnotVector));
		if(tp){
			if(tp[i]){
				MGBPointSeq bptp(n, maxSdim);
				for(int j = 0; j < n; j++){
					MGVector tpatj=tp[i]->eval(tau(j));
					bptp.store_at(j, tpatj);
					tpleni+=tpatj.len();
				}
				tp[i]=new MGLBRep(tau, bptp, mixKnotVector);
				tpleni/=double(n);
			}
		}
		tplen+=tpleni;
	}
	tplen/=double(nCrv);

	//曲線列のノットベクトルを精度を保持して削除を行う
	::remove_knot_curves(rtnBrepList,tp,tplen);
	//std::cout<<rtnBrepList[0]->knot_vector()<<endl;//////////////
	return rtnBrepList;
}

MGPvector<MGLBRep> rebuild_knot(
	const MGPvector<MGCurve>& brepl,	//入力曲線列
	int ordr,		//指定オーダー
	MGLBRep**	tp	//接続面
){
	size_t n=brepl.size();
	std::vector<const MGCurve*> curves(n);
	for(size_t i=0; i<n; i++)
		curves[i]=brepl[i];
	return rebuild_knot(curves,ordr,tp);
}

///Get average tangent length.
///The average means the average of start, mid, and end point's tangent.
double MGCurve::get_average_tangent_length()const{
	double ts=param_s(), te=param_e();
	double tm=(ts+te)*.5;
	MGVector v0=eval(ts,1);
	MGVector v1=eval(tm,1);
	MGVector v2=eval(te,1);
	return (v0.len()+v1.len()+v2.len())/3.;
}

///Rebuild this curve.
std::unique_ptr<MGCurve> MGCurve::rebuild(
	int how_rebuild,
		//intdicates how rebuild be done.
		// =0: no approximation(only parameter change)
		// =1: if this is rational spline(MGRLBRep), reconstructed with new knot configuration
		//     as rational spline(MGRLBRep).
		//     Otherwise approximated by non-rational spline(MGLBRep) with new knot configuration.
		// =2: approximated by non-rational spline(MGLBRep) with new knot configuration,
		//     if this is rational spline. If this is not rational spline, same as =1.
	int parameter_normalization,
		//Indicates how the parameter normalization be done:
		//=0: no parameter normalization.
		//=1: normalize to range=(0., 1.);
		//=2: normalize to make the average length of the 1st derivative 
		//    is as equal to 1. as possible.
		//=3: specify parameter range in param_range.
	double tol,	///<tolerance allowed for the approximation
		///When tol<=0., MGTolerance::line_zero() will be employed.
	int ordr,	///<order of the new MGLBRep, >=4 is recommended.
		///When order=0 is input, the original order is unchanged if this curve is
		///MGLBRep or MGRLBRep. Otherwise order is set to 4.
	const double* param_range
)const{
	if(parameter_normalization>=3 && param_range==0)
		parameter_normalization=2;

	double tol_old, err=tol;
	if(tol>0.)
		tol_old=MGTolerance::set_line_zero(tol);
	else
		err=MGTolerance::line_zero();

	int neworder=ordr;
	if(neworder==0){
		neworder=order();
	}
	int reparaType=parameter_normalization;
	int reparaType2=reparaType;
		
	double tspan[2]={0., 1.};
	double& t1=tspan[0];
	double& t2=tspan[1];
	if(reparaType>=2)
		reparaType2=1;
	if(reparaType==3){
		t1=param_range[0];
		t2=param_range[1];
		if(t1>=t2){
			t2=t1+1.;
		}
	}

	std::unique_ptr<MGCurve> curvenew;
	if(how_rebuild){
		const MGRLBRep* rlb=dynamic_cast<const MGRLBRep*>(this);
		if(rlb && how_rebuild==1){
			std::unique_ptr<MGRLBRep> rlb2=rlb->rebuild_with_new_knot_configuration(err,reparaType2);
			curvenew.reset(rlb2.release());
		}else{
			MGLBRep* lbrep=new MGLBRep;
			approximate_as_LBRep(*lbrep,neworder,reparaType2);
			curvenew.reset(lbrep);
		}
	}else{
		curvenew=std::unique_ptr<MGCurve>(clone());
		if(reparaType)
			curvenew->change_range(0.,1.);
	}
	if(reparaType==2)
		t2=curvenew->get_average_tangent_length();

	if(reparaType>=2)
		curvenew->change_range(t1,t2);

	if(tol>0.)
		MGTolerance::set_line_zero(tol_old);
	curvenew->copy_appearance(*this);
	return curvenew;
}

///Approximate this curve as a MGLBRep curve
///within the tolerance MGTolerance::line_zero().
///When parameter_normalization=0, reparameterization will not be done, and
///the evaluation at the same parameter has the same values before and after
///of approximate_as_LBRep.
void MGCurve::approximate_as_LBRep(
	MGLBRep& lb,	///<Approximated LBRep will be set.
	int ordr,		///<new order
	int parameter_normalization,
		//Indicates how the parameter normalization be done:
		//=0: no parameter normalization.
		//=1: normalize to range=(0., 1.);
		//=2: normalize to make the average length of the 1st derivative 
		//    is as equal to 1. as possible.
	bool neglectMutli
)const{
	lb.copy_appearance(*this);

	const MGKnotVector& t=knot_vector();
	double ts=param_s(), te=param_e();
	int k=t.locate(ts)+1, n=t.locate(te)+1;
		//k!=t.order() or n!=t.bdim() of t occurs when this is a trimmedcurve.
	int km1=t.order()-1, start=k;
	int index, multi_found;
	MGInterval pspan(ts,te);
	MGKnotVector& tnew=lb.knot_vector();
	int norder=ordr ? ordr:4;

	do{	//Approximation by dividing to parts of continuity>=C0.
		multi_found=t.locate_multi(start,km1,index);//Locate C0 continuity point.
		if(start==k){//For the 1st span.
			approximate_as_LBRep2(lb,norder,start-1,index,neglectMutli);//1st approximation.
			//std::cout<<lb<<std::endl;/////************
			if(ts>lb.param_s() || te<lb.param_e())
				lb.limit(pspan);

			double te2=lb.param_e();
			if(parameter_normalization){
				double vlen=lb.eval(te2,1).len();
				double tsnew=tnew.param_s();
				tnew-=tsnew;
				tnew*=vlen;
			}
			//std::cout<<lb<<endl;
		}else{//For the span from the 2nd.
			MGLBRep lbt;
			approximate_as_LBRep2(lbt,norder,start-1,index,neglectMutli);//from the 2nd approximation.
			//std::cout<<lbt<<endl;
			if(te<lbt.param_e())
				lbt.limit(pspan);
			int which=2;
			int cn=0;
			if(parameter_normalization){
				double ratio;
				cn=lb.continuity(lbt,which,ratio);
			}
			lb.connect(cn,which,lbt);
		}
		start=index+multi_found;
	}while(index<n);
	if(parameter_normalization==1)
		tnew.change_range(0.,1.);
	lb.update_mark();
}

#define INCREASE_NUM 10
///Get data points for approximate_as_LBRep2.
void MGCurve::data_points_for_approximate_as_LBRep2(
	int is, int ie,//approximation parameter range, from knot_vector()[is] to [ie].
	MGKnotVector& t,//New knot configuration will be output.
				//t's order is input. other information of t will be updated.
	MGNDDArray& tau,//Data point for t will be output.
	bool neglectMulti///<Indicates if multiple knots be kept.
		///< true: multiplicity is removed.
		///< false: multiplicity is kept.
)const{
	const MGKnotVector& told=knot_vector();
	int kold=told.order();
	int knew=t.order();//new order.
	int kdiff=knew-kold;
	int nnew=knew;//nnew will be new B-Rep dimension.
	if(neglectMulti){
		tau.update_from_knot(told);
		tau.change_number(told.bdim()*INCREASE_NUM);
		t=MGKnotVector(tau,knew);
		return;
	}

	int index, multi_found,multi_new;
	int start=is+1;
	do{	//Count the new B-Rep dimension.
		multi_found=told.locate_multi(start,1,index);//Get the next multiplicity.
		multi_new=multi_found+kdiff;
		if(multi_new<1 || multi_found==1)
			multi_new=1;

		nnew+=multi_new+INCREASE_NUM;
		start+=multi_found;
	}while(index<ie);
	nnew-=multi_new;
	t.size_change(knew,nnew);
	
	const double& ts=told(is);
	int i=0;//is the index of t to store the data t(i).
	for(; i<knew; i++)
		t(i)=ts;
	double incNp1=double(INCREASE_NUM+1);

	start=is+1;
	do{	//Count the new bspline rep dimension.
		multi_found=told.locate_multi(start,1,index);//Get the next multiplicity.
		multi_new=multi_found+kdiff;
		if(multi_new<1 || multi_found==1)
			multi_new=1;

		double dif=(told(start)-told(start-1))/incNp1;
		for(int j=0; j<INCREASE_NUM; j++,i++)
			t(i)=t(i-1)+dif;
		for(int j=0; j<multi_new; j++,i++)
			t(i)=told(index);
		start+=multi_found;
	}while(index<ie);
	const double& te=told(ie);
	for(int j=0; j<knew-multi_new; j++)
		t(i++)=te;
	//std::cout<<t<<std::endl;///////********
	tau.update_from_knot(t);assert(tau.length()==nnew);
}

//Approximate this curve as a MGLBRep curve from knot_vector[is] to [ie].
//This is an internal program of MGLBRep constructor.
void MGCurve::approximate_as_LBRep2(
	MGLBRep& lb,		//Approximated LBRep will be set.
	int ordr,		//new order
	int is, int ie,//approximation parameter range, from knot_vector()[is] to [ie].
	bool neglectMulti///<Indicates if multiple knots be kept.
		///< true: multiplicity is removed.
		///< false: multiplicity is kept.
)const{
	//精度十分の曲線を生成する
	MGNDDArray tau;
	MGKnotVector t(ordr);
	data_points_for_approximate_as_LBRep2(is,ie,t,tau,neglectMulti);
	MGBPointSeq bp1;
	eval_line(tau,bp1);

	//制御点を生成する
	lb=MGLBRep(tau, bp1, t);
	lb.remove_knot();
}