C++ Transform3D::SetTranslation方法代码示例

RDGeom::Transform3D *computeCanonicalTransform(const Conformer &conf,
                                               const RDGeom::Point3D *center,
                                               bool normalizeCovar,
                                               bool ignoreHs) {
  RDGeom::Point3D origin;
  if (!center) {
    origin = computeCentroid(conf, ignoreHs);
  } else {
    origin = (*center);
  }
  RDNumeric::DoubleSymmMatrix *covMat =
      computeCovarianceMatrix(conf, origin, normalizeCovar, ignoreHs);
  // find the eigen values and eigen vectors for the covMat
  RDNumeric::DoubleMatrix eigVecs(3, 3);
  RDNumeric::DoubleVector eigVals(3);
  // if we have a single atom system we don't need to do anyhting other than
  // setting translation
  // translation
  unsigned int nAtms = conf.getNumAtoms();
  RDGeom::Transform3D *trans = new RDGeom::Transform3D;

  // set the translation
  origin *= -1.0;
  // trans->SetTranslation(origin);
  // if we have a single atom system we don't need to do anyhting setting
  // translation is sufficient
  if (nAtms > 1) {
    RDNumeric::EigenSolvers::powerEigenSolver(3, *covMat, eigVals, eigVecs,
                                              conf.getNumAtoms());
    // deal with zero eigen value systems
    unsigned int i, j, dim = 3;
    for (i = 0; i < 3; ++i) {
      if (fabs(eigVals.getVal(i)) < EIGEN_TOLERANCE) {
        dim--;
      }
    }
    CHECK_INVARIANT(dim >= 1, "");
    if (dim < 3) {
      RDGeom::Point3D first(eigVecs.getVal(0, 0), eigVecs.getVal(0, 1),
                            eigVecs.getVal(0, 2));
      if (dim == 1) {
        // pick an arbitrary eigen vector perpendicular to the first vector
        RDGeom::Point3D second(first.getPerpendicular());
        eigVecs.setVal(1, 0, second.x);
        eigVecs.setVal(1, 1, second.y);
        eigVecs.setVal(1, 2, second.z);
        if (eigVals.getVal(0) > 1.0) {
          eigVals.setVal(1, 1.0);
        } else {
          eigVals.setVal(1, eigVals.getVal(0) / 2.0);
        }
      }
      RDGeom::Point3D second(eigVecs.getVal(1, 0), eigVecs.getVal(1, 1),
                             eigVecs.getVal(1, 2));
      // pick the third eigen vector perpendicular to the first two
      RDGeom::Point3D third = first.crossProduct(second);
      eigVecs.setVal(2, 0, third.x);
      eigVecs.setVal(2, 1, third.y);
      eigVecs.setVal(2, 2, third.z);
      if (eigVals.getVal(1) > 1.0) {
        eigVals.setVal(2, 1.0);
      } else {
        eigVals.setVal(2, eigVals.getVal(1) / 2.0);
      }
    }
    // now set the transformation
    for (i = 0; i < 3; ++i) {
      for (j = 0; j < 3; ++j) {
        trans->setVal(i, j, eigVecs.getVal(i, j));
      }
    }
  }  // end of multiple atom system
  trans->TransformPoint(origin);
  trans->SetTranslation(origin);
  delete covMat;

  return trans;
}

double AlignPoints(const RDGeom::Point3DConstPtrVect &refPoints,
                   const RDGeom::Point3DConstPtrVect &probePoints,
                   RDGeom::Transform3D &trans, const DoubleVector *weights,
                   bool reflect, unsigned int maxIterations) {
  unsigned int npt = refPoints.size();
  PRECONDITION(npt == probePoints.size(), "Mismatch in number of points");
  trans.setToIdentity();
  const DoubleVector *wts;
  double wtsSum;
  bool ownWts;
  if (weights) {
    PRECONDITION(npt == weights->size(), "Mismatch in number of points");
    wts = weights;
    wtsSum = _sumOfWeights(*wts);
    ownWts = false;
  } else {
    wts = new DoubleVector(npt, 1.0);
    wtsSum = static_cast<double>(npt);
    ownWts = true;
  }

  RDGeom::Point3D rptSum = _weightedSumOfPoints(refPoints, *wts);
  RDGeom::Point3D pptSum = _weightedSumOfPoints(probePoints, *wts);

  double rptSumLenSq = _weightedSumOfLenSq(refPoints, *wts);
  double pptSumLenSq = _weightedSumOfLenSq(probePoints, *wts);

  double covMat[3][3];

  // compute the co-variance matrix
  _computeCovarianceMat(refPoints, probePoints, *wts, covMat);
  if (ownWts) {
    delete wts;
    wts = 0;
  }
  if (reflect) {
    rptSum *= -1.0;
    reflectCovMat(covMat);
  }

  // convert the covariance matrix to a 4x4 matrix that needs to be diagonalized
  double quad[4][4];
  _covertCovMatToQuad(covMat, rptSum, pptSum, wtsSum, quad);

  // get the eigenVecs and eigenVals for the matrix
  double eigenVecs[4][4], eigenVals[4];
  jacobi(quad, eigenVals, eigenVecs, maxIterations);

  // get the quaternion
  double quater[4];
  quater[0] = eigenVecs[0][0];
  quater[1] = eigenVecs[1][0];
  quater[2] = eigenVecs[2][0];
  quater[3] = eigenVecs[3][0];

  trans.SetRotationFromQuaternion(quater);
  if (reflect) {
    // put the flip in the rotation matrix
    trans.Reflect();
  }
  // compute the SSR value
  double ssr = eigenVals[0] - (pptSum.lengthSq() + rptSum.lengthSq()) / wtsSum +
               rptSumLenSq + pptSumLenSq;

  if ((ssr < 0.0) && (fabs(ssr) < TOLERANCE)) {
    ssr = 0.0;
  }
  if (reflect) {
    rptSum *= -1.0;
  }

  // set the translation
  trans.TransformPoint(pptSum);
  RDGeom::Point3D move = rptSum;
  move -= pptSum;
  move /= wtsSum;
  trans.SetTranslation(move);
  return ssr;
}