/********************************************************************/
/* Copyright (c) 2019 System fugen G.K. and Yuzi Mizuno          */
/* All rights reserved.                                             */
/********************************************************************/
#include "StdAfx.h"
#include "cskernel/Blgi2d.h"
#include "mg/Tolerance.h"
#include "mg/Vector.h"
#include "mg/Unit_vector.h"
#include "mg/Position.h"
#include "mg/Transf.h"
#include "mg/Position_list.h"
#include "mg/Straight.h"
#include "mg/LBRep.h"
#include "mg/Plane.h"
#include "mg/SBRep.h"


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

// MGSBRep5.cpp
//
// Implements Surface B-Representation class MGSBRep.

//Compute ratio at s of the span (s0, s1).
//Ratio computation of linear case.
void get_ratio2(
	double s0,double s1,double s,
	double& r0, double& r1//r0:ratio at s0 side, r1:ration at s1 side will be output.
){
	r0=(s1-s)/(s1-s0);
	r1=1.-r0;
}

double cfunc(double x){
		double x42=x*x;	x42*=4.;
		return x42-x*x42;
}
//Ratio computation of cubic case.
void get_ratio3(
	double s0,double s1,double s,
	double& r0, double& r1//r0:ratio at s0 side, r1:ration at s1 side will be output.
){
	double x=(s1-s)/(s1-s0);
	if(x<=.5){
		r0=cfunc(x); r1=1.-r0;
	}else{
		r1=cfunc(1.-x); r0=1.-r1;
	}	
}

void get_angle(
	const MGLBRep& tie,	//tie curve of rail0 and 1.
	bool& is_straight_tie0,
	bool& is_straight_tie1,
		//boolian that indicates if tie's tangent vector at the start(is_straight_tie0) or
		//end(is_straight_tie1) point is parallel to the vector from the start point
		//to the end point of tie will be output.
	double& theta0,
	double& theta1
		//angles between the surface's normal and the tie's principal normal
		//at the tie's start(theta0) or end(theta1) point will be output.
){
	double s0=tie.param_s(), s1=tie.param_e();
	MGUnit_vector E1=tie.eval(s1)-tie.eval(s0);

//at the start of the tie.
	MGUnit_vector T0,N0,B0;
	double crvtr0,torsn0;
	tie.Frenet_frame(s0,T0,N0,B0,crvtr0,torsn0);
//	std::cout<<T0<<N0<<B0<<",k="<<crvtr0<<",t="<<torsn0<<endl;
	MGUnit_vector E20=T0*E1;//is the normal of the virtual surface.
	is_straight_tie0=false;
	if(MGRZero(crvtr0)){
		if(E1.parallel(T0)) is_straight_tie0=true;
		else B0=T0*E1;
	}
	theta0=E20.angle(B0);

//at the end of the tie.
	MGUnit_vector T1,N1,B1;
	double crvtr1,torsn1;
	tie.Frenet_frame(s1,T1,N1,B1,crvtr1,torsn1);
//	std::cout<<T1<<N1<<B1<<",k="<<crvtr1<<",t="<<torsn1<<endl;
	MGUnit_vector E21=E1*T1;//is the normal of the virtual surface.
	is_straight_tie1=false;
	if(MGRZero(crvtr1)){
		if(E1.parallel(T1)) is_straight_tie0=true;
		else B1=T1*E1;
	}
	theta1=E21.angle(B1);
}

void get_tie(
	double t,			//param of rail0 and 1(the same).
	const MGLBRep& rail0,
	const MGLBRep& rail1,
	const MGLBRep& tie0,	//initial tie(tie at end point).
	const MGPosition& Q0,	//is tie0.start_point().
	const MGPosition& Q1,	//is tie0.end_point().
	bool parallel0, bool parallel1,
		//indicate if tie0's tangent vector at the start(parallel0) or end(parallel1) point
		//is parallel to the vector from Q0 to Q1.
	double angle_s, double angle_e,
		//are angles at the tie0's start(angle_s) or end(angle_e) point between the surface's normal
		//and the tie0's principal normal.
	MGLBRep& tie
){
	MGPosition C=rail0.eval(t), D=rail1.eval(t);
	tie=tie0*MGTransf(Q0,Q1,C,D);
	if(!parallel0 || !parallel1){
		double theta0, theta1;
		bool para_ui0, para_ui1;
		get_angle(tie,para_ui0,para_ui1,theta0,theta1);
		//std::cout<<endl<<"Befor theta="<<theta0<<","<<theta1<<endl;
		double dif0, dif1;
		if(parallel0 || para_ui0) dif0=0.; else dif0=theta0-angle_s;
		if(parallel1 || para_ui1) dif1=0.; else dif1=theta1-angle_e;
		double dif=(dif0+dif1)*.5;//std::cout<<" dif0="<<dif0<<",dif1="<<dif1<<",dif="<<dif<<endl;
		MGTransf tr; tr.set_rotate_3D(D-C,dif,C);
		tie*=tr;
		get_angle(tie,para_ui0,para_ui1,theta0,theta1);
		//std::cout<<"After theta="<<theta0<<","<<theta1<<endl;
	}
}

