论文：Linear Angle Based Parameterization(2007)的代码片段

class LinABF::Private::Implement :
	delegate public Abstract {
private:
	using HINT_PROC = DEF<void (const INDEX &)> ;
	typedef  Eigen::SparseMatrix<double,Eigen::RowMajor> SpMatR ;
	typedef  Eigen::SparseMatrix<double> SpMat ;
	const Array<Vector<VALXA>>& mMeshVertices ;
	const Array<ARRAY3<INDEX>>& mMeshFaces ; 
	BitSet<> mBoundary ;


public:
	implicit Implement () = default ;

    explicit Implement (const Array<Vector<VALXA>>& mMeshVertex ,const Array<ARRAY3<INDEX>>& mMeshFace) 
		: mMeshVertices(mMeshVertex) ,mMeshFaces (mMeshFace){}


	void solve (Array<ARRAY2<VALXA>>& mMeshVertexUV )  {
		
		Eigen::MatrixXd solution = normal_equation_setup () ;
		// mMeshVertexUV = least_square_calculate_texture_uv (solution) ;
		// mMeshVertexUV = calculate_texture_uv (solution) ;
		mMeshVertexUV = lscm_compute (solution) ;

		
		
	}

private:
	BitSet<> find_boundary () {
		BitSet<> ret = BitSet<> (mMeshVertices.length()) ;

		Set<ARRAY2<INDEX>> mEdage = Set<ARRAY2<INDEX>> (mMeshFaces.length() * 3) ;
		ArrayList<INDEX> mEdage_info = ArrayList<INDEX> (mEdage.size()) ;

		for (auto&& i : _RANGE_(0 ,mMeshFaces.length())) {
			const auto& i11x = mMeshFaces[i] ;
			for (auto &&j : _RANGE_(0 ,3)) {
				const auto& j11x = (j+1) % 3 ;
				ARRAY2<INDEX> j12x = ARRAY2<INDEX> {VAR_NONE ,VAR_NONE} ;

				if (i11x[j] < i11x[j11x] ) {
					j12x[0] = i11x[j] ;
					j12x[1] = i11x[j11x] ;
				}
				else {
					j12x[1] = i11x[j] ;
					j12x[0] = i11x[j11x] ;
				}

				INDEX ix = mEdage.map (j12x) ;
				if (ix == VAR_NONE) {
					ix = mEdage_info.insert () ;
					mEdage_info[ix] = 0 ;
					mEdage.add (j12x ,ix) ;
				}
				mEdage_info[ix]++ ;
			}
		}

		for (auto &&i : mEdage) {
			const auto i11x = mEdage.at(i) ;
			const auto i12x = mEdage.map_get(i11x) ;
			_DYNAMIC_ASSERT_(i12x != VAR_NONE) ;
			if (mEdage_info[i12x] >= 2) 
				continue ;
			
			ret[i[0]] = TRUE ;
			ret[i[1]] = TRUE ;
		}
		return _MOVE_(ret) ;
	}

	Array<VALXA> calculate_triangle_angle () {
		Array<VALXA> ret = Array<VALXA> (mMeshFaces.length() * 3) ;
		for (auto &&i : _RANGE_(0 , mMeshFaces.length())) {
			const auto i10x = i * 3 ;
			const auto& i11x = mMeshFaces[i] ;

			const auto i12x = (mMeshVertices[i11x[1]] - mMeshVertices[i11x[0]]).normalize () ;
			const auto i13x = (mMeshVertices[i11x[2]] - mMeshVertices[i11x[1]]).normalize () ;
			const auto i14x = (mMeshVertices[i11x[0]] - mMeshVertices[i11x[2]]).normalize () ;
			
			ret[i10x + 0] = MathProc::arccos(i12x * i14x * VALXA(-1.0)) ;
			ret[i10x + 1] = MathProc::arccos(i13x * i12x * VALXA(-1.0)) ;
			ret[i10x + 2] = MathProc::arccos(i14x * i13x * VALXA(-1.0)) ;
		}
		return _MOVE_ (ret) ;
	}

	Array<INDEX> internal_vertex_index_sorting (const BitSet<>& boundary) {
		Array<INDEX> ret = Array<INDEX> (mMeshVertices.length()) ;
		INDEX ix = 0 ;
		for (auto&& i : _RANGE_(0 , mMeshVertices.length())) {
			ret[i] = VAR_NONE ;
			if (boundary[i])
				continue ;
			ret[i] = ix ;
			ix++ ;
		}
		return _MOVE_ (ret) ;
	}

	/* 在论文 4.2. Choice of initial estimation (第5页) */
	Array<VALXA> choice_of_initial_estimation (const BitSet<>& boundary) {
		auto ang123 = calculate_triangle_angle () ;

		const auto vlen = mMeshVertices.length () ;

		std::vector<Eigen::Triplet<double, int>> Mang_T (mMeshFaces.length() * 3) ;
		for (auto &&i : _RANGE_(0 , mMeshFaces.length())) {
			const auto i11x = i * 3 ;
			const auto& i12x = mMeshFaces[i] ;
			Eigen::Triplet<double, int> t1(i12x[0], i12x[1], ang123[i11x + 0]) ;
			Eigen::Triplet<double, int> t2(i12x[1], i12x[2], ang123[i11x + 1]) ;
			Eigen::Triplet<double, int> t3(i12x[2], i12x[0], ang123[i11x + 2]) ;
			Mang_T.push_back (t1) ; 
			Mang_T.push_back (t2) ; 
			Mang_T.push_back (t3) ; 
		}

		SpMatR Mang = Eigen::SparseMatrix<double,Eigen::RowMajor>(vlen,vlen) ;
		Mang.setFromTriplets (Mang_T.begin () , Mang_T.end()) ;

		// 计算环角
		Eigen::VectorXd sumringangles(vlen);
		for(int i = 0; i < vlen; i++) {
			sumringangles(i) = Mang.row(i).sum();
		}
		
		// 优化初始角度值
		const auto _2PI = VALXA (2) * M_PI ;
		Array<VALXA> ang123opt = Array<VALXA> (ang123.length()) ;
		for (auto &&i : _RANGE_(0 , mMeshFaces.length())) {
			const auto i11x = i * 3 ;
			const auto& i12x = mMeshFaces[i] ;

			
			const auto i14x = _2PI * ang123[i11x + 0] * MathProc::inverse (sumringangles[i12x[0]]) ;
			const auto i15x = _2PI * ang123[i11x + 1] * MathProc::inverse (sumringangles[i12x[1]]) ;
			const auto i16x = _2PI * ang123[i11x + 2] * MathProc::inverse (sumringangles[i12x[2]]) ;
			ang123opt[i11x + 0] = i14x ;
			ang123opt[i11x + 1] = i15x ;
			ang123opt[i11x + 2] = i16x ;
		}

		// 计算亏损的角度 
		Eigen::VectorXd rhsk = Eigen::VectorXd::Constant(vlen, _2PI) - sumringangles ;
		BitSet<> icha = BitSet<> (vlen) ;
		// 对于角度缺陷大于 1 的顶点环，建议切换到从该特定环的指数映射中获得的角度。
		for (auto &&i : _RANGE_(0 , vlen)) {
			if (MathProc::abs(rhsk[i]) > 1.0) {
				icha[i] = TRUE ;
			}
		}
		
		icha = icha - boundary ;

		for (auto &&i : _RANGE_(0 , mMeshFaces.length())) {
			const auto i11x = i * 3 ;
			const auto& i12x = mMeshFaces[i] ;
			for (auto &&j : _RANGE_(0 ,3)) {
				if (icha[i12x[j]]) {
					ang123[i11x + j] = ang123opt[i11x + j] ;
				}
			}
		}

		return _MOVE_(ang123) ;
	}


	/**	在论文 3. Reformulation and linearization (第3页右边)
	 *  
	 * 	Vertex consistency　(等式５)
	*/
	void vertex_consistency (const Array<VALXA>& mAngles ,
							 const BitSet<>& boundary ,
							 const Array<INDEX>& mInternal_index ,
							 SpMat& Ver_consistency,
							 Eigen::VectorXd& Ver_consistency_rhs 
							) 
	{
		
		const auto r1x = mMeshFaces.length() ;
		std::vector<Eigen::Triplet<double, int>> coefficients (r1x * 3);

		for (auto &&i : _RANGE_(0 ,r1x)) {
			const auto i11x = i * 3 ;
			const auto& i12x = mMeshFaces[i] ;
			for (auto &&j : _RANGE_(0 ,3)) {
				if (boundary[i12x[j]])
					continue ;
				INDEX tem_row = mInternal_index[i12x[j]] ;
				_DYNAMIC_ASSERT_(tem_row != VAR_NONE) ;
				Eigen::Triplet<double, int> t1 (tem_row ,(i11x + j) ,1.0) ;
				coefficients.push_back(t1) ;
			}
		}
		
		int Ver_rows = mMeshVertices.length() - boundary.length() ;
		int Ver_cols = mAngles.length() ;

		Ver_consistency = Eigen::SparseMatrix<double> (Ver_rows, Ver_cols) ;
		Ver_consistency.setFromTriplets( coefficients.begin(), coefficients.end()) ;

		
		Eigen::VectorXd mtemPA = Eigen::VectorXd::Zero(mAngles.length());
		for (auto &&i : _RANGE_(0 ,mAngles.length())) {
			mtemPA[i] = mAngles[i] ;
		}

		Eigen::VectorXd tempAangle = Ver_consistency * mtemPA ;
		Eigen::VectorXd two_pi = Eigen::VectorXd::Constant(tempAangle.rows(), 2 * M_PI) ;
		Ver_consistency_rhs = two_pi - tempAangle ;

	}
	

	/**	在论文 3. Reformulation and linearization (第3页右边)
	 *  
	 * 	Triangle consistency (等式６)
	*/
	void triangle_consistency( const Array<VALXA>& mAngles ,
							   SpMat& Tri_consistency ,
							   Eigen::VectorXd& Tri_consistency_rhs 
							 ) 
	{

		const auto r1x = mMeshFaces.length() ;
		Tri_consistency_rhs = Eigen::VectorXd (r1x) ;
		std::vector<Eigen::Triplet<double, int>> coefficients (r1x * 3);
		for (auto &&i : _RANGE_(0 , r1x)) {
			const auto i11x = i * 3 ;
			Eigen::Triplet<double, int> t1 (i ,(i11x + 0) ,1.0) ;
			Eigen::Triplet<double, int> t2 (i ,(i11x + 1) ,1.0) ;
			Eigen::Triplet<double, int> t3 (i ,(i11x + 2) ,1.0) ;
			coefficients.push_back (t1) ;
			coefficients.push_back (t2) ;
			coefficients.push_back (t3) ;
			Tri_consistency_rhs[i] = VALXA(M_PI) - (mAngles[i11x + 0] + mAngles[i11x + 1] + mAngles[i11x + 2]) ;
		}
		
		Tri_consistency = Eigen::SparseMatrix<double> (r1x,mAngles.length()) ;
		Tri_consistency.setFromTriplets (coefficients.begin () , coefficients.end()) ;

	}

	
	/**	在论文 3. Reformulation and linearization (第3页右边)
	 *  
	 * 	Wheel consistency　(等式１０)
	*/
	void wheel_consistency (const Array<VALXA>& mAngles ,
							const BitSet<>& boundary ,
							const Array<INDEX>& mInternal_index ,
							SpMat& Wheel_consistency_L ,
							Eigen::VectorXd& Wheel_consistency_rhs 
							)
	{
		const auto r1x = mMeshFaces.length() ;
		const auto r2x = mMeshVertices.length() ;

		// 左边等式
		std::vector<Eigen::Triplet<double, int>> coefficients (r1x * 3) ;
		for (auto &&i  : _RANGE_(0 ,r1x)) {
			const auto i11x = i * 3 ;
			const auto& i12x = mMeshFaces[i] ;

			for (auto &&j : _RANGE_(0 , 3)) {
				if (boundary[i12x[j]])
					continue ;
				const auto j15x = (j + 1) % 3 ;
				const auto j16x = (j + 2) % 3 ;

				INDEX mTemp_row = mInternal_index[i12x[j]] ;
				_DYNAMIC_ASSERT_ (mTemp_row != VAR_NONE) ;
	
				const auto j20x = MathProc::tan(mAngles[i11x + j15x]) ;
				const auto j21x = VALXA(1.0) * MathProc::inverse (j20x) ;
				Eigen::Triplet<double, int> t1 (mTemp_row ,(i11x + j15x) ,j21x) ;
				coefficients.push_back (t1) ;
				
				const auto j22x = MathProc::tan(mAngles[i11x + j16x]) ;
				const auto j23x = VALXA(-1.0) * MathProc::inverse (j22x) ;
				Eigen::Triplet<double, int> t2 (mTemp_row ,(i11x + j16x) ,j23x) ;
				coefficients.push_back (t2) ;
			}
		}
		
		const auto Wh_rows = mMeshVertices.length() - boundary.length() ;
		const auto Wh_cols = mAngles.length() ;

		Wheel_consistency_L = Eigen::SparseMatrix<double> (Wh_rows, Wh_cols) ;
		Wheel_consistency_L.setFromTriplets( coefficients.begin(), coefficients.end()) ;
		coefficients.clear () ;

		for (auto &&i  : _RANGE_(0 ,r1x)) {
			const auto i11x = i * 3 ;
			const auto& i12x = mMeshFaces[i] ;

			for (auto &&j : _RANGE_(0 , 3)) {
				if (boundary[i12x[j]])
					continue ;
				const auto j15x = (j + 1) % 3 ;
				const auto j16x = (j + 2) % 3 ;
	
				const auto j20x = MathProc::sin(mAngles[i11x + j16x]) ;
				const auto j21x = VALXA(1.0) * MathProc::log (j20x) ;
				Eigen::Triplet<double, int> t1 (i12x[j] ,(i11x + j16x) ,j21x) ;
				coefficients.push_back (t1) ;
				
				const auto j22x = MathProc::sin(mAngles[i11x + j15x]) ;
				const auto j23x = VALXA(-1.0) * MathProc::log (j22x) ;
				Eigen::Triplet<double, int> t2 (i12x[j] ,(i11x + j15x) ,j23x) ;
				coefficients.push_back (t2) ;
			}
		}

		SpMatR Wheel_consistency_R = SpMatR (r2x, Wh_cols) ;
		Wheel_consistency_R.setFromTriplets( coefficients.begin(), coefficients.end()) ;
		
		Eigen::VectorXd sumCC(r2x) ;
		for (auto &&i : _RANGE_(0 ,r2x)) {
			sumCC(i) = Wheel_consistency_R.row(i).sum() ;
		}

		Wheel_consistency_rhs = Eigen::VectorXd::Zero (Wh_rows) ;
		for (auto &&i : _RANGE_(0 ,r2x)) {
			if (boundary[i])
				continue ;
			const auto& i11x = mInternal_index[i] ;
			_DYNAMIC_ASSERT_(i11x != VAR_NONE) ;
			Wheel_consistency_rhs[i11x] = sumCC(i); 
		}

	}


	/** 4.1. Normal equation setup (第４页)
	 * loss function = sum(1 / ai^2 · (ai* - ai)^2)
	 * ri = ei / ai
	 * 
	 * minimize \\r\\^2 
	 * std:Cr = b
	 * 
	 * C = A·Da
	 * Da = diag(ai)
	 * 
	 * 
	 * r = Ct · x
	 * r = Da · r
	 * a* = a + r
	*/
	Eigen::MatrixXd normal_equation_setup () {
		auto boundary = find_boundary () ;
		mBoundary = boundary ;
		auto mInternal_index = internal_vertex_index_sorting (boundary) ;

		// Choice of initial estimation
		auto mAngles = choice_of_initial_estimation (boundary) ;

		// Vertex consistency
		SpMat Ver_consistency ;
		Eigen::VectorXd Ver_consistency_rhs ;
		vertex_consistency (mAngles ,boundary ,mInternal_index ,Ver_consistency ,Ver_consistency_rhs) ;

		// Triangle consistency
		SpMat Tri_consistency ;
		Eigen::VectorXd Tri_consistency_rhs ;
		triangle_consistency (mAngles ,Tri_consistency ,Tri_consistency_rhs) ;

		// Wheel consistency
		SpMat Wheel_consistency_L ;
		Eigen::VectorXd Wheel_consistency_rhs ;
		wheel_consistency (mAngles ,boundary ,mInternal_index ,Wheel_consistency_L ,Wheel_consistency_rhs) ;

		int Tri_rows = Tri_consistency.rows () ;
		int Ver_rows = Ver_consistency.rows () ;
		int Wh_rows = Wheel_consistency_L.rows () ;

		// A 矩阵的维度(nt + 2 · ni )　x (3 · nt)
		int A_rows = Tri_rows + Ver_rows + Wh_rows ;
		int A_cols = mAngles.length () ;
		SpMatR A = SpMatR (A_rows ,A_cols) ;
		A.middleRows(0, Tri_rows) = Tri_consistency ;
		A.middleRows(Tri_rows, Ver_rows) = Ver_consistency ;
		A.middleRows(Tri_rows + Ver_rows, Wh_rows) = Wheel_consistency_L ;

		// Dα = diag (αi) 
		std::vector<Eigen::Triplet<double, int>> coefficients (A_cols) ;
		for (auto &&i : _RANGE_ (0 , A_cols)) {
			Eigen::Triplet<double, int> t1 (i ,i ,mAngles[i]) ;
			coefficients.push_back (t1) ;
		}
		SpMat Da = SpMat (A_cols ,A_cols) ;
		Da.setFromTriplets(coefficients.begin() ,coefficients.end()) ; 
		coefficients.clear() ;

		// matrix C = A * Dα
		SpMat C = A * Da ;
		SpMat Ct = C.transpose () ;


		// b 
		INDEX rhs1 = INDEX (Tri_consistency_rhs.size()) ;
		INDEX rhs2 = INDEX (Ver_consistency_rhs.size()) ;
		INDEX rhs3 = INDEX (Wheel_consistency_rhs.size()) ;

		Eigen::VectorXd rhs = Eigen::VectorXd::Zero (rhs1 + rhs2 + rhs3) ;
		INDEX iz = 0 ;
		for (auto &&i : _RANGE_(0 ,rhs1)) {
			rhs[iz] = Tri_consistency_rhs[i] ;
			iz++ ;
		}
		for (auto &&i : _RANGE_(0 , rhs2)) {
			rhs[iz] = Ver_consistency_rhs[i] ;
			iz++ ;
		}
		for (auto &&i : _RANGE_(0 , rhs3)) {
			rhs[iz] = Wheel_consistency_rhs[i] ;
			iz++ ;
		}
		
		// r = e / a (等式13)
		// Solve:
		// Eigen::CholmodSupernodalLLT< Eigen::SparseMatrix<double>> solver (C * Ct) ;
		// Eigen::VectorXd r = Ct * solver.solve (rhs) ;
		Eigen::SimplicialLLT<Eigen::SparseMatrix<double>> solver(C * Ct) ;
		Eigen::VectorXd r = Ct * solver.solve (rhs) ;
		r = Da * r ;
		if(solver.info() != Eigen::Success) {
			
			Singleton<ConsoleService>::instance ().error (_PCSTR_ ("normal_equation_setup fail")) ;
			Singleton<ConsoleService>::instance ().error (_PCSTR_ ("info ") ,solver.info()) ;
		}
		Eigen::MatrixXd solution (mMeshFaces.length() , 3) ;
		for (auto &&i : _RANGE_(0 ,mMeshFaces.length())) {
			for (auto &&j : _RANGE_(0 ,3)) {
				const auto j11x = i * 3 + j ;
				solution (i ,j) = mAngles[j11x] + r[j11x] ;
			}
		}

		return solution ;
	}


	/**  最小二乘重建算法
	 * 	设一个三角形T(p1,p2,p3)
	 *  1. ||p1p3|| / ||p1p2|| = sin(a1) / sin(a1)
	 *  得
	 * 	p1p3 = sin(a1) / sin(a2) * {{cos(a0) ,sin(a0)} ,{-sin(a0) ,cos(a0)}}
	 * 
	 * 	因此给顶一个三角形的两个顶点位置和３个角度，第３个顶点的位置就可以有上面的公式计算出来
	 * 
	 * 第二步： 使用贪狼算法先确定两个顶点的位置．其他顶点位置就根据上面的公式依次计算
	 * 
	 * 第三步：	重构上面的公式
	 * 
	 *			Mt (pj - pk) + pi - pj = 0 
	 * 				 sin(a1)  { cos(a0)    sin(a0) }
	 * 			Mt = -------  {				       }
	 *				 cos(a2)  {-sin(a0)    cos(a0) }
	 * 
	 *　第四步： 上述方程构成一个线性系统,未知数为顶点在平面上的x,y坐标,线性系统的系数有三角形的角度计算出来
	 * 
	 *  第五步: 因为三角形的角度确定后,只具有缩放和旋转的自由度.因此需要确定两个顶点的坐标来去除线性系统里面的自由度,从而得到唯一的解
	 * 
	 *  第六步:　假设p(nv-1) 和 p(nv) 到（０，０）和（１，0），那么去除变量后，线性系统公式方程数量是2nf，而未知数变量的数量是2(nv - 2).
	 * 
	 *  第七步：　把第三步的矩阵写成：　AP = 0 ;其中矩阵Ａ的大是2nf x 2(nv - 2) ,P是又顶点的x,y坐标构成
	 * 
	*/
	Array<ARRAY2<VALXA>> calculate_texture_uv (const Eigen::MatrixXd& solution) {
		
		const auto r1x = mMeshFaces.length() ;
		const auto r2x = mMeshVertices.length() ;

		Eigen::VectorXd coef (r1x) ;
		Eigen::VectorXd a (r1x) ;   // G * cos
		Eigen::VectorXd b (r1x) ;	// G * sin
		for (auto &&i : _RANGE_(0 ,r1x)) {
			
			const auto i11x = MathProc::sin (solution(i,1)) ;
			const auto i12x = MathProc::sin (solution(i,2)) ;
			coef[i] = i11x * MathProc::inverse (i12x) ;

			a[i] = coef[i] * MathProc::cos(solution(i,0)) ;
			b[i] = coef[i] * MathProc::sin(solution(i,0)) ;

		}
		
		std::vector<Eigen::Triplet<double, int>> MV_T ;
		for (auto &&i : _RANGE_(0 , r1x)) {
			const auto& i11x = mMeshFaces[i][0] ;
			const auto& i12x = mMeshFaces[i][1] ;
			const auto& i13x = mMeshFaces[i][2] ;

			if switch_once (TRUE) {
				Eigen::Triplet<double, int> t1 (i ,i11x ,1.0 - a(i)) ;
				Eigen::Triplet<double, int> t2 (i ,i11x + r2x ,b(i)) ;

				Eigen::Triplet<double, int> t3 (i ,i12x ,a(i)) ;
				Eigen::Triplet<double, int> t4 (i ,i12x + r2x ,-b(i)) ;

				Eigen::Triplet<double, int> t5 (i ,i13x ,-1.0) ;

				MV_T.push_back (t1) ;
				MV_T.push_back (t2) ;
				MV_T.push_back (t3) ;
				MV_T.push_back (t4) ;
				MV_T.push_back (t5) ;
			}
			
			INDEX ix = i + r1x ;
			if switch_once (TRUE) {
				Eigen::Triplet<double, int> t1 (ix ,i11x ,-b(i)) ;
				Eigen::Triplet<double, int> t2 (ix, i11x + r2x	,1.0 - a(i)) ;

				Eigen::Triplet<double, int> t3 (ix ,i12x ,b(i)) ;
				Eigen::Triplet<double, int> t4 (ix ,i12x + r2x ,a(i)) ;

				Eigen::Triplet<double, int> t5 (ix ,i13x + r2x ,-1.0) ;

				MV_T.push_back (t1) ;
				MV_T.push_back (t2) ;
				MV_T.push_back (t3) ;
				MV_T.push_back (t4) ;
				MV_T.push_back (t5) ;
			}

		}

		INDEX MV_rows = 2 * r1x ;
		INDEX MV_cols = 2 * r2x ;

		// (xxxxxxx yyyyyyy)
		SpMatR MV = Eigen::SparseMatrix<double,Eigen::RowMajor> (MV_rows, MV_cols) ;
		MV.setFromTriplets (MV_T.begin(), MV_T.end()) ;
		SpMatR MVt = MV.transpose () ;
		MV = MVt * MV ;

		// 固定顶点条件
		INDEX a1 = mMeshFaces[0][0] ;
		INDEX a2 = mMeshFaces[0][1] ;
		ARRAY4<INDEX> pinned = ARRAY4<INDEX> { a1, a2, a1 + r2x, a2 + r2x } ;
		

		SpMatR MV_pinned (4 , 2 * r2x) ;
		for (auto &&i : _RANGE_(0 ,4)) {
			MV_pinned.insert (i ,pinned[i]) = 1 ;
		}
		for (auto &&i : _RANGE_(0 ,4)) {
			MV.middleRows (pinned[i] ,1) = MV_pinned.row(i) ;
		}

		Eigen::VectorXd uvc = Eigen::VectorXd::Zero (2 * r2x) ;
		uvc[pinned[1]] = 1.0 ;
	

		MVt = MV.transpose();
		// Eigen::CholmodSupernodalLLT< Eigen::SparseMatrix<Scalar> > postion_solver( MVt * MV );
		// Eigen::VectorXd uv = postion_solver.solve( MVt * uvc );
		Eigen::SimplicialLLT<Eigen::SparseMatrix<double>> solver (MVt * MV) ;
		Eigen::VectorXd uv = solver.solve( MVt * uvc );
		if(solver.info() != Eigen::Success) {
			
			Singleton<ConsoleService>::instance ().error (_PCSTR_ ("calculate_texture_uv fail")) ;
			Singleton<ConsoleService>::instance ().error (_PCSTR_ ("info ") ,solver.info()) ;
		}
		// 保存最终UV坐标
		Array<ARRAY2<VALXA>> tex_coord = Array<ARRAY2<VALXA>> (r2x) ;
		VALXA dx = 0 ;
		VALXA minx = VALXA_INF ;
		VALXA maxx = -VALXA_INF ;

		VALXA dy = 0 ;
		VALXA miny = VALXA_INF ;
		VALXA maxy = -VALXA_INF ;

		for (auto &&i : _RANGE_(0, r2x)) {
			tex_coord[i] = ARRAY2<VALXA> {uv[i] ,uv(i + r2x)} ;
			minx = MathProc::minof (minx ,tex_coord[i][0]) ;
			maxx = MathProc::maxof (maxx ,tex_coord[i][0]) ;
			miny = MathProc::minof (miny ,tex_coord[i][1]) ;
			maxy = MathProc::maxof (maxy ,tex_coord[i][1]) ;
		}
		dx = maxx - minx ;
		dy = maxy - miny ;

		for (auto &&i : _RANGE_(0, tex_coord.length())) {
			tex_coord[i][0] = (tex_coord[i][0] - minx) * MathProc::inverse (dx) ;
			tex_coord[i][1] = (tex_coord[i][1] - miny) * MathProc::inverse (dy) ;

		}
		return _MOVE_ (tex_coord) ;
	}
	
	Array<ARRAY2<VALXA>> lscm_compute (const Eigen::MatrixXd& solution) {

		using Triplet = Eigen::Triplet<double, int> ;
		using Solver = Eigen::SparseLU<Eigen::SparseMatrix<double>, Eigen::COLAMDOrdering<int>> ;
		using DenseMatrix = Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic>;

		INDEX a1 = mMeshFaces[0][0] ;
		INDEX a2 = mMeshFaces[0][1] ;

		if switch_once (TRUE) {
			for (auto &&i : _RANGE_ (0 ,mMeshFaces.length())) {
				const auto& i11x = mMeshFaces[i] ;
				INDEX mtempa1 = 0 ;
				ARRAY2<INDEX> mTempVe = ARRAY2<INDEX> {VAR_NONE ,VAR_NONE} ;
				for (auto &&j : _RANGE_(0 ,3)) {
					if(mBoundary[i11x[j]]){
						mTempVe[mtempa1] = i11x[j] ;
						mtempa1 ++ ;
					}
					if (mtempa1 == 2)
						break ;
				}
				
				if (mtempa1 == 2) {
					a1 = mTempVe[0] ;
					a2 = mTempVe[1] ;
					break ;
				}
				
			}
		}

		ARRAY2<VALXA> p1 = ARRAY2<VALXA> {0 ,0} ;
		ARRAY2<VALXA> p2 = ARRAY2<VALXA> {0 ,0} ;

		const auto r1x = mMeshVertices[a2] - mMeshVertices[a1] ;
		// const auto r2x = r1x.magnitude () * 2 ;
		const auto r3x = r1x.normalize () ;
		
		if switch_once (TRUE) {
			INDEX max_id = VAR_NONE ;
			VALXA temV = -VALXA_INF ;
			for (auto &&i : _RANGE_(0 ,3)) {
				if (temV < r3x[i]) {
					temV = r3x[i] ;
					max_id = i ;
				}

			}
			if (max_id == 0) {
				p2[0] = 1 ;
			}
			else {
				p2[1] = 1 ;
			}
		}

		
		auto numFaces = mMeshFaces.length () ;
		auto numVerts = mMeshVertices.length () ;
		auto numFixed = 2 ;
        auto numFree = numVerts - numFixed ;

		std::map<std::size_t, std::size_t> freeIdxTable;
		for (auto &&i : _RANGE_(0 ,mMeshVertices.length())) {
			if (i == a1 || i == a2) 
				continue ;
			auto newIdx = freeIdxTable.size() ;
            freeIdxTable[i] = newIdx ;
		}

		std::vector<Triplet> tripletsB;
        tripletsB.emplace_back (0, 0, p1[0]) ;
        tripletsB.emplace_back (1, 0, p1[1]) ;
        tripletsB.emplace_back (2, 0, p2[0]) ;
        tripletsB.emplace_back (3, 0, p2[1]) ;
        SpMat bFixed(2 * numFixed, 1) ;
        bFixed.reserve(tripletsB.size());
        bFixed.setFromTriplets(tripletsB.begin(), tripletsB.end());

		std::vector<Triplet> tripletsA ;
        tripletsB.clear () ;
		for (auto &&i : _RANGE_(0 ,mMeshFaces.length())) {
			auto e0 = 0 ;
			auto e1 = 1 ;
			auto e2 = 2 ;

			auto sin0 = MathProc::sin(solution(i ,e0)) ;
            auto sin1 = MathProc::sin(solution(i ,e1)) ;
            auto sin2 = MathProc::sin(solution(i ,e2)) ;

			std::vector<double> sins{sin0, sin1, sin2};
            auto sinMaxElem = std::max_element(sins.begin(), sins.end()) ;
            auto sinMaxIdx = std::distance(sins.begin(), sinMaxElem) ;

			if (sinMaxIdx == 0) {
                auto temp = e0;
                e0 = e1;
                e1 = e2;
                e2 = temp;
                sin0 = sins[1];
                sin1 = sins[2];
                sin2 = sins[0];
            } 
			else if (sinMaxIdx == 1) {
                auto temp = e2;
                e2 = e1;
                e1 = e0;
                e0 = temp;
                sin0 = sins[2];
                sin1 = sins[0];
                sin2 = sins[1];
            }

			auto c0 = MathProc::cos(solution(i ,e0)) ;

			auto ratio = (sin2 == VALXA(0.0)) ? VALXA(1) : sin1 / sin2 ;
            auto cosine = c0 * ratio ;
            auto sine = sin0 * ratio ;

			auto row = 2 * i ;

			// v1
            if (mMeshFaces[i][e0] == a1) {
                tripletsB.emplace_back(row, 0, cosine - VALXA(1));
                tripletsB.emplace_back(row, 1, -sine);
                tripletsB.emplace_back(row + 1, 0, sine);
                tripletsB.emplace_back(row + 1, 1, cosine - VALXA(1));
            } 
			else if (mMeshFaces[i][e0] == a2) {
                tripletsB.emplace_back(row, 2, cosine - VALXA(1));
                tripletsB.emplace_back(row, 3, -sine);
                tripletsB.emplace_back(row + 1, 2, sine);
                tripletsB.emplace_back(row + 1, 3, cosine - VALXA(1));
            } 
			else {
                auto freeIdx = freeIdxTable.at(mMeshFaces[i][e0]);
                tripletsA.emplace_back(row, 2 * freeIdx, cosine - VALXA(1));
                tripletsA.emplace_back(row, 2 * freeIdx + 1, -sine);
                tripletsA.emplace_back(row + 1, 2 * freeIdx, sine);
                tripletsA.emplace_back(row + 1, 2 * freeIdx + 1, cosine - VALXA(1));
            }

			// v2
            if (mMeshFaces[i][e1] == a1) {
                tripletsB.emplace_back(row, 0, -cosine);
                tripletsB.emplace_back(row, 1, sine);
                tripletsB.emplace_back(row + 1, 0, -sine);
                tripletsB.emplace_back(row + 1, 1, -cosine);
            } 
			else if (mMeshFaces[i][e1] == a2) {
                tripletsB.emplace_back(row, 2, -cosine);
                tripletsB.emplace_back(row, 3, sine);
                tripletsB.emplace_back(row + 1, 2, -sine);
                tripletsB.emplace_back(row + 1, 3, -cosine);
            } 
			else {
                auto freeIdx = freeIdxTable.at(mMeshFaces[i][e1]);
                tripletsA.emplace_back(row, 2 * freeIdx, -cosine);
                tripletsA.emplace_back(row, 2 * freeIdx + 1, sine);
                tripletsA.emplace_back(row + 1, 2 * freeIdx, -sine);
                tripletsA.emplace_back(row + 1, 2 * freeIdx + 1, -cosine);
            }
			
			// v3
			if (mMeshFaces[i][e2] == a1) {
                tripletsB.emplace_back(row, 0, VALXA(1));
                tripletsB.emplace_back(row + 1, 1, VALXA(1));
            } 
			else if (mMeshFaces[i][e2] == a2) {
                tripletsB.emplace_back(row, 2, VALXA(1));
                tripletsB.emplace_back(row + 1, 3, VALXA(1));
            } 
			else {
                auto freeIdx = freeIdxTable.at(mMeshFaces[i][e2]);
                tripletsA.emplace_back(row, 2 * freeIdx, VALXA(1));
                tripletsA.emplace_back(row + 1, 2 * freeIdx + 1, VALXA(1));
            }
		}
		
		SpMat A(2 * numFaces, 2 * numFree);
        A.reserve(tripletsA.size());
        A.setFromTriplets(tripletsA.begin(), tripletsA.end());

        SpMat bFree(2 * numFaces, 2 * numFixed);
        bFree.reserve(tripletsB.size());
        bFree.setFromTriplets(tripletsB.begin(), tripletsB.end());

        // 从自由矩阵和固定矩阵计算rhs
        SpMat b = bFree * bFixed * - 1;

		SpMat AtA = A.transpose() * A;
        AtA.makeCompressed();
		
        Solver solver ;
        solver.compute(AtA);
        if (solver.info() != Eigen::ComputationInfo::Success) {
            Singleton<ConsoleService>::instance ().error (_PCSTR_ ("lscm_compute fail")) ;
			Singleton<ConsoleService>::instance ().error (_PCSTR_ ("info ") ,solver.info()) ;
        }

        // Setup Atb
        SpMat Atb = A.transpose() * b ;

        // Solve AtAx = AtAb
       DenseMatrix x = solver.solve(Atb);


		Array<ARRAY2<VALXA>> tex_coord = Array<ARRAY2<VALXA>> (mMeshVertices.length()) ;
		VALXA dx = 0 ;
		VALXA minx = VALXA_INF ;
		VALXA maxx = -VALXA_INF ;

		VALXA dy = 0 ;
		VALXA miny = VALXA_INF ;
		VALXA maxy = -VALXA_INF ;

		for (auto &&i : _RANGE_(0, mMeshVertices.length())) {
			if (i == a1 || i == a2) {
                continue;
            }
			auto newIdx = 2 * freeIdxTable.at(i);
            const auto i11x  = x(newIdx, 0);
            const auto i12x  = x(newIdx + 1, 0);
			tex_coord[i] = ARRAY2<VALXA> {i11x ,i12x} ;
			minx = MathProc::minof (minx ,i11x) ;
			maxx = MathProc::maxof (maxx ,i11x) ;
			miny = MathProc::minof (miny ,i12x) ;
			maxy = MathProc::maxof (maxy ,i12x) ;
		}
		tex_coord[a1] = p1 ;
		tex_coord[a2] = p2 ;
		minx = MathProc::minof (minx ,p1[0]) ;
		maxx = MathProc::maxof (maxx ,p1[0]) ;
		miny = MathProc::minof (miny ,p1[1]) ;
		maxy = MathProc::maxof (maxy ,p1[1]) ;

		minx = MathProc::minof (minx ,p2[0]) ;
		maxx = MathProc::maxof (maxx ,p2[0]) ;
		miny = MathProc::minof (miny ,p2[1]) ;
		maxy = MathProc::maxof (maxy ,p2[1]) ;

		dx = maxx - minx ;
		dy = maxy - miny ;

		for (auto &&i : _RANGE_(0, tex_coord.length())) {
			tex_coord[i][0] = (tex_coord[i][0] - minx) * MathProc::inverse (dx) ;
			tex_coord[i][1] = (tex_coord[i][1] - miny) * MathProc::inverse (dy) ;

		}
		return _MOVE_ (tex_coord) ;
	}

	/**	LSCM 算法
	 * (一)	假设三维模型的每个三角形有一个局部正交坐标系，那么三角形的３个顶点的局部坐标为(x1 ,y1) ,(x2 ,y2) (x3 ,y3),
	 * 		法向是沿着Ｚ轴方向．相邻两个三角形的方向一至．
	 * (二) 假设U:(x ,y) --> (u ,v) 是从三维空间到二维平面的映射，那么在局部坐标系里，柯西－黎曼方程离散化为：
	 * 		aX / au - i * aX / av = 0
	 * (三) 如果Ｘ用复数表示为Ｘ＝x + iy ,让Ｕ＝u + iv.那么根据反函数导数定理得出：
	 * 		aU / ax + i * aU / ay = 0
	*/
} ;