C++ JacobiSVD::compute方法代码示例

bool MatrixXr_pseudoInverse(const MatrixXr &a, MatrixXr &a_pinv, double epsilon) {

    // see : http://en.wikipedia.org/wiki/Moore-Penrose_pseudoinverse#The_general_case_and_the_SVD_method
    if ( a.rows()<a.cols() ) return false;

    // SVD
    Eigen::JacobiSVD<MatrixXr> svdA;
    svdA.compute(a,Eigen::ComputeThinU|Eigen::ComputeThinV);
    MatrixXr vSingular = svdA.singularValues();

    // Build a diagonal matrix with the Inverted Singular values
    // The pseudo inverted singular matrix is easy to compute :
    // is formed by replacing every nonzero entry by its reciprocal (inversing).
    VectorXr vPseudoInvertedSingular(svdA.matrixV().cols(),1);

    for (int iRow =0; iRow<vSingular.rows(); iRow++) {
        if(fabs(vSingular(iRow))<=epsilon) vPseudoInvertedSingular(iRow,0)=0.;
        else vPseudoInvertedSingular(iRow,0)=1./vSingular(iRow);
    }

    // A little optimization here
    MatrixXr mAdjointU = svdA.matrixU().adjoint().block(0,0,vSingular.rows(),svdA.matrixU().adjoint().cols());

    // Pseudo-Inversion : V * S * U'
    a_pinv = (svdA.matrixV() *  vPseudoInvertedSingular.asDiagonal()) * mAdjointU;

    return true;
}

void ctms_decompositions()
{
  const int maxSize = 16;
  const int size    = 12;

  typedef Eigen::Matrix<Scalar,
                        Eigen::Dynamic, Eigen::Dynamic,
                        0,
                        maxSize, maxSize> Matrix;

  typedef Eigen::Matrix<Scalar,
                        Eigen::Dynamic, 1,
                        0,
                        maxSize, 1> Vector;

  typedef Eigen::Matrix<std::complex<Scalar>,
                        Eigen::Dynamic, Eigen::Dynamic,
                        0,
                        maxSize, maxSize> ComplexMatrix;

  const Matrix A(Matrix::Random(size, size));
  const ComplexMatrix complexA(ComplexMatrix::Random(size, size));
  const Matrix saA = A.adjoint() * A;

  // Cholesky module
  Eigen::LLT<Matrix>  LLT;  LLT.compute(A);
  Eigen::LDLT<Matrix> LDLT; LDLT.compute(A);

  // Eigenvalues module
  Eigen::HessenbergDecomposition<ComplexMatrix> hessDecomp;        hessDecomp.compute(complexA);
  Eigen::ComplexSchur<ComplexMatrix>            cSchur(size);      cSchur.compute(complexA);
  Eigen::ComplexEigenSolver<ComplexMatrix>      cEigSolver;        cEigSolver.compute(complexA);
  Eigen::EigenSolver<Matrix>                    eigSolver;         eigSolver.compute(A);
  Eigen::SelfAdjointEigenSolver<Matrix>         saEigSolver(size); saEigSolver.compute(saA);
  Eigen::Tridiagonalization<Matrix>             tridiag;           tridiag.compute(saA);

  // LU module
  Eigen::PartialPivLU<Matrix> ppLU; ppLU.compute(A);
  Eigen::FullPivLU<Matrix>    fpLU; fpLU.compute(A);

  // QR module
  Eigen::HouseholderQR<Matrix>        hQR;  hQR.compute(A);
  Eigen::ColPivHouseholderQR<Matrix>  cpQR; cpQR.compute(A);
  Eigen::FullPivHouseholderQR<Matrix> fpQR; fpQR.compute(A);

  // SVD module
  Eigen::JacobiSVD<Matrix> jSVD; jSVD.compute(A, ComputeFullU | ComputeFullV);
}

int main() {

  std::ifstream file;
  file.open("SVD_benchmark");
  if (!file) 
  {
    CGAL_TRACE_STREAM << "Error loading file!\n";
    return 0;
  }

  int ite = 200000;
  Eigen::JacobiSVD<Eigen::Matrix3d> svd;
  Eigen::Matrix3d u, v, cov, r;         
  Eigen::Vector3d w;   

  int matrix_idx = rand()%200;
  for (int i = 0; i < matrix_idx; i++)
  {
    for (int j = 0; j < 3; j++)
    {
      for (int k = 0; k < 3; k++)
      {
        file >> cov(j, k);
      }
    }
  }


  CGAL::Timer task_timer; 

  CGAL_TRACE_STREAM << "Start SVD decomposition...";
  task_timer.start();
  for (int i = 0; i < ite; i++)
  {
    
    svd.compute( cov, Eigen::ComputeFullU | Eigen::ComputeFullV );
    u = svd.matrixU(); v = svd.matrixV(); w = svd.singularValues();
    r = v*u.transpose();
  }
  task_timer.stop();
  file.close();

  CGAL_TRACE_STREAM << "done: " << task_timer.time() << "s\n";

  return 0;
}

