1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
|
VectorXd MyMesh::meanValueCoordinates(Vector3d p, bool &postiveity) { postiveity = true; VectorXd mvc_coor = VectorXd::Zero(this->numberOfVertices);
Matrix<double,3,Dynamic> pro_vertices;
pro_vertices.resize(3, this->numberOfVertices); for (int i = 0; i < this->numberOfVertices; ++i) { double r = (this->vertices.col(i) - p).norm(); pro_vertices(0, i) = (this->vertices(0, i) - p(0)) / r; pro_vertices(1, i) = (this->vertices(1, i) - p(1)) / r; pro_vertices(2, i) = (this->vertices(2, i) - p(2)) / r; }
Matrix<double, 3, 3> T; Vector3d n0, n1, n2; double beta0, beta1, beta2; for (int i = 0; i < this->numberOfFaces; ++i) { int i0, i1, i2; i0 = this->faces(0, i); i1 = this->faces(1, i); i2 = this->faces(2, i);
T.col(0) = pro_vertices.col(i0); T.col(1) = pro_vertices.col(i1); T.col(2) = pro_vertices.col(i2);
n0 = pro_vertices.col(i0).cross(pro_vertices.col(i1)).normalized(); n1 = pro_vertices.col(i1).cross(pro_vertices.col(i2)).normalized(); n2 = pro_vertices.col(i2).cross(pro_vertices.col(i0)).normalized();
beta0 = acos(pro_vertices.col(i0).dot(pro_vertices.col(i1))); beta1 = acos(pro_vertices.col(i1).dot(pro_vertices.col(i2))); beta2 = acos(pro_vertices.col(i2).dot(pro_vertices.col(i0)));
mvc_coor[i0] += (beta1 + beta0*n0.dot(n1) + beta2*n2.dot(n1)) / (2 * pro_vertices.col(i0).dot(n1)); mvc_coor[i1] += (beta2 + beta1*n1.dot(n2) + beta0*n0.dot(n2)) / (2 * pro_vertices.col(i1).dot(n2)); mvc_coor[i2] += (beta0 + beta2*n2.dot(n0) + beta1*n1.dot(n0)) / (2 * pro_vertices.col(i2).dot(n0)); } for (int i = 0; i < this->numberOfVertices; ++i) { mvc_coor[i] /= (this->vertices.col(i) - p).norm(); } mvc_coor /= mvc_coor.sum(); postiveity = mvc_coor.minCoeff() >= 0; return mvc_coor; }
|