//Given two curves(peri0 and peri1), get a curve(peri2) that is peri1's slide
//along peri0
//from the peri0's start point to the end point if StoE=true,
//from the peri0's end point to the start point if StoE=false,
//Here, periN=perimeters2[peri[N]].
void get_perixxx(
	int peri[3],
	bool StoE,
	std::unique_ptr<MGLBRep> perimeters2[4]
){
	int n0=peri[0], n1=peri[1], n2=peri[2];
	const MGLBRep& peri0=*(perimeters2[n0]);
	const MGLBRep& peri1=*(perimeters2[n1]);
	perimeters2[n2].reset(peri1.clone());
	MGLBRep& peri2=*(perimeters2[n2]);
	double ts=peri0.param_s(), te=peri0.param_e();
	MGVector V0, V1;
	MGPosition P0, P1;
	double t0 = te, t1 = ts;
	if (StoE) {
		t0 = ts;
		t1 = te;
	}
	V0=peri0.eval(t0, 1);
	P0=peri0.eval(t0);
	V1=peri0.eval(t1, 1);
	P1=peri0.eval(t1);
	V0.set_unit(); V1.set_unit();
	MGMatrix M; M.set_rotate(V0,V1);
	peri2-=P0;
	peri2*=M;
	peri2+=P1;
}

//Given two curves(peri0 and 3), get a curve(peri1) that is peri3's slide along peri0
//from the peri0's start point to the end point.
//peri0's start point and peri3's start point must coincide.
void get_peri031(
	std::unique_ptr<MGLBRep> perimeters2[4]
){
	int peri[3]={0, 3, 1};
	get_perixxx(peri,true, perimeters2);
}
//Given two curves(peri0 and 1), get a curve(peri3) that is peri1's slide along peri0
//from the peri0's end point to the start point.
//peri0's end point and peri1's start point must coincide.
void get_peri013(
	std::unique_ptr<MGLBRep> perimeters2[4]
){
	int peri[3]={0, 1, 3};
	get_perixxx(peri,false, perimeters2);
}
//Given two curves(peri1 and 2), get a curve(peri0) that is peri2's slide along peri1
//from the peri1's end point to the start point.
//peri1's end point and peri2's end point must coincide.
void get_peri120(
	std::unique_ptr<MGLBRep> perimeters2[4]
){
	int peri[3]={1,2,0};
	get_perixxx(peri,false, perimeters2);
}
//Given two curves(peri3 and 2), get a curve(peri0) that is peri2's slide along peri3
//from the peri3's end point to the start point.
//peri3's end point and peri2's start point must coincide.
void get_peri320(
	std::unique_ptr<MGLBRep> perimeters2[4]
){
	int peri[3]={3,2,0};
	get_perixxx(peri,false, perimeters2);
}

