Page MenuHome
Differential D5078 Diff 15940 intern/dualcon/intern/octree.cpp

Changeset View

Standalone View

View Options

intern/dualcon/intern/octree.cpp

Show First 20 Lines • Show All 2,174 Lines • ▼ Show 20 Lines for (int i = 0; i < 3; i++) {
if (getFaceEdgeNum(&node->leaf, i * 2)) { if (getFaceEdgeNum(&node->leaf, i * 2)) {
nface++; nface++;
} }
} }
} }
} }
/* from http://eigen.tuxfamily.org/bz/show_bug.cgi?id=257 */ /* from http://eigen.tuxfamily.org/bz/show_bug.cgi?id=257 */
template<typename _Matrix_Type_> static void pseudoInverse(const Eigen::Matrix3f &a, Eigen::Matrix3f &result, float tolerance)
void pseudoInverse(const _Matrix_Type_ &a,
_Matrix_Type_ &result,
double epsilon = std::numeric_limits<typename _Matrix_Type_::Scalar>::epsilon())
{ {
Eigen::JacobiSVD<_Matrix_Type_> svd = a.jacobiSvd(Eigen::ComputeFullU | Eigen::ComputeFullV); Eigen::JacobiSVD<Eigen::Matrix3f> svd = a.jacobiSvd(Eigen::ComputeFullU | Eigen::ComputeFullV);
typename _Matrix_Type_::Scalar tolerance = epsilon * std::max(a.cols(), a.rows()) *
svd.singularValues().array().abs().maxCoeff();
result = svd.matrixV() * result = svd.matrixV() *
_Matrix_Type_((svd.singularValues().array().abs() > tolerance) Eigen::Vector3f((svd.singularValues().array().abs() > tolerance)
.select(svd.singularValues().array().inverse(), 0)) .select(svd.singularValues().array().inverse(), 0))
.asDiagonal() * .asDiagonal() *
svd.matrixU().adjoint(); svd.matrixU().adjoint();
} }
static void solve_least_squares(const float halfA[], static void solve_least_squares(const float halfA[],
const float b[], const float b[],
const float midpoint[], const float midpoint[],
float rvalue[]) float rvalue[])
{ {
/* calculate pseudo-inverse */ /* calculate pseudo-inverse */
Eigen::MatrixXf A(3, 3), pinv(3, 3); Eigen::Matrix3f A, pinv;
A << halfA[0], halfA[1], halfA[2], halfA[1], halfA[3], halfA[4], halfA[2], halfA[4], halfA[5]; A << halfA[0], halfA[1], halfA[2], halfA[1], halfA[3], halfA[4], halfA[2], halfA[4], halfA[5];
pseudoInverse(A, pinv); pseudoInverse(A, pinv, 0.1f);
Eigen::Vector3f b2(b), mp(midpoint), result; Eigen::Vector3f b2(b), mp(midpoint), result;
b2 = b2 + A * -mp; b2 = b2 + A * -mp;
result = pinv * b2 + mp; result = pinv * b2 + mp;
for (int i = 0; i < 3; i++) for (int i = 0; i < 3; i++)
rvalue[i] = result(i); rvalue[i] = result(i);
} }
▲ Show 20 Lines • Show All 683 Lines • Show Last 20 Lines