// Derived from code by Yohann Solaro ( http://listengine.tuxfamily.org/lists.tuxfamily.org/eigen/2010/01/msg00187.html )
// see : http://en.wikipedia.org/wiki/Moore-Penrose_pseudoinverse#The_general_case_and_the_SVD_method
Eigen::MatrixXd pinv( const Eigen::MatrixXd &b, double rcond )
{
  // TODO: Figure out why it wants fewer rows than columns
  // if ( a.rows()<a.cols() )
  // return false;
  bool flip = false;
  Eigen::MatrixXd a;
  if( a.rows() < a.cols() )
  {
    a = b.transpose();
    flip = true;
  }
  else
    a = b;
  // SVD
  Eigen::JacobiSVD<Eigen::MatrixXd> svdA;
  svdA.compute( a, Eigen::ComputeFullU | Eigen::ComputeThinV );
  Eigen::JacobiSVD<Eigen::MatrixXd>::SingularValuesType vSingular = svdA.singularValues();
  // Build a diagonal matrix with the Inverted Singular values
  // The pseudo inverted singular matrix is easy to compute :
  // is formed by replacing every nonzero entry by its reciprocal (inversing).
  Eigen::VectorXd vPseudoInvertedSingular( svdA.matrixV().cols() );

  for (int iRow=0; iRow<vSingular.rows(); iRow++)
  {
    if ( fabs(vSingular(iRow)) <= rcond )
    {
      vPseudoInvertedSingular(iRow)=0.;
    }
    else
      vPseudoInvertedSingular(iRow)=1./vSingular(iRow);
  }
  // A little optimization here
  Eigen::MatrixXd mAdjointU = svdA.matrixU().adjoint().block( 0, 0, vSingular.rows(), svdA.matrixU().adjoint().cols() );
  // Pseudo-Inversion : V * S * U'
  Eigen::MatrixXd a_pinv = (svdA.matrixV() * vPseudoInvertedSingular.asDiagonal()) * mAdjointU;
  if( flip )
  {
    a = a.transpose();
    a_pinv = a_pinv.transpose();
  }
  return a_pinv;
}

IGL_INLINE void igl::min_quad_dense_precompute(
  const Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>& A,
  const Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>& Aeq,    
  const bool use_lu_decomposition,
  Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>& S)
{
  typedef Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> Mat;
        // This threshold seems to matter a lot but I'm not sure how to
        // set it
  const T treshold = igl::FLOAT_EPS;
  //const T treshold = igl::DOUBLE_EPS;

  const int n = A.rows();
  assert(A.cols() == n);
  const int m = Aeq.rows();
  assert(Aeq.cols() == n);

  // Lagrange multipliers method:
  Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> LM(n + m, n + m);
  LM.block(0, 0, n, n) = A;
  LM.block(0, n, n, m) = Aeq.transpose();
  LM.block(n, 0, m, n) = Aeq;
  LM.block(n, n, m, m).setZero();

  Mat LMpinv;
  if(use_lu_decomposition)
  {
    // if LM is close to singular, use at your own risk :)
    LMpinv = LM.inverse();
  }else
  {
    // use SVD
    typedef Eigen::Matrix<T, Eigen::Dynamic, 1> Vec; 
    Vec singValues;
    Eigen::JacobiSVD<Mat> svd;
    svd.compute(LM, Eigen::ComputeFullU | Eigen::ComputeFullV );
    const Mat& u = svd.matrixU();
    const Mat& v = svd.matrixV();
    const Vec& singVals = svd.singularValues();

    Vec pi_singVals(n + m);
    int zeroed = 0;
    for (int i=0; i<n + m; i++)
    {
      T sv = singVals(i, 0);
      assert(sv >= 0);      
                 // printf("sv: %lg ? %lg\n",(double) sv,(double)treshold);
      if (sv > treshold) pi_singVals(i, 0) = T(1) / sv;
      else 
      {
        pi_singVals(i, 0) = T(0);
        zeroed++;
      }
    }

    printf("min_quad_dense_precompute: %i singular values zeroed (threshold = %e)\n", zeroed, treshold);
    Eigen::DiagonalMatrix<T, Eigen::Dynamic> pi_diag(pi_singVals);

    LMpinv = v * pi_diag * u.transpose();
  }
  S = LMpinv.block(0, 0, n, n + m);

  //// debug:
  //mlinit(&g_pEngine);
  //
  //mlsetmatrix(&g_pEngine, "A", A);
  //mlsetmatrix(&g_pEngine, "Aeq", Aeq);
  //mlsetmatrix(&g_pEngine, "LM", LM);
  //mlsetmatrix(&g_pEngine, "u", u);
  //mlsetmatrix(&g_pEngine, "v", v);
  //MatrixXd svMat = singVals;
  //mlsetmatrix(&g_pEngine, "singVals", svMat);
  //mlsetmatrix(&g_pEngine, "LMpinv", LMpinv);
  //mlsetmatrix(&g_pEngine, "S", S);

  //int hu = 1;
}

