C++ MatrixXd::isIdentity方法代码示例

/**
 * Constructor using std::vector (for python exposition purposes)
 *
 * @param[in] p base::problem to be rotated
 * @param[in] rotation std::vector<std::vector<double> > expressing the problem rotation
 *
 * @see problem::base constructors.
 */
rotated::rotated(const base &p,
				 const std::vector<std::vector<double> > &rotation):
		base_meta(
		 p,
		 p.get_dimension(),
		 p.get_i_dimension(),
		 p.get_f_dimension(),
		 p.get_c_dimension(),
		 p.get_ic_dimension(),
		 p.get_c_tol()),
	m_Rotate(),m_normalize_translation(), m_normalize_scale()
{
	if(!(rotation.size()==get_dimension())){
			pagmo_throw(value_error,"The input matrix dimensions seem incorrect");
	}
	if(p.get_i_dimension()>0){
		pagmo_throw(value_error,"Input problem has an integer dimension. Cannot rotate it.");
	}
	m_Rotate.resize(rotation.size(),rotation.size());
	for (base::size_type i = 0; i < rotation.size(); ++i) {
		if(!(rotation.size()==rotation[i].size())){
			pagmo_throw(value_error,"The input matrix seems not to be square");
		}
		for (base::size_type j = 0; j < rotation[i].size(); ++j) {
			m_Rotate(i,j) = rotation[i][j];
		}
	}
	m_InvRotate = m_Rotate.transpose();
	
	Eigen::MatrixXd check = m_InvRotate * m_Rotate;
	if(!check.isIdentity(1e-5)){
		pagmo_throw(value_error,"The input matrix seems not to be orthonormal (to a tolerance of 1e-5)");
	}
	configure_new_bounds();
}

//.........这里部分代码省略.........
      input_header.size(3) == (int) Math::SH::NforL (Math::SH::LforN (input_header.size(3)))) {
    CONSOLE ("SH series detected, performing apodised PSF reorientation");
    fod_reorientation = true;

    Eigen::MatrixXd directions_az_el;
    opt = get_options ("directions");
    if (opt.size())
      directions_az_el = load_matrix (opt[0][0]);
    else
      directions_az_el = DWI::Directions::electrostatic_repulsion_300();
    Math::SH::spherical2cartesian (directions_az_el, directions_cartesian);

    // load with SH coeffients contiguous in RAM
    stride = Stride::contiguous_along_axis (3, input_header);
  }

  // Modulate FODs
  bool modulate = false;
  if (get_options ("modulate").size()) {
    modulate = true;
    if (!fod_reorientation)
      throw Exception ("modulation can only be performed with FOD reorientation");
  }

  // Rotate/Flip gradient directions if present
  if (linear && input_header.ndim() == 4 && !warp && !fod_reorientation) {
    try {
      auto grad = DWI::get_DW_scheme (input_header);
      if (input_header.size(3) == (ssize_t) grad.rows()) {
        INFO ("DW gradients detected and will be reoriented");
        Eigen::MatrixXd rotation = linear_transform.linear();
        Eigen::MatrixXd test = rotation.transpose() * rotation;
        test = test.array() / test.diagonal().mean();
        if (!test.isIdentity (0.001))
        WARN ("the input linear transform contains shear or anisotropic scaling and "
              "therefore should not be used to reorient diffusion gradients");
        if (replace)
          rotation = linear_transform.linear() * input_header.transform().linear().inverse();
        for (ssize_t n = 0; n < grad.rows(); ++n) {
          Eigen::Vector3 grad_vector = grad.block<1,3>(n,0);
          grad.block<1,3>(n,0) = rotation * grad_vector;
        }
        DWI::set_DW_scheme (output_header, grad);
      }
    }
    catch (Exception& e) {
      e.display (2);
    }
  }

  // Interpolator
  int interp = 2;  // cubic
  opt = get_options ("interp");
  if (opt.size()) {
    interp = opt[0][0];
    if (!warp && !template_header)
      WARN ("interpolator choice ignored since the input image will not be regridded");
  }

  // Out of bounds value
  float out_of_bounds_value = 0.0;
  opt = get_options ("nan");
  if (opt.size()) {
    out_of_bounds_value = NAN;
    if (!warp && !template_header)
      WARN ("Out of bounds value ignored since the input image will not be regridded");