//Generalized ruled surface construction.
//Build a surface by ruled surface method. That is, costruct a surface by sliding
//the blending curve(a tie curve of the rails) of perimeters 3 and 1 along
//perimeter 0 and 2(0 and 2 make the rail).
//Or by sliding the blending curve of perimeter 0 and 2
//along perimeter 3 and 1(3 and 1 make the rail).
void MGSBRep::buildGeneralizedRuledSurface(
	const std::unique_ptr<MGLBRep> perimeters[4],
	//«EόXg(vmin,umax,vmax,uminΜ,ΣΤ0,1,2,3Μ)
	//perimeters must be the same knot configuration along u(perimeter 0 and 2)
	//and along v(perimeter 3 and1).
	bool along_u	//indicates which perimeters make a rail.
			//if true, perimeter 0 and 2, if false, perimeter 3 and 1 make rail.
){
	const MGLBRep* peris[4];
	extractConstPointerVec(perimeters, perimeters+4, peris);
	buildGeneralizedRuledSurface(peris, along_u);
}
void MGSBRep::buildGeneralizedRuledSurface(const MGLBRep * perimeters[4], bool along_u){
	invalidateBox();

	const MGLBRep& peri0 = *(perimeters[0]);
	const MGLBRep& peri1 = *(perimeters[1]);
	const MGLBRep& peri2 = *(perimeters[2]);
	const MGLBRep& peri3 = *(perimeters[3]);
	int sd = peri0.sdim();
	int nu = peri0.bdim(), nv = peri1.bdim();
	int num1 = nu-1, nvm1 = nv-1;
	int len = nu; if(len<nv) len = nv;
	int one = 1;

	double u0 = peri0.param_s(), u1 = peri0.param_e();
	double uspan = u1-u0;
	double v0 = peri1.param_s(), v1 = peri1.param_e();
	double vspan = v1-v0;
	const MGKnotVector& uknot = peri0.knot_vector();
	const MGKnotVector& vknot = peri1.knot_vector();
	int orderu = uknot.order(), orderv = vknot.order();
	int order = orderu; if(order<orderv) order = orderv;

	double* work = new double[len+len*(2*order-1)];
	double* q = work+len;

	MGPosition P[4] =//Four corner point.
	{peri3.start_point(), peri1.start_point(),peri1.end_point(), peri3.end_point()};

	double alpha0, alpha1, alpha2, alpha3;
	bool parallel0, parallel1, parallel2, parallel3;
	m_uknot = uknot; m_vknot = vknot;
	m_surface_bcoef.resize(nu, nv, sd);
	double r0, r1;

	if(along_u){
		MGSPointSeq sp(nu, nv, sd);//Temporal spoint.
		int usize = sp.capacity_u();
		MGNDDArray utau; utau.buildByKnotVector(uknot);
		get_angle(peri3, parallel0, parallel3, alpha0, alpha3);
		get_angle(peri1, parallel1, parallel2, alpha1, alpha2);

		sp.store_BP_along_v_at(0, peri3.line_bcoef());
		for(int i = 1; i<num1; i++){
			double ui = utau[i];
			MGLBRep tiei0;
			get_tie(ui, peri0, peri2, peri3, P[0], P[3], parallel0, parallel3, alpha0, alpha3, tiei0);
			MGLBRep tiei1;
			get_tie(ui, peri0, peri2, peri1, P[1], P[2], parallel1, parallel2, alpha1, alpha2, tiei1);
			get_ratio2(u0, u1, ui, r0, r1);//***********
			sp.store_BP_along_v_at(i, tiei0.line_bcoef()*r0+tiei1.line_bcoef()*r1);
		}
		sp.store_BP_along_v_at(num1, peri1.line_bcoef());
		int error = 2;
		for(int k = 0; k<sd; k++){
			for(int j = 0; j<nv; j++)
				blgi2d_(&error, utau.data(), sp.data(0, j, k), m_uknot.data(),
					orderu, nu, one, usize, one, work, q, &m_surface_bcoef(0, j, k));
			if(error!=1) break;
		}
	} else{
		MGSPointSeq sp(nv, nu, sd);//Temporal spoint.
		int vsize = sp.capacity_u();
		MGNDDArray vtau; vtau.buildByKnotVector(vknot);
		get_angle(peri0, parallel0, parallel1, alpha0, alpha1);
		get_angle(peri2, parallel3, parallel2, alpha3, alpha2);

		sp.store_BP_along_v_at(0, peri0.line_bcoef());
		for(int j = 1; j<nvm1; j++){
			double vj = vtau[j];
			MGLBRep tiej0;
			get_tie(vj, peri3, peri1, peri0, P[0], P[1], parallel0, parallel1, alpha0, alpha1, tiej0);
			MGLBRep tiej1;
			get_tie(vj, peri3, peri1, peri2, P[3], P[2], parallel3, parallel2, alpha3, alpha2, tiej1);
			get_ratio2(v0, v1, vj, r0, r1);//********
			sp.store_BP_along_v_at(j, tiej0.line_bcoef()*r0+tiej1.line_bcoef()*r1);
		}
		sp.store_BP_along_v_at(nvm1, peri2.line_bcoef());
		int error = 2;
		for(int k = 0; k<sd; k++){
			blgi2d_(&error, vtau.data(), sp.data(0, 0, k), m_vknot.data(),
				orderv, nv, nu, vsize, nu, work, q, &m_surface_bcoef(0, 0, k));
			if(error!=1) break;
		}
	}
	delete[]  work;
	copy_appearance(*(perimeters[0]));
}

