Commit 7db26a5d by Kenneth Vanhoey

### Corrections to matrix creation for lightfield metric computation.

```Optimization of matrix calculation (exploit fact that 2 of them are
symmetric)```
parent 2b452e36
 ... ... @@ -40,8 +40,6 @@ // Eigen includes #include #define CONST_VAL -5212368.54127 // random value /*! \namespace CGoGN * \brief namespace for all elements composing the CGoGN library */ ... ... @@ -563,12 +561,22 @@ public: static std::vector coefsFromTensors(Geom::Tensor3d* T) ; /*! * \brief method to complete a symmetric tensor that was * only filled in its first half * \brief method to complete a symmetric matrix that was * only filled in its first half (line >= column) * * \param T the tensor to fill * \param M the matrix to fill */ static void fillTensor(Geom::Tensor3d& T) ; static void fillSymmetricMatrix(Eigen::MatrixXd& M) ; /*! * \brief method to rotate a tensor representing a polynomial light field * * \param T the tensor representing a polynomial light field * \param R the 3x3 matrix representing a rotation in (u,v,1)-space * * \return a new rotated tensor representing a polynomial light field. */ static Geom::Tensor3d rotate(const Geom::Tensor3d& T, const Geom::Matrix33d& R) ; private: // Double computation is crucial for stability ... ... @@ -595,7 +603,7 @@ private: */ bool optimize(std::vector& coefs) const ; Geom::Tensor3d rotate(const Geom::Tensor3d& T, const Geom::Matrix33d& R) ; // Geom::Tensor3d rotate(const Geom::Tensor3d& T, const Geom::Matrix33d& R) ; /*! * \brief method to build a rotate matrix (rotation in tangent plane) ... ... @@ -639,6 +647,9 @@ private: * \return the integral of product of monomes */ Eigen::MatrixXd buildIntegralMatrix_C(const REAL& alpha, unsigned int size) ; Eigen::MatrixXd buildLowerLeftIntegralMatrix_A(const REAL& alpha, unsigned int size) ; Eigen::MatrixXd buildLowerLeftIntegralMatrix_C(const REAL& alpha, unsigned int size) ; } ; } // Utils ... ...
 ... ... @@ -344,18 +344,6 @@ QuadricHF::QuadricHF(const Geom::Tensor3d* T, const REAL& gamma, const REA Eigen::MatrixXd integ_b = buildIntegralMatrix_B(alpha,nbcoefs) ; // Parameterized integral matrix b Eigen::MatrixXd integ_c = buildIntegralMatrix_C(alpha,nbcoefs) ; // Parameterized integral matrix c // for (unsigned int i = 0 ; i < nbcoefs ; ++i) // for (unsigned int j = 0 ; j < nbcoefs ; ++j) // std::cout << "(" << i << "," << j << ")=" << m_A(i,j) << std::endl ; // // for (unsigned int i = 0 ; i < nbcoefs ; ++i) // for (unsigned int j = 0 ; j < nbcoefs ; ++j) // std::cout << "(" << i << "," << j << ")=" << integ_b(i,j) << std::endl ; // // for (unsigned int i = 0 ; i < nbcoefs ; ++i) // for (unsigned int j = 0 ; j < nbcoefs ; ++j) // std::cout << "(" << i << "," << j << ")=" << integ_c(i,j) << std::endl ; for (unsigned int c = 0 ; c < 3 ; ++c) { Eigen::VectorXd v ; v.resize(nbcoefs) ; ... ... @@ -369,11 +357,7 @@ QuadricHF::QuadricHF(const Geom::Tensor3d* T, const REAL& gamma, const REA template QuadricHF::~QuadricHF() { //delete m_A ; //for (unsigned int c = 0 ; c < 3 ; ++c) // delete m_b[c] ; } {} template void ... ... @@ -462,17 +446,6 @@ QuadricHF::operator /= (const REAL& v) return *this ; } // REAL // QuadricHF::operator() (const VECNp& v) const // { // VECN hv ; // for (unsigned int i = 0 ; i < 3 ; ++i) // hv[i] = v[i] ; // // return evaluate(v) ; //} template REAL QuadricHF::operator() (const std::vector& coefs) const ... ... @@ -484,16 +457,6 @@ template bool QuadricHF::findOptimizedCoefs(std::vector& coefs) { // unsigned int nbcoefs = 6 ; // for (unsigned int i = 0 ; i < nbcoefs ; ++i) // for (unsigned int j = 0 ; j < nbcoefs ; ++j) // std::cout << "A(" << i << "," << j << ")=" << m_A(i,j) << std::endl ; // // for (unsigned int i = 0 ; i < nbcoefs ; ++i) // std::cout << "b(" << i << ")=" << m_b[0][i] << std::endl ; // // std::cout << "c=" << m_c[0] << std::endl ; return optimize(coefs) ; } ... ... @@ -523,7 +486,7 @@ QuadricHF::optimize(std::vector& coefs) const { coefs.resize(m_b[0].size()) ; if (fabs(m_A.determinant()) < 1e-6 ) if (fabs(m_A.determinant()) < 1e-10 ) return false ; Eigen::MatrixXd Ainv = m_A.inverse() ; ... ... @@ -578,171 +541,183 @@ QuadricHF::rotate(const Geom::Tensor3d& T, const Geom::Matrix33d& R) return Tp ; } template void QuadricHF::fillSymmetricMatrix(Eigen::MatrixXd& M) { assert(M.cols() == M.rows() || !"QuadricHF::fillSymmetricMatrix: matrix to fill should be a square matrix") ; for (unsigned int c = 0 ; c < M.cols() ; ++c) for (unsigned int l = 0 ; l < c ; ++l) M( l, c) = M (c, l) ; } template Eigen::MatrixXd QuadricHF::buildIntegralMatrix_A(const REAL& alpha, unsigned int size) { Eigen::MatrixXd res = buildLowerLeftIntegralMatrix_A(alpha,size) ; fillSymmetricMatrix(res) ; return res ; } template Eigen::MatrixXd QuadricHF::buildIntegralMatrix_C(const REAL& alpha, unsigned int size) { Eigen::MatrixXd res = buildLowerLeftIntegralMatrix_C(alpha,size) ; fillSymmetricMatrix(res) ; return res ; } template Eigen::MatrixXd QuadricHF::buildLowerLeftIntegralMatrix_A(const REAL& alpha, unsigned int size) { Eigen::MatrixXd A(size,size) ; A( 0 , 0 ) = 2*(M_PI-alpha) ; A( 0 , 1 ) = M_PI*sin(alpha)/2.0 ; A( 0 , 2 ) = 0 ; A( 0 , 3 ) = 0 ; A( 0 , 4 ) = -(sin(2*alpha)+2*alpha-2*M_PI)/3.0 ; A( 0 , 5 ) = 2.0*(M_PI-alpha)/3.0 ; A( 1 , 0 ) = M_PI*sin(alpha)/2.0 ; A( 1 , 1 ) = -(sin(2*alpha)+2*alpha-2*M_PI)/3.0 ; A( 1 , 2 ) = 0 ; A( 1 , 3 ) = 0 ; A( 1 , 4 ) = -M_PI*(pow(sin(alpha),3)-3*sin(alpha))/8.0 ; A( 1 , 5 ) = M_PI*sin(alpha)/8.0 ; A( 2 , 0 ) = 0 ; A( 2 , 1 ) = 0 ; A( 2 , 2 ) = 2.0*(M_PI-alpha)/3.0 ; A( 2 , 3 ) = M_PI*sin(alpha)/8.0 ; A( 2 , 4 ) = 0 ; A( 2 , 5 ) = 0 ; A( 3 , 0 ) = 0 ; A( 3 , 1 ) = 0 ; A( 3 , 2 ) = M_PI*sin(alpha)/8.0 ; A( 3 , 3 ) = -(sin(2*alpha)+2*alpha-2*M_PI)/1.5E+1 ; A( 3 , 4 ) = 0 ; A( 3 , 5 ) = 0 ; A( 4 , 0 ) = -(sin(2*alpha)+2*alpha-2*M_PI)/3.0 ; A( 4 , 1 ) = -M_PI*(pow(sin(alpha),3)-3*sin(alpha))/8.0 ; A( 4 , 2 ) = 0 ; A( 4 , 3 ) = 0 ; A( 4 , 4 ) = -(sin(4*alpha)+8*sin(2*alpha)+12*alpha-12*M_PI)/3.0E+1 ; A( 4 , 5 ) = -(sin(2*alpha)+2*alpha-2*M_PI)/1.5E+1 ; A( 5 , 0 ) = 2.0*(M_PI-alpha)/3.0 ; A( 5 , 1 ) = M_PI*sin(alpha)/8.0 ; A( 5 , 2 ) = 0 ; A( 5 , 3 ) = 0 ; A( 5 , 4 ) = -(sin(2*alpha)+2*alpha-2*M_PI)/1.5E+1 ; A( 5 , 5 ) = 2.0*(M_PI-alpha)/5.0 ; A( 0 , 0 ) = 2*(M_PI-alpha) ; A( 1 , 0 ) = M_PI*sin(alpha)/2.0 ; A( 1 , 1 ) = -(sin(2*alpha)+2*alpha-2*M_PI)/3.0 ; A( 2 , 0 ) = 0 ; A( 2 , 1 ) = 0 ; A( 2 , 2 ) = 2.0*(M_PI-alpha)/3.0 ; A( 3 , 0 ) = 0 ; A( 3 , 1 ) = 0 ; A( 3 , 2 ) = M_PI*sin(alpha)/8.0 ; A( 3 , 3 ) = -(sin(2*alpha)+2*alpha-2*M_PI)/1.5E+1 ; A( 4 , 0 ) = -(sin(2*alpha)+2*alpha-2*M_PI)/3.0 ; A( 4 , 1 ) = -M_PI*(pow(sin(alpha),3)-3*sin(alpha))/8.0 ; A( 4 , 2 ) = 0 ; A( 4 , 3 ) = 0 ; A( 4 , 4 ) = -(sin(4*alpha)+8*sin(2*alpha)+12*alpha-12*M_PI)/3.0E+1 ; A( 5 , 0 ) = 2.0*(M_PI-alpha)/3.0 ; A( 5 , 1 ) = M_PI*sin(alpha)/8.0 ; A( 5 , 2 ) = 0 ; A( 5 , 3 ) = 0 ; A( 5 , 4 ) = -(sin(2*alpha)+2*alpha-2*M_PI)/1.5E+1 ; A( 5 , 5 ) = 2.0*(M_PI-alpha)/5.0 ; if (size < 7) return A ; A( 6 , 0 ) = 0 ; A( 6 , 1 ) = 0 ; A( 6 , 2 ) = -(sin(2*alpha)+2*alpha-2*M_PI)/1.5E+1 ; A( 6 , 3 ) = -M_PI*(pow(sin(alpha),3)-3*sin(alpha))/4.8E+1 ; A( 6 , 4 ) = 0 ; A( 6 , 5 ) = 0 ; A( 6 , 6 ) = -(sin(4*alpha)+8*sin(2*alpha)+12*alpha-12*M_PI)/2.1E+2 ; A( 6 , 7 ) = 0 ; A( 6 , 8 ) = 0 ; A( 6 , 9 ) = -(sin(2*alpha)+2*alpha-2*M_PI)/3.5E+1 ; A( 7 , 0 ) = M_PI*sin(alpha)/8.0 ; A( 7 , 1 ) = -(sin(2*alpha)+2*alpha-2*M_PI)/1.5E+1 ; A( 7 , 2 ) = 0 ; A( 7 , 3 ) = 0 ; A( 7 , 4 ) = -M_PI*(pow(sin(alpha),3)-3*sin(alpha))/4.8E+1 ; A( 7 , 5 ) = M_PI*sin(alpha)/1.6E+1 ; A( 7 , 6 ) = 0 ; A( 7 , 7 ) = -(sin(2*alpha)+2*alpha-2*M_PI)/3.5E+1 ; A( 7 , 8 ) = -(sin(4*alpha)+8*sin(2*alpha)+12*alpha-12*M_PI)/2.1E+2 ; A( 7 , 9 ) = 0 ; A( 8 , 0 ) = -M_PI*(pow(sin(alpha),3)-3*sin(alpha))/8.0 ; A( 8 , 1 ) = -(sin(4*alpha)+8*sin(2*alpha)+12*alpha-12*M_PI)/3.0E+1 ; A( 8 , 2 ) = 0 ; A( 8 , 3 ) = 0 ; A( 8 , 4 ) = M_PI*(3*pow(sin(alpha),5)-10*pow(sin(alpha),3)+15*sin(alpha))/4.8E+1 ; A( 8 , 5 ) = -M_PI*(pow(sin(alpha),3)-3*sin(alpha))/4.8E+1 ; A( 8 , 6 ) = 0 ; A( 8 , 7 ) = -(sin(4*alpha)+8*sin(2*alpha)+12*alpha-12*M_PI)/2.1E+2 ; A( 8 , 8 ) = -(9*sin(4*alpha)-4*pow(sin(2*alpha),3)+48*sin(2*alpha)+60*alpha-60*M_PI )/2.1E+2 ; A( 8 , 9 ) = 0 ; A( 9 , 0 ) = 0 ; A( 9 , 1 ) = 0 ; A( 9 , 2 ) = 2.0*(M_PI-alpha)/5.0 ; A( 9 , 3 ) = M_PI*sin(alpha)/1.6E+1 ; A( 9 , 4 ) = 0 ; A( 9 , 5 ) = 0 ; A( 9 , 6 ) = -(sin(2*alpha)+2*alpha-2*M_PI)/3.5E+1 ; A( 9 , 7 ) = 0 ; A( 9 , 8 ) = 0 ; A( 9 , 9 ) = 2.0*(M_PI-alpha)/7.0 ; A( 6 , 0 ) = 0 ; A( 6 , 1 ) = 0 ; A( 6 , 2 ) = -(sin(2*alpha)+2*alpha-2*M_PI)/1.5E+1 ; A( 6 , 3 ) = -M_PI*(pow(sin(alpha),3)-3*sin(alpha))/4.8E+1 ; A( 6 , 4 ) = 0 ; A( 6 , 5 ) = 0 ; A( 6 , 6 ) = -(sin(4*alpha)+8*sin(2*alpha)+12*alpha-12*M_PI)/2.1E+2 ; A( 7 , 0 ) = M_PI*sin(alpha)/8.0 ; A( 7 , 1 ) = -(sin(2*alpha)+2*alpha-2*M_PI)/1.5E+1 ; A( 7 , 2 ) = 0 ; A( 7 , 3 ) = 0 ; A( 7 , 4 ) = -M_PI*(pow(sin(alpha),3)-3*sin(alpha))/4.8E+1 ; A( 7 , 5 ) = M_PI*sin(alpha)/1.6E+1 ; A( 7 , 6 ) = 0 ; A( 7 , 7 ) = -(sin(2*alpha)+2*alpha-2*M_PI)/3.5E+1 ; A( 8 , 0 ) = -M_PI*(pow(sin(alpha),3)-3*sin(alpha))/8.0 ; A( 8 , 1 ) = -(sin(4*alpha)+8*sin(2*alpha)+12*alpha-12*M_PI)/3.0E+1 ; A( 8 , 2 ) = 0 ; A( 8 , 3 ) = 0 ; A( 8 , 4 ) = M_PI*(3*pow(sin(alpha),5)-10*pow(sin(alpha),3)+15*sin(alpha))/4.8E+1 ; A( 8 , 5 ) = -M_PI*(pow(sin(alpha),3)-3*sin(alpha))/4.8E+1 ; A( 8 , 6 ) = 0 ; A( 8 , 7 ) = -(sin(4*alpha)+8*sin(2*alpha)+12*alpha-12*M_PI)/2.1E+2 ; A( 8 , 8 ) = -(9*sin(4*alpha)-4*pow(sin(2*alpha),3)+48*sin(2*alpha)+60*alpha-60*M_PI )/2.1E+2 ; A( 9 , 0 ) = 0 ; A( 9 , 1 ) = 0 ; A( 9 , 2 ) = 2.0*(M_PI-alpha)/5.0 ; A( 9 , 3 ) = M_PI*sin(alpha)/1.6E+1 ; A( 9 , 4 ) = 0 ; A( 9 , 5 ) = 0 ; A( 9 , 6 ) = -(sin(2*alpha)+2*alpha-2*M_PI)/3.5E+1 ; A( 9 , 7 ) = 0 ; A( 9 , 8 ) = 0 ; A( 9 , 9 ) = 2.0*(M_PI-alpha)/7.0 ; if (size < 11) return A ; A( 10 , 0 ) = -(sin(2*alpha)+2*alpha-2*M_PI)/1.5E+1 ; A( 10 , 1 ) = -M_PI*(pow(sin(alpha),3)-3*sin(alpha))/4.8E+1 ; A( 10 , 2 ) = 0 ; A( 10 , 3 ) = 0 ; A( 10 , 4 ) = -(sin(4*alpha)+8*sin(2*alpha)+12*alpha-12*M_PI)/2.1E+2 ; A( 10 , 5 ) = -(sin(2*alpha)+2*alpha-2*M_PI)/3.5E+1 ; A( 10 , 6 ) = 0 ; A( 10 , 7 ) = -M_PI*(pow(sin(alpha),3)-3*sin(alpha))/1.28E+2 ; A( 10 , 8 ) = M_PI*(3*pow(sin(alpha),5)-10*pow(sin(alpha),3)+15*sin(alpha))/3.84E+2 ; A( 10 , 9 ) = 0 ; A( 10 , 10 ) = -(sin(4*alpha)+8*sin(2*alpha)+12*alpha-12*M_PI)/6.3E+2 ; A( 10 , 11 ) = 0 ; A( 10 , 12 ) = 0 ; A( 10 , 13 ) = -(9*sin(4*alpha)-4*pow(sin(2*alpha),3)+48*sin(2*alpha)+60*alpha-60*M_PI )/1.89E+3 ; A( 10 , 14 ) = -(sin(2*alpha)+2*alpha-2*M_PI)/6.3E+1 ; A( 11 , 0 ) = 0 ; A( 11 , 1 ) = 0 ; A( 11 , 2 ) = -M_PI*(pow(sin(alpha),3)-3*sin(alpha))/4.8E+1 ; A( 11 , 3 ) = -(sin(4*alpha)+8*sin(2*alpha)+12*alpha-12*M_PI)/2.1E+2 ; A( 11 , 4 ) = 0 ; A( 11 , 5 ) = 0 ; A( 11 , 6 ) = M_PI*(3*pow(sin(alpha),5)-10*pow(sin(alpha),3)+15*sin(alpha))/3.84E+2 ; A( 11 , 7 ) = 0 ; A( 11 , 8 ) = 0 ; A( 11 , 9 ) = -M_PI*(pow(sin(alpha),3)-3*sin(alpha))/1.28E+2 ; A( 11 , 10 ) = 0 ; A( 11 , 11 ) = -(9*sin(4*alpha)-4*pow(sin(2*alpha),3)+48*sin(2*alpha)+60*alpha-60*M_PI )/1.89E+3 ; A( 11 , 12 ) = -(sin(4*alpha)+8*sin(2*alpha)+12*alpha-12*M_PI)/6.3E+2 ; A( 11 , 13 ) = 0 ; A( 11 , 14 ) = 0 ; A( 12 , 0 ) = 0 ; A( 12 , 1 ) = 0 ; A( 12 , 2 ) = M_PI*sin(alpha)/1.6E+1 ; A( 12 , 3 ) = -(sin(2*alpha)+2*alpha-2*M_PI)/3.5E+1 ; A( 12 , 4 ) = 0 ; A( 12 , 5 ) = 0 ; A( 12 , 6 ) = -M_PI*(pow(sin(alpha),3)-3*sin(alpha))/1.28E+2 ; A( 12 , 7 ) = 0 ; A( 12 , 8 ) = 0 ; A( 12 , 9 ) = 5.0*M_PI*sin(alpha)/1.28E+2 ; A( 12 , 10 ) = 0 ; A( 12 , 11 ) = -(sin(4*alpha)+8*sin(2*alpha)+12*alpha-12*M_PI)/6.3E+2 ; A( 12 , 12 ) = -(sin(2*alpha)+2*alpha-2*M_PI)/6.3E+1 ; A( 12 , 13 ) = 0 ; A( 12 , 14 ) = 0 ; A( 13 , 0 ) = -(sin(4*alpha)+8*sin(2*alpha)+12*alpha-12*M_PI)/3.0E+1 ; A( 13 , 1 ) = M_PI*(3*pow(sin(alpha),5)-10*pow(sin(alpha),3)+15*sin(alpha))/4.8E+1 ; A( 13 , 2 ) = 0 ; A( 13 , 3 ) = 0 ; A( 13 , 4 ) = -(9*sin(4*alpha)-4*pow(sin(2*alpha),3)+48*sin(2*alpha)+60*alpha-60*M_PI )/2.1E+2 ; A( 13 , 5 ) = -(sin(4*alpha)+8*sin(2*alpha)+12*alpha-12*M_PI)/2.1E+2 ; A( 13 , 6 ) = 0 ; A( 13 , 7 ) = M_PI*(3*pow(sin(alpha),5)-10*pow(sin(alpha),3)+15*sin(alpha))/3.84E+2 ; A( 13 , 8 ) = -M_PI*(5*pow(sin(alpha),7)-21*pow(sin(alpha),5)+35*pow(sin(alpha),3)-35*sin(alpha))/1.28E+2 ; A( 13 , 9 ) = 0 ; A( 13 , 10 ) = -(9*sin(4*alpha)-4*pow(sin(2*alpha),3)+48*sin(2*alpha)+60*alpha-60*M_PI )/1.89E+3 ; A( 13 , 11 ) = 0 ; A( 13 , 12 ) = 0 ; A( 13 , 13 ) = -(3*sin(8*alpha)+168*sin(4*alpha)-128*pow(sin(2*alpha),3)+768*sin(2*alpha)+840*alpha-840*M_PI)/3.78E+3 ; A( 13 , 14 ) = -(sin(4*alpha)+8*sin(2*alpha)+12*alpha-12*M_PI)/6.3E+2 ; A( 14 , 0 ) = 2.0*(M_PI-alpha)/5.0 ; A( 14 , 1 ) = M_PI*sin(alpha)/1.6E+1 ; A( 14 , 2 ) = 0 ; A( 14 , 3 ) = 0 ; A( 14 , 4 ) = -(sin(2*alpha)+2*alpha-2*M_PI)/3.5E+1 ; A( 14 , 5 ) = 2.0*(M_PI-alpha)/7.0 ; A( 14 , 6 ) = 0 ; A( 14 , 7 ) = 5.0*M_PI*sin(alpha)/1.28E+2 ; A( 14 , 8 ) = -M_PI*(pow(sin(alpha),3)-3*sin(alpha))/1.28E+2 ; A( 14 , 9 ) = 0 ; A( 14 , 10 ) = -(sin(2*alpha)+2*alpha-2*M_PI)/6.3E+1 ; A( 14 , 11 ) = 0 ; A( 14 , 12 ) = 0 ; A( 14 , 13 ) = -(sin(4*alpha)+8*sin(2*alpha)+12*alpha-12*M_PI)/6.3E+2 ; A( 14 , 14 ) = 2.0*(M_PI-alpha)/9.0 ; A( 10 , 0 ) = -(sin(2*alpha)+2*alpha-2*M_PI)/1.5E+1 ; A( 10 , 1 ) = -M_PI*(pow(sin(alpha),3)-3*sin(alpha))/4.8E+1 ; A( 10 , 2 ) = 0 ; A( 10 , 3 ) = 0 ; A( 10 , 4 ) = -(sin(4*alpha)+8*sin(2*alpha)+12*alpha-12*M_PI)/2.1E+2 ; A( 10 , 5 ) = -(sin(2*alpha)+2*alpha-2*M_PI)/3.5E+1 ; A( 10 , 6 ) = 0 ; A( 10 , 7 ) = -M_PI*(pow(sin(alpha),3)-3*sin(alpha))/1.28E+2 ; A( 10 , 8 ) = M_PI*(3*pow(sin(alpha),5)-10*pow(sin(alpha),3)+15*sin(alpha))/3.84E+2 ; A( 10 , 9 ) = 0 ; A( 10 , 10 ) = -(sin(4*alpha)+8*sin(2*alpha)+12*alpha-12*M_PI)/6.3E+2 ; A( 11 , 0 ) = 0 ; A( 11 , 1 ) = 0 ; A( 11 , 2 ) = -M_PI*(pow(sin(alpha),3)-3*sin(alpha))/4.8E+1 ; A( 11 , 3 ) = -(sin(4*alpha)+8*sin(2*alpha)+12*alpha-12*M_PI)/2.1E+2 ; A( 11 , 4 ) = 0 ; A( 11 , 5 ) = 0 ; A( 11 , 6 ) = M_PI*(3*pow(sin(alpha),5)-10*pow(sin(alpha),3)+15*sin(alpha))/3.84E+2 ; A( 11 , 7 ) = 0 ; A( 11 , 8 ) = 0 ; A( 11 , 9 ) = -M_PI*(pow(sin(alpha),3)-3*sin(alpha))/1.28E+2 ; A( 11 , 10 ) = 0 ; A( 11 , 11 ) = -(9*sin(4*alpha)-4*pow(sin(2*alpha),3)+48*sin(2*alpha)+60*alpha-60*M_PI )/1.89E+3 ; A( 12 , 0 ) = 0 ; A( 12 , 1 ) = 0 ; A( 12 , 2 ) = M_PI*sin(alpha)/1.6E+1 ; A( 12 , 3 ) = -(sin(2*alpha)+2*alpha-2*M_PI)/3.5E+1 ; A( 12 , 4 ) = 0 ; A( 12 , 5 ) = 0 ; A( 12 , 6 ) = -M_PI*(pow(sin(alpha),3)-3*sin(alpha))/1.28E+2 ; A( 12 , 7 ) = 0 ; A( 12 , 8 ) = 0 ; A( 12 , 9 ) = 5.0*M_PI*sin(alpha)/1.28E+2 ; A( 12 , 10 ) = 0 ; A( 12 , 11 ) = -(sin(4*alpha)+8*sin(2*alpha)+12*alpha-12*M_PI)/6.3E+2 ; A( 12 , 12 ) = -(sin(2*alpha)+2*alpha-2*M_PI)/6.3E+1 ; A( 13 , 0 ) = -(sin(4*alpha)+8*sin(2*alpha)+12*alpha-12*M_PI)/3.0E+1 ; A( 13 , 1 ) = M_PI*(3*pow(sin(alpha),5)-10*pow(sin(alpha),3)+15*sin(alpha))/4.8E+1 ; A( 13 , 2 ) = 0 ; A( 13 , 3 ) = 0 ; A( 13 , 4 ) = -(9*sin(4*alpha)-4*pow(sin(2*alpha),3)+48*sin(2*alpha)+60*alpha-60*M_PI )/2.1E+2 ; A( 13 , 5 ) = -(sin(4*alpha)+8*sin(2*alpha)+12*alpha-12*M_PI)/2.1E+2 ; A( 13 , 6 ) = 0 ; A( 13 , 7 ) = M_PI*(3*pow(sin(alpha),5)-10*pow(sin(alpha),3)+15*sin(alpha))/3.84E+2 ; A( 13 , 8 ) = -M_PI*(5*pow(sin(alpha),7)-21*pow(sin(alpha),5)+35*pow(sin(alpha),3)-35*sin(alpha))/1.28E+2 ; A( 13 , 9 ) = 0 ; A( 13 , 10 ) = -(9*sin(4*alpha)-4*pow(sin(2*alpha),3)+48*sin(2*alpha)+60*alpha-60*M_PI )/1.89E+3 ; A( 13 , 11 ) = 0 ; A( 13 , 12 ) = 0 ; A( 13 , 13 ) = -(3*sin(8*alpha)+168*sin(4*alpha)-128*pow(sin(2*alpha),3)+768*sin(2*alpha)+840*alpha-840*M_PI)/3.78E+3 ; A( 14 , 0 ) = 2.0*(M_PI-alpha)/5.0 ; A( 14 , 1 ) = M_PI*sin(alpha)/1.6E+1 ; A( 14 , 2 ) = 0 ; A( 14 , 3 ) = 0 ; A( 14 , 4 ) = -(sin(2*alpha)+2*alpha-2*M_PI)/3.5E+1 ; A( 14 , 5 ) = 2.0*(M_PI-alpha)/7.0 ; A( 14 , 6 ) = 0 ; A( 14 , 7 ) = 5.0*M_PI*sin(alpha)/1.28E+2 ; A( 14 , 8 ) = -M_PI*(pow(sin(alpha),3)-3*sin(alpha))/1.28E+2 ; A( 14 , 9 ) = 0 ; A( 14 , 10 ) = -(sin(2*alpha)+2*alpha-2*M_PI)/6.3E+1 ; A( 14 , 11 ) = 0 ; A( 14 , 12 ) = 0 ; A( 14 , 13 ) = -(sin(4*alpha)+8*sin(2*alpha)+12*alpha-12*M_PI)/6.3E+2 ; A( 14 , 14 ) = 2.0*(M_PI-alpha)/9.0 ; return A ; } ... ... @@ -753,362 +728,421 @@ QuadricHF::buildIntegralMatrix_B(const REAL& alpha, unsigned int size) { Eigen::MatrixXd B(size,size) ; B( 0 , 0 ) = 2*(M_PI-alpha) ; B( 0 , 1 ) = M_PI*(sin(2*alpha)+sin(alpha))/2.0 ; B( 0 , 2 ) = 0 ; B( 0 , 3 ) = 0 ; B( 0 , 4 ) = 4.0*((sin(2*alpha)+2*alpha)/4.0-(sin(4*alpha)+4*alpha-2*M_PI) /4.0)/3.0 ; B( 0 , 5 ) = 2.0*(M_PI-alpha)/3.0 ; B( 1 , 0 ) = M_PI*sin(alpha)/2.0 ; B( 1 , 1 ) = 4.0*(sin(alpha)/4.0-(sin(3*alpha)+(2*alpha-2*M_PI)*cos(alpha) )/4.0)/3.0 ; B( 1 , 2 ) = 0 ; B( 1 , 3 ) = 0 ; B( 1 , 4 ) = 3.0*M_PI*((sin(5*alpha)+3*sin(3*alpha)+6*sin(alpha))/1.2E+1+sin( 2*alpha)/3.0)/8.0 ; B( 1 , 5 ) = M_PI*sin(alpha)/8.0 ; B( 2 , 0 ) = 0 ; B( 2 , 1 ) = 0 ; B( 2 , 2 ) = 2.0*(M_PI-alpha)/3.0 ; B( 2 , 3 ) = M_PI*(sin(2*alpha)+sin(alpha))/8.0 ; B( 2 , 4 ) = 0 ; B( 2 , 5 ) = 0 ; B( 3 , 0 ) = 0 ; B( 3 , 1 ) = 0 ; B( 3 , 2 ) = M_PI*sin(alpha)/8.0 ; B( 3 , 3 ) = 4.0*(sin(alpha)/4.0-(sin(3*alpha)+(2*alpha-2*M_PI)*cos(alpha) )/4.0)/1.5E+1 ; B( 3 , 4 ) = 0 ; B( 3 , 5 ) = 0 ; B( 4 , 0 ) = -(sin(2*alpha)+2*alpha-2*M_PI)/3.0 ; B( 4 , 1 ) = 3.0*M_PI*((sin(4*alpha)+6*sin(2*alpha))/1.2E+1+sin(alpha)/3.0 )/8.0 ; B( 4 , 2 ) = 0 ; B( 4 , 3 ) = 0 ; B( 4 , 4 ) = 1.6E+1*(5.0*sin(2*alpha)/3.2E+1-(sin(6*alpha)+4*sin(4*alpha)+4* sin(2*alpha)+(4*alpha-4*M_PI)*cos(2*alpha)+8*alpha-8*M_PI)/3.2E+1 )/1.5E+1 ; B( 4 , 5 ) = -(sin(2*alpha)+2*alpha-2*M_PI)/1.5E+1 ; B( 5 , 0 ) = 2.0*(M_PI-alpha)/3.0 ; B( 5 , 1 ) = M_PI*(sin(2*alpha)+sin(alpha))/8.0 ; B( 5 , 2 ) = 0 ; B( 5 , 3 ) = 0 ; B( 5 , 4 ) = 4.0*((sin(2*alpha)+2*alpha)/4.0-(sin(4*alpha)+4*alpha-2*M_PI) /4.0)/1.5E+1 ; B( 5 , 5 ) = 2.0*(M_PI-alpha)/5.0 ; if (size < 7) return B ; B( 6 , 0 ) = 0 ; B( 6 , 1 ) = 0 ; B( 6 , 2 ) = -(sin(2*alpha)+2*alpha-2*M_PI)/1.5E+1 ; B( 6 , 3 ) = M_PI*((sin(4*alpha)+6*sin(2*alpha))/1.2E+1+sin(alpha)/3.0)/1.6E+1 ; B( 6 , 4 ) = 0 ; B( 6 , 5 ) = 0 ; B( 6 , 6 ) = 1.6E+1*(5.0*sin(2*alpha)/3.2E+1-(sin(6*alpha)+4*sin(4*alpha)+4* sin(2*alpha)+(4*alpha-4*M_PI)*cos(2*alpha)+8*alpha-8*M_PI)/3.2E+1 )/1.05E+2 ; B( 6 , 7 ) = 0 ; B( 6 , 8 ) = 0 ; B( 6 , 9 ) = -(sin(2*alpha)+2*alpha-2*M_PI)/3.5E+1 ; B( 7 , 0 ) = M_PI*sin(alpha)/8.0 ; B( 7 , 1 ) = 4.0*(sin(alpha)/4.0-(sin(3*alpha)+(2*alpha-2*M_PI)*cos(alpha) )/4.0)/1.5E+1 ; B( 7 , 2 ) = 0 ; B( 7 , 3 ) = 0 ; B( 7 , 4 ) = M_PI*((sin(5*alpha)+3*sin(3*alpha)+6*sin(alpha))/1.2E+1+sin(2*alpha )/3.0)/1.6E+1 ; B( 7 , 5 ) = M_PI*sin(alpha)/1.6E+1 ; B( 7 , 6 ) = 0 ; B( 7 , 7 ) = 4.0*(sin(alpha)/4.0-(sin(3*alpha)+(2*alpha-2*M_PI)*cos(alpha) )/4.0)/3.5E+1 ; B( 7 , 8 ) = 1.6E+1*((3*sin(3*alpha)+6*sin(alpha))/3.2E+1-(sin(7*alpha)+2*sin(5 *alpha)+6*sin(3*alpha)+(12*alpha-12*M_PI)*cos(alpha))/3.2E+1)/1.05E+2 ; B( 7 , 9 ) = 0 ; B( 8 , 0 ) = -M_PI*(pow(sin(alpha),3)-3*sin(alpha))/8.0 ; B( 8 , 1 ) = 1.6E+1*(5.0*sin(alpha)/3.2E+1-(sin(5*alpha)+6*sin(3*alpha)+2*sin(alpha)+(12*alpha-12*M_PI)*cos(alpha))/3.2E+1)/1.5E+1 ; B( 8 , 2 ) = 0 ; B( 8 , 3 ) = 0 ; B( 8 , 4 ) = 5.0*M_PI*((3*sin(7*alpha)+15*sin(5*alpha)+55*sin(3*alpha)+75*sin (alpha))/2.4E+2+sin(2*alpha)/5.0)/1.6E+1 ; B( 8 , 5 ) = -M_PI*(pow(sin(alpha),3)-3*sin(alpha))/4.8E+1 ; B( 8 , 6 ) = 0 ; B( 8 , 7 ) = 1.6E+1*(5.0*sin(alpha)/3.2E+1-(sin(5*alpha)+6*sin(3*alpha)+2*sin(alpha)+(12*alpha-12*M_PI)*cos(alpha))/3.2E+1)/1.05E+2 ; B( 8 , 8 ) = 3.2E+1*((29*sin(3*alpha)+45*sin(alpha))/3.84E+2-(2*sin(9*alpha)+9*sin(7*alpha)+27*sin(5*alpha)+54*sin(3*alpha)+(12*alpha-12*M_PI)*cos(3*alpha)+18*sin(alpha)+(108*alpha-108*M_PI)*cos(alpha))/3.84E+2)/3.5E+1 ; B( 8 , 9 ) = 0 ; B( 9 , 0 ) = 0 ; B( 9 , 1 ) = 0 ; B( 9 , 2 ) = 2.0*(M_PI-alpha)/5.0 ; B( 9 , 3 ) = M_PI*(sin(2*alpha)+sin(alpha))/1.6E+1 ; B( 9 , 4 ) = 0 ; B( 9 , 5 ) = 0 ; B( 9 , 6 ) = 4.0*((sin(2*alpha)+2*alpha)/4.0-(sin(4*alpha)+4*alpha-2*M_PI) /4.0)/3.5E+1 ; B( 9 , 7 ) = 0 ; B( 9 , 8 ) = 0 ; B( 9 , 9 ) = 2.0*(M_PI-alpha)/7.0 ; if (size < 11) return B ; B( 10 , 0 ) = -(sin(2*alpha)+2*alpha-2*M_PI)/1.5E+1 ; B( 10 , 1 ) = M_PI*((sin(4*alpha)+6*sin(2*alpha))/1.2E+1+sin(alpha)/3.0)/1.6E+1 ; B( 10 , 2 ) = 0 ; B( 10 , 3 ) = 0 ; B( 10 , 4 ) = 1.6E+1*(5.0*sin(2*alpha)/3.2E+1-(sin(6*alpha)+4*sin(4*alpha)+4* sin(2*alpha)+(4*alpha-4*M_PI)*cos(2*alpha)+8*alpha-8*M_PI)/3.2E+1 )/1.05E+2 ; B( 10 , 5 ) = -(sin(2*alpha)+2*alpha-2*M_PI)/3.5E+1 ; B( 10 , 6 ) = 0 ; B( 10 , 7 ) = 3.0*M_PI*((sin(4*alpha)+6*sin(2*alpha))/1.2E+1+sin(alpha)/3.0 )/1.28E+2 ; B( 10 , 8 ) = 5.0*M_PI*((3*sin(8*alpha)+10*sin(6*alpha)+30*sin(4*alpha)+90*sin (2*alpha))/2.4E+2+(7*sin(3*alpha)+15*sin(alpha))/6.0E+1)/1.28E+2 ; B( 10 , 9 ) = 0 ; B( 10 , 10 ) = 1.6E+1*(5.0*sin(2*alpha)/3.2E+1-(sin(6*alpha)+4*sin(4*alpha)+4* sin(2*alpha)+(4*alpha-4*M_PI)*cos(2*alpha)+8*alpha-8*M_PI)/3.2E+1 )/3.15E+2 ; B( 10 , 11 ) = 0 ; B( 10 , 12 ) = 0 ; B( 10 , 13 ) = 3.2E+1*((7*sin(4*alpha)+30*sin(2*alpha))/1.92E+2-(sin(10*alpha)+3* sin(8*alpha)+9*sin(6*alpha)+24*sin(4*alpha)+18*sin(2*alpha)+(24 *alpha-24*M_PI)*cos(2*alpha)+36*alpha-36*M_PI)/1.92E+2)/3.15E+2 ; B( 10 , 14 ) = -(sin(2*alpha)+2*alpha-2*M_PI)/6.3E+1 ; B( 11 , 0 ) = 0 ; B( 11 , 1 ) = 0 ; B( 11 , 2 ) = -M_PI*(pow(sin(alpha),3)-3*sin(alpha))/4.8E+1 ; B( 11 , 3 ) = 1.6E+1*(5.0*sin(alpha)/3.2E+1-(sin(5*alpha)+6*sin(3*alpha)+2*sin(alpha)+(12*alpha-12*M_PI)*cos(alpha))/3.2E+1)/1.05E+2 ; B( 11 , 4 ) = 0 ; B( 11 , 5 ) = 0 ; B( 11 , 6 ) = 5.0*M_PI*((3*sin(7*alpha)+15*sin(5*alpha)+55*sin(3*alpha)+75*sin (alpha))/2.4E+2+sin(2*alpha)/5.0)/1.28E+2 ; B( 11 , 7 ) = 0 ; B( 11 , 8 ) = 0 ; B( 11 , 9 ) = -M_PI*(pow(sin(alpha),3)-3*sin(alpha))/1.28E+2 ; B( 11 , 10 ) = 0 ; B( 11 , 11 ) = 3.2E+1*((29*sin(3*alpha)+45*sin(alpha))/3.84E+2-(2*sin(9*alpha)+9*sin(7*alpha)+27*sin(5*alpha)+54*sin(3*alpha)+(12*alpha-12*M_PI)* cos(3*alpha)+18*sin(alpha)+(108*alpha-108*M_PI)*cos(alpha))/3.84E+2)/3.15E+2 ; B( 11 , 12 ) = 1.6E+1*(5.0*sin(alpha)/3.2E+1-(sin(5*alpha)+6*sin(3*alpha)+2*sin(alpha)+(12*alpha-12*M_PI)*cos(alpha))/3.2E+1)/3.15E+2 ; B( 11 , 13 ) = 0 ; B( 11 , 14 ) = 0 ; B( 12 , 0 ) = 0 ; B( 12 , 1 ) = 0 ; B( 12 , 2 ) = M_PI*sin(alpha)/1.6E+1 ; B( 12 , 3 ) = 4.0*(sin(alpha)/4.0-(sin(3*alpha)+(2*alpha-2*M_PI)*cos(alpha) )/4.0)/3.5E+1 ; B( 12 , 4 ) = 0 ; B( 12 , 5 ) = 0 ; B( 12 , 6 ) = 3.0*M_PI*((sin(5*alpha)+3*sin(3*alpha)+6*sin(alpha))/1.2E+1+sin( 2*alpha)/3.0)/1.28E+2 ; B( 12 , 7 ) = 0 ; B( 12 , 8 ) = 0 ; B( 12 , 9 ) = 5.0*M_PI*sin(alpha)/1.28E+2 ; B( 12 , 10 ) = 0 ; B( 12 , 11 ) = 1.6E+1*((3*sin(3*alpha)+6*sin(alpha))/3.2E+1-(sin(7*alpha)+2*sin(5 *alpha)+6*sin(3*alpha)+(12*alpha-12*M_PI)*cos(alpha))/3.2E+1)/3.15E+2 ; B( 12 , 12 ) = 4.0*(sin(alpha)/4.0-(sin(3*alpha)+(2*alpha-2*M_PI)*cos(alpha) )/4.0)/6.3E+1 ; B( 12 , 13 ) = 0 ; B( 12 , 14 ) = 0 ; B( 13 , 0 ) = -(sin(4*alpha)+8*sin(2*alpha)+12*alpha-12*M_PI)/3.0E+1 ; B( 13 , 1 ) = 5.0*M_PI*((3*sin(6*alpha)+20*sin(4*alpha)+95*sin(2*alpha))/2.4E+2 + sin(alpha)/5.0)/1.6E+1 ; B( 13 , 2 ) = 0 ; B( 13 , 3 ) = 0 ; B( 13 , 4 ) = 3.2E+1*(1.1E+1*sin(2*alpha)/9.6E+1-(sin(8*alpha)+6*sin(6*alpha)+21 *sin(4*alpha)+24*sin(2*alpha)+(24*alpha-24*M_PI)*cos(2*alpha)+36 *alpha-36*M_PI)/1.92E+2)/3.5E+1 ; B( 13 , 5 ) = -(sin(4*alpha)+8*sin(2*alpha)+12*alpha-12*M_PI)/2.1E+2 ; B( 13 , 6 ) = 0 ; B( 13 , 7 ) = 5.0*M_PI*((3*sin(6*alpha)+20*sin(4*alpha)+95*sin(2*alpha))/2.4E+2+sin(alpha)/5.0)/1.28E+2 ; B( 13 , 8 ) = 3.5E+1*M_PI*((5*sin(10*alpha)+28*sin(8*alpha)+91*sin(6*alpha)+280*sin(4*alpha)+630*sin(2*alpha))/2.24E+3+(13*sin(3*alpha)+21*sin(alpha))/1.4E+2)/1.28E+2 ; B( 13 , 9 ) = 0 ; B( 13 , 10 ) = 3.2E+1*(1.1E+1*sin(2*alpha)/9.6E+1-(sin(8*alpha)+6*sin(6*alpha)+21 *sin(4*alpha)+24*sin(2*alpha)+(24*alpha-24*M_PI)*cos(2*alpha)+36 *alpha-36*M_PI)/1.92E+2)/3.15E+2 ; B( 13 , 11 ) = 0 ; B( 13 , 12 ) = 0 ; B( 13 , 13 ) = 2.56E+2*((103*sin(4*alpha)+352*sin(2*alpha))/3.072E+3-(3*sin(12*alpha)+16*sin(10*alpha)+52*sin(8*alpha)+144*sin(6*alpha)+324*sin(4*alpha)+(24*alpha-24*M_PI)*cos(4*alpha)+288*sin(2*alpha)+(384*alpha-384*M_PI)*cos(2*alpha)+432*alpha-432*M_PI)/3.072E+3)/3.15E+2 ; B( 13 , 14 ) = -(sin(4*alpha)+8*sin(2*alpha)+12*alpha-12*M_PI)/6.3E+2 ; B( 14 , 0 ) = 2.0*(M_PI-alpha)/5.0 ; B( 14 , 1 ) = M_PI*(sin(2*alpha)+sin(alpha))/1.6E+1 ; B( 14 , 2 ) = 0 ; B( 14 , 3 ) = 0 ; B( 14 , 4 ) = 4.0*((sin(2*alpha)+2*alpha)/4.0-(sin(4*alpha)+4*alpha-2*M_PI) /4.0)/3.5E+1 ; B( 14 , 5 ) = 2.0*(M_PI-alpha)/7.0 ; B( 14 , 6 ) = 0 ; B( 14 , 7 ) = 5.0*M_PI*(sin(2*alpha)+sin(alpha))/1.28E+2 ; B( 14 , 8 ) = 3.0*M_PI*(-(pow(sin(2*alpha),3)-3*sin(2*alpha))/3.0-(pow(sin(alpha),3)-3*sin(alpha))/3.0)/1.28E+2 ; B( 14 , 9 ) = 0 ; B( 14 , 10 ) = 4.0*((sin(2*alpha)+2*alpha)/4.0-(sin(4*alpha)+4*alpha-2*M_PI) /4.0)/6.3E+1 ; B( 14 , 11 ) = 0 ; B( 14 , 12 ) = 0 ; B( 14 , 13 ) = 1.6E+1*((sin(4*alpha)+8*sin(2*alpha)+12*alpha)/3.2E+1-(sin(8*alpha )+8*sin(4*alpha)+24*alpha-12*M_PI)/3.2E+1)/3.15E+2 ; B( 14 , 14 ) = 2.0*(M_PI-alpha)/9.0 ; B( 0 , 0 ) = 2*(M_PI-alpha) ; B( 0 , 1 ) = M_PI*(sin(2*alpha)+sin(alpha))/2.0 ; B( 0 , 2 ) = 0 ; B( 0 , 3 ) = 0 ; B( 0 , 4 ) = 4.0*((sin(2*alpha)+2*alpha)/4.0-(sin(4*alpha)+4*alpha-2*M_PI) /4.0)/3.0 ; B( 0 , 5 ) = 2.0*(M_PI-alpha)/3.0 ; B( 1 , 0 ) = M_PI*sin(alpha)/2.0 ; B( 1 , 1 ) = 4.0*(sin(alpha)/4.0-(sin(3*alpha)+(2*alpha-2*M_PI)*cos(alpha) )/4.0)/3.0 ; B( 1 , 2 ) = 0 ; B( 1 , 3 ) = 0 ; B( 1 , 4 ) = 3.0*M_PI*((sin(5*alpha)+3*sin(3*alpha)+6*sin(alpha))/1.2E+1+sin( 2*alpha)/3.0)/8.0 ; B( 1 , 5 ) = M_PI*sin(alpha)/8.0 ; B( 2 , 0 ) = 0 ; B( 2 , 1 ) = 0 ; B( 2 , 2 ) = 2.0*(M_PI-alpha)/3.0 ; B( 2 , 3 ) = M_PI*(sin(2*alpha)+sin(alpha))/8.0 ; B( 2 , 4 ) = 0 ; B( 2 , 5 ) = 0 ; B( 3 , 0 ) = 0 ; B( 3 , 1 ) = 0 ; B( 3 , 2 ) = M_PI*sin(alpha)/8.0 ; B( 3 , 3 ) = 4.0*(sin(alpha)/4.0-(sin(3*alpha)+(2*alpha-2*M_PI)*cos(alpha) )/4.0)/1.5E+1 ; B( 3 , 4 ) = 0 ; B( 3 , 5 ) = 0 ; B( 4 , 0 ) = -(sin(2*alpha)+2*alpha-2*M_PI)/3.0 ; B( 4 , 1 ) = 3.0*M_PI*((sin(4*alpha)+6*sin(2*alpha))/1.2E+1+sin(alpha)/3.0 )/8.0 ; B( 4 , 2 ) = 0 ; B( 4 , 3 ) = 0 ; B( 4 , 4 ) = 1.6E+1*(5.0*sin(2*alpha)/3.2E+1-(sin(6*alpha)+4*sin(4*alpha)+4* sin(2*alpha)+(4*alpha-4*M_PI)*cos(2*alpha)+8*alpha-8*M_PI)/3.2E+1 )/1.5E+1 ; B( 4 , 5 ) = -(sin(2*alpha)+2*alpha-2*M_PI)/1.5E+1 ;