Types::Transform PlaneToPlaneCalibration::estimateTransform(const std::vector<PlanePair> & plane_pair_vector)
{
  const int size = plane_pair_vector.size();

  Eigen::MatrixXd normals_1(3, size);
  Eigen::MatrixXd normals_2(3, size);
  Eigen::VectorXd distances_1(size);
  Eigen::VectorXd distances_2(size);

  for (int i = 0; i < size; ++i)
  {
    const Types::Plane &plane_1 = plane_pair_vector[i].plane_1_;
    const Types::Plane &plane_2 = plane_pair_vector[i].plane_2_;

    if (plane_1.offset() > 0)
    {
      normals_1.col(i) = -plane_1.normal();
      distances_1(i) = plane_1.offset();
    }
    else
    {
      normals_1.col(i) = plane_1.normal();
      distances_1(i) = -plane_1.offset();
    }

    if (plane_2.offset() > 0)
    {
      normals_2.col(i) = -plane_2.normal();
      distances_2(i) = plane_2.offset();
    }
    else
    {
      normals_2.col(i) = plane_2.normal();
      distances_2(i) = -plane_2.offset();
    }
  }

//  std::cout << normals_1 << std::endl;
//  std::cout << distances_1.transpose() << std::endl;
//  std::cout << normals_2 << std::endl;
//  std::cout << distances_2.transpose() << std::endl;

  Eigen::Matrix3d USV = normals_2 * normals_1.transpose();
  Eigen::JacobiSVD<Eigen::Matrix3d> svd;
  svd.compute(USV, Eigen::ComputeFullU | Eigen::ComputeFullV);

  Types::Pose pose;
  pose.translation() = (normals_1 * normals_1.transpose()).inverse() * normals_1 * (distances_1 - distances_2);
  pose.linear() = svd.matrixV() * svd.matrixU().transpose();

//  // Point-Plane Constraints
//
//  // Initial system (Eq. 10)
//
//  Eigen::MatrixXd system(size * 3, 13);
//  for (int i = 0; i < size; ++i)
//  {
//    const double d = plane_pair_vector[i].plane_1_.offset();
//    const Eigen::Vector3d n = plane_pair_vector[i].plane_1_.normal();
//
//    const Point3d x = -plane_pair_vector[i].plane_2_.normal() * plane_pair_vector[i].plane_2_.offset();
//
//    Eigen::Matrix3d X(Eigen::Matrix3d::Random());
//    while (X.determinant() < 1e-5)
//      X = Eigen::Matrix3d::Random();
//
//    const Eigen::Vector3d & n2 = plane_pair_vector[i].plane_2_.normal();
////    X.col(1)[2] = -n2.head<2>().dot(X.col(1).head<2>()) / n2[2];
////    X.col(2)[2] = -n2.head<2>().dot(X.col(2).head<2>()) / n2[2];
//    X.col(1) -= n2 * n2.dot(X.col(1));
//    X.col(2) -= n2 * n2.dot(X.col(2));
//
//    X.col(0) = x;
//    X.col(1) += x;
//    X.col(2) += x;
//
//    for (int j = 0; j < 3; ++j)
//    {
//      const Point3d & x = X.col(j);
//
//      system.row(i + size * j) << d + n[0] * x[0] + n[1] * x[1] + n[2] * x[2],  // q_0^2
//                                  2 * n[2] * x[1] - 2 * n[1] * x[2],            // q_0 * q_1
//                                  2 * n[0] * x[2] - 2 * n[2] * x[0],            // q_0 * q_2
//                                  2 * n[1] * x[0] - 2 * n[0] * x[1],            // q_0 * q_3
//                                  d + n[0] * x[0] - n[1] * x[1] - n[2] * x[2],  // q_1^2
//                                  2 * n[0] * x[1] + 2 * n[1] * x[0],            // q_1 * q_2
//                                  2 * n[0] * x[2] + 2 * n[2] * x[0],            // q_1 * q_3
//                                  d - n[0] * x[0] + n[1] * x[1] - n[2] * x[2],  // q_2^2
//                                  2 * n[1] * x[2] + 2 * n[2] * x[1],            // q_2 * q_3
//                                  d - n[0] * x[0] - n[1] * x[1] + n[2] * x[2],  // q_3^2
//                                  n[0], n[1], n[2]; // q'*q*t
//    }
//  }
//
//  //std::cout << system << std::endl;
//
//  // Gaussian reduction
//  for (int k = 0; k < 3; ++k)
//    for (int i = k + 1; i < size * 3; ++i)
//      system.row(i) = system.row(i) - system.row(k) * system.row(i)[10 + k] / system.row(k)[10 + k];
//.........这里部分代码省略.........

//.........这里部分代码省略.........
			}
		}
	}
	
	/// STEP #4: Compute new grid velocities

	//Temporary variables for plasticity and force calculation
	//We need one set of variables for each thread that will be running
	eigen_matrix3 def_elastic, def_plastic, energy, svd_u, svd_v;
	Eigen::JacobiSVD<eigen_matrix3, Eigen::NoQRPreconditioner> svd;
	eigen_vector3 svd_e;
	matrix3  HDK_def_plastic, HDK_def_elastic, HDK_energy;
	freal* data_dp = HDK_def_plastic.data();
	freal* data_de = HDK_def_elastic.data();
	freal* data_energy = HDK_energy.data();
	//Map Eigen matrices to HDK matrices
	Eigen::Map<eigen_matrix3> data_dp_map(data_dp);
	Eigen::Map<eigen_matrix3> data_de_map(data_de);
	Eigen::Map<eigen_matrix3> data_energy_map(data_energy);	

	//Compute force at each particle and transfer to Eulerian grid
	//We use "nvel" to hold the grid force, since that variable is not in use
	for (GA_Iterator it(gdp_in->getPointRange()); !it.atEnd(); it.advance()){
		int pid = it.getOffset();
		
		//Apply plasticity to deformation gradient, before computing forces
		//We need to use the Eigen lib to do the SVD; transfer houdini matrices to Eigen matrices
		HDK_def_plastic = p_Fp.get(pid);
		HDK_def_elastic = p_Fe.get(pid);
		def_plastic = Eigen::Map<eigen_matrix3>(data_dp);
		def_elastic = Eigen::Map<eigen_matrix3>(data_de);
		
		//Compute singular value decomposition (uev*)
		svd.compute(def_elastic, Eigen::ComputeFullV | Eigen::ComputeFullU);
		svd_e = svd.singularValues();
		svd_u = svd.matrixU();
		svd_v = svd.matrixV();
		//Clamp singular values
		for (int i=0; i<3; i++){
			if (svd_e[i] < CRIT_COMPRESS) 
				svd_e[i] = CRIT_COMPRESS;
			else if (svd_e[i] > CRIT_STRETCH)
				svd_e[i] = CRIT_STRETCH;
		}
		//Put SVD back together for new elastic and plastic gradients
		def_plastic = svd_v * svd_e.asDiagonal().inverse() * svd_u.transpose() * def_elastic * def_plastic;
		svd_v.transposeInPlace();
		def_elastic = svd_u * svd_e.asDiagonal() * svd_v;
		
		//Now compute the energy partial derivative (which we use to get force at each grid node)
		energy = 2*mu*(def_elastic - svd_u*svd_v)*def_elastic.transpose();
		//Je is the determinant of def_elastic (equivalent to svd_e.prod())
		freal Je = svd_e.prod(),
			contour = lambda*Je*(Je-1),
			jp = def_plastic.determinant(),
			particle_vol = p_volume.get(pid);
		for (int i=0; i<3; i++)
			energy(i,i) += contour;
		energy *=  particle_vol * exp(HARDENING*(1-jp));
		
		//Transfer Eigen matrices back to HDK
		data_dp_map = def_plastic;
		data_de_map = def_elastic;
		data_energy_map = energy;
		
		p_Fp.set(pid,HDK_def_plastic);