//Construct 4 perimeters, given at least two of the four.
//Input perimeters may have different knot configuration.
//In this case they will be updated so as to have the same configuration.
//Function's return value indicates which perimeter(s) was missing:
// 10: all of the 4 were input(and knot configurations were updated to have the same).
//  0: only perimeter 0 was missing.
//  1: only perimeter 1 was missing.
//  2: only perimeter 2 was missing.
//  3: only perimeter 3 was missing.
//  4: perimeter 2 and 3 were missing.
//  5: perimeter 1 and 3 were missing.
//  6: perimeter 1 and 2 were missing.
//  7: perimeter 0 and 3 were missing.
//  8: perimeter 0 and 2 were missing.
//  9: perimeter 0 and 1 were missing.
// -2: less than 2 perimeters were provided.
//The missing perimeters are made by straight lines in the cases of 5 or 8.
int construct_perimeters(
	const MGCurve* perimeters[4],
		//«EόXg(vmin,umax,vmax,uminΜ,ΣΤ0,1,2,3Μ). Let i be the perimeter number,
		//and the data is missing, perimeter[i] must be null. If perimeter 3 data is missing,
		//perimeters[3] is nullptr. If perimeter 2 and 3 data are missing, 
		//perimeters[2]=perimeters[3]=nullptr.
		//When perimeters were not the same knot configuration along u(perimeter 0 and 2)
		//or along v(perimeter 3 and1), they will be rebuild to have the same configuration.
		//At least 2 perimeters must be provided.
	std::unique_ptr<MGLBRep> perimeters2[4]
			//new perimeters will be output.
){
//1. Classify which perimeters are input.
	int type;
	if(!perimeters[0]){//0 missing
		if(!perimeters[3]){//3 missing
			if(!perimeters[1] || !perimeters[2]) return -2;
			type=7;
		}else{//3 provided
			if(perimeters[1]){//1 provided
				if(perimeters[2]) type=0;
				else              type=8;
			}else{
				if(!perimeters[2]) return -2;
				type=9;
			}
		}
	}else{//0 provided
		if(perimeters[1]){//1 provided
			if(perimeters[2]){//2 provided
				if(perimeters[3]) type=10;
				else type=3;
			}else{
				if(perimeters[3]) type=2;
				else type=4;
			}
		}else{//1 missing
			if(perimeters[2]){//2 provided
				if(perimeters[3]) type=1;
				else type=5;
			}else{
				if(perimeters[3]) type=6;
				else return -2;
			}
		}
	}

//2. knot configuration adjustment.
	std::vector<const MGCurve*> peri2(2);
	if(type==10 || type==1 || type==3 || type==5){//when 0 and 2 provided.
		peri2[0]=perimeters[0];
		peri2[1]=perimeters[2];
		std::vector<UniqueLBRep> lb2=rebuildAsSameKnotVector(peri2);
		perimeters2[0].reset(lb2[0].release());
		perimeters2[2].reset(lb2[1].release());
	}
	if(type==10 || type==0 || type==2 || type==8){//when 1 and 3 provided.
		peri2[0]=perimeters[1];
		peri2[1]=perimeters[3];
		std::vector<UniqueLBRep> lb2=rebuildAsSameKnotVector(peri2);
		perimeters2[1].reset(lb2[0].release());
		perimeters2[3].reset(lb2[1].release());
	}
	if(type==10) return 10;

//3. process when parallel two perimeters are provided and the other two are missing.
	if(type==5){
		perimeters2[1].reset(new MGLBRep(
			MGStraight(perimeters2[2]->end_point(),perimeters2[0]->end_point())));
		perimeters2[3].reset(new MGLBRep(
			MGStraight(perimeters2[2]->start_point(),perimeters2[0]->start_point())));
		double te=(perimeters2[1]->param_e()+perimeters2[3]->param_e())*.5;
		perimeters2[1]->change_range(0., te);
		perimeters2[3]->change_range(0., te);
		return 5;
	}
	if(type==8){
		perimeters2[0].reset(new MGLBRep(
			MGStraight(perimeters2[1]->start_point(),perimeters2[3]->start_point())));
		perimeters2[2].reset(new MGLBRep(
			MGStraight(perimeters2[1]->end_point(),perimeters2[3]->end_point())));
		double te=(perimeters2[0]->param_e()+perimeters2[2]->param_e())*.5;
		perimeters2[0]->change_range(0., te);
		perimeters2[2]->change_range(0., te);
		return 8;
	}

//4. approximate the perimeter that has not the opposite one.
	for(int i=0; i<4; i++){
		if(!perimeters2[i].get()){
			if(perimeters[i])
				perimeters2[i].reset(new MGLBRep(*(perimeters[i])));
		}
	}

//5. Adjust the common corner points.
	MGLBRep& l0=*(perimeters2[0]);
	MGLBRep& l1=*(perimeters2[1]);
	MGLBRep& l2=*(perimeters2[2]);
	MGLBRep& l3=*(perimeters2[3]);

	MGPosition P0, P1, P2, P3;
	if(type==1 || type==2 || type==6) P0=(l0.start_point()+l3.start_point())*.5;
	if(type==2 || type==3 || type==4) P1=(l0.end_point()+l1.start_point())*.5;
	if(type==0 || type==3 || type==7) P2=(l1.end_point()+l2.end_point())*.5;
	if(type==0 || type==1 || type==9) P3=(l3.end_point()+l2.start_point())*.5;

	double fixp[2];
	if(!P0.is_null()){//P0
		fixp[0]=l0.param_e(); l0.move(2,l0.param_s(),P0,fixp);
		fixp[0]=l3.param_e(); l3.move(2,l3.param_s(),P0,fixp);
	}
	if(!P1.is_null()){//P1
		fixp[0]=l0.param_s(); l0.move(2,l0.param_e(),P1,fixp);
		fixp[0]=l1.param_e(); l1.move(2,l1.param_s(),P1,fixp);
	}
	if(!P2.is_null()){//P2
		fixp[0]=l1.param_s(); l1.move(2,l1.param_e(),P2,fixp);
		fixp[0]=l2.param_s(); l2.move(2,l2.param_e(),P2,fixp);
	}
	if(!P3.is_null()){//P3
		fixp[0]=l2.param_e(); l2.move(2,l2.param_s(),P3,fixp);
		fixp[0]=l3.param_s(); l3.move(2,l3.param_e(),P3,fixp);
	}

//6. process when non-parallel two perimeters are provided.
	int type2=type;
	if(type==4){
		get_peri013(perimeters2);type2=2;
	}else if(type==6){
		get_peri031(perimeters2);type2=2;
	}else if(type==7){
		get_peri120(perimeters2);type2=3;
	}else if(type==9){
		get_peri320(perimeters2);type2=1;
	}
	//std::cout<<perimeters2;

//7. Now 3 of the 4 perimeters are constructed, construct the missing one.
	bool para0,para1;
	double theta0, theta1;
	perimeters2[type2].reset(new MGLBRep);
	MGLBRep& peri=*(perimeters2[type2]);//peri to be constructed.

	MGLBRep& peri_opo=*(perimeters2[(type2+2)%4]);//opposide perimeter of peri.
	int is;
	if(type2<=1) is=(type2+3)%4; else is=(type2+1)%4;
	MGLBRep& peri_s=*(perimeters2[is]);//starting side perimeter of peri.
	MGLBRep& peri_e=*(perimeters2[(is+2)%4]);//ending side perimeter of peri.

	get_angle(peri_opo,para0,para1,theta0,theta1);
	MGPosition Q0=peri_opo.start_point(), Q1=peri_opo.end_point();
	double t;
	if(type2==0 || type2==3) t=peri_s.param_s(); else t=peri_s.param_e();
	get_tie(t,peri_s,peri_e,peri_opo,Q0,Q1,para0,para1,theta0,theta1,peri);
	return type;
}
