|
|
|
|
|
#include "nodeSampler.h" |
|
|
#include "tools.h" |
|
|
#include "OmpHelper.h" |
|
|
#include "geodesic_algorithm_exact.h" |
|
|
#include "nanoflann.h" |
|
|
|
|
|
namespace svr |
|
|
{ |
|
|
|
|
|
static auto square = [](const Scalar argu) { return argu * argu; }; |
|
|
static auto cube = [](const Scalar argu) { return argu * argu * argu; }; |
|
|
static auto max = [](const Scalar lhs, const Scalar rhs) { return lhs > rhs ? lhs : rhs; }; |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
Scalar nodeSampler::SampleAndConstuctAxis(Mesh &mesh, Scalar sampleRadiusRatio, sampleAxis axis) |
|
|
{ |
|
|
|
|
|
m_meshVertexNum = mesh.n_vertices(); |
|
|
m_meshEdgeNum = mesh.n_edges(); |
|
|
m_mesh = & mesh; |
|
|
|
|
|
|
|
|
for (size_t i = 0; i < m_meshEdgeNum; ++i) |
|
|
{ |
|
|
OpenMesh::EdgeHandle eh = mesh.edge_handle(i); |
|
|
Scalar edgeLen = mesh.calc_edge_length(eh); |
|
|
m_averageEdgeLen += edgeLen; |
|
|
} |
|
|
m_averageEdgeLen /= m_meshEdgeNum; |
|
|
|
|
|
|
|
|
m_sampleRadius = sampleRadiusRatio * m_averageEdgeLen; |
|
|
|
|
|
|
|
|
std::vector<size_t> vertexReorderedAlongAxis(m_meshVertexNum); |
|
|
size_t vertexIdx = 0; |
|
|
std::generate(vertexReorderedAlongAxis.begin(), vertexReorderedAlongAxis.end(), [&vertexIdx]() -> size_t { return vertexIdx++; }); |
|
|
std::sort(vertexReorderedAlongAxis.begin(), vertexReorderedAlongAxis.end(), [&mesh, axis](const size_t &lhs, const size_t &rhs) -> bool { |
|
|
size_t lhsIdx = lhs; |
|
|
size_t rhsIdx = rhs; |
|
|
OpenMesh::VertexHandle vhl = mesh.vertex_handle(lhsIdx); |
|
|
OpenMesh::VertexHandle vhr = mesh.vertex_handle(rhsIdx); |
|
|
Mesh::Point vl = mesh.point(vhl); |
|
|
Mesh::Point vr = mesh.point(vhr); |
|
|
return vl[axis] > vr[axis]; |
|
|
}); |
|
|
|
|
|
|
|
|
size_t firstVertexIdx = vertexReorderedAlongAxis[0]; |
|
|
VertexNodeIdx.resize(m_meshVertexNum); |
|
|
VertexNodeIdx.setConstant(-1); |
|
|
VertexNodeIdx[firstVertexIdx] = 0; |
|
|
size_t cur_node_idx = 0; |
|
|
|
|
|
m_vertexGraph.resize(m_meshVertexNum); |
|
|
VectorX weight_sum = VectorX::Zero(m_meshVertexNum); |
|
|
|
|
|
for (auto &vertexIdx : vertexReorderedAlongAxis) |
|
|
{ |
|
|
if(VertexNodeIdx[vertexIdx] < 0 && m_vertexGraph.at(vertexIdx).empty()) |
|
|
{ |
|
|
m_nodeContainer.emplace_back(cur_node_idx, vertexIdx); |
|
|
VertexNodeIdx[vertexIdx] = cur_node_idx; |
|
|
|
|
|
std::vector<size_t> neighbor_verts; |
|
|
geodesic::GeodesicAlgorithmExact geoalg(&mesh, vertexIdx, m_sampleRadius); |
|
|
geoalg.propagate(vertexIdx, neighbor_verts); |
|
|
for(size_t i = 0; i < neighbor_verts.size(); i++) |
|
|
{ |
|
|
int neighIdx = neighbor_verts[i]; |
|
|
Scalar geodist = mesh.data(mesh.vertex_handle(neighIdx)).geodesic_distance; |
|
|
if(geodist < m_sampleRadius) |
|
|
{ |
|
|
Scalar weight = std::pow(1-std::pow(geodist/m_sampleRadius, 2), 3); |
|
|
m_vertexGraph.at(neighIdx).emplace(std::pair<int, Scalar>(cur_node_idx, weight)); |
|
|
weight_sum[neighIdx] += weight; |
|
|
} |
|
|
} |
|
|
cur_node_idx++; |
|
|
} |
|
|
} |
|
|
|
|
|
m_nodeGraph.resize(cur_node_idx); |
|
|
for (auto &vertexIdx : vertexReorderedAlongAxis) |
|
|
{ |
|
|
for(auto &node: m_vertexGraph[vertexIdx]) |
|
|
{ |
|
|
size_t nodeIdx = node.first; |
|
|
for(auto &neighNode: m_vertexGraph[vertexIdx]) |
|
|
{ |
|
|
size_t neighNodeIdx = neighNode.first; |
|
|
if(nodeIdx != neighNodeIdx) |
|
|
{ |
|
|
m_nodeGraph.at(nodeIdx).emplace(std::pair<int, Scalar>(neighNodeIdx, 1.0)); |
|
|
} |
|
|
} |
|
|
m_vertexGraph.at(vertexIdx).at(nodeIdx) /= weight_sum[vertexIdx]; |
|
|
} |
|
|
} |
|
|
return m_sampleRadius; |
|
|
} |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
Scalar nodeSampler::SampleAndConstuct(Mesh &mesh, Scalar sampleRadiusRatio, const Matrix3X & src_points) |
|
|
{ |
|
|
|
|
|
m_meshVertexNum = mesh.n_vertices(); |
|
|
m_meshEdgeNum = mesh.n_edges(); |
|
|
m_mesh = & mesh; |
|
|
|
|
|
|
|
|
|
|
|
for (size_t i = 0; i < m_meshEdgeNum; ++i) |
|
|
{ |
|
|
OpenMesh::EdgeHandle eh = mesh.edge_handle(i); |
|
|
Scalar edgeLen = mesh.calc_edge_length(eh); |
|
|
m_averageEdgeLen += edgeLen; |
|
|
} |
|
|
m_averageEdgeLen /= m_meshEdgeNum; |
|
|
|
|
|
|
|
|
m_sampleRadius = sampleRadiusRatio * m_averageEdgeLen; |
|
|
|
|
|
|
|
|
Vector3 means = src_points.rowwise().mean(); |
|
|
Matrix3X demean_src_points = src_points.colwise() - means; |
|
|
Matrix33 covariance = demean_src_points * demean_src_points.transpose(); |
|
|
|
|
|
Eigen::SelfAdjointEigenSolver<Matrix33> eigensolver(covariance); |
|
|
Vector3 eigen_values = eigensolver.eigenvalues(); |
|
|
int max_idx = 2; |
|
|
if(eigen_values[0] > eigen_values[max_idx]) max_idx = 0; |
|
|
if(eigen_values[1] > eigen_values[max_idx]) max_idx = 1; |
|
|
|
|
|
Vector3 main_axis = eigensolver.eigenvectors().col(max_idx); |
|
|
|
|
|
VectorX projection = demean_src_points.transpose() * main_axis; |
|
|
|
|
|
projection_indices.resize(projection.size()); |
|
|
projection_indices.setLinSpaced(projection.size(), 0, projection.size() - 1); |
|
|
|
|
|
auto compare = [&projection](int i, int j) { |
|
|
return projection(i) < projection(j); |
|
|
}; |
|
|
|
|
|
std::sort(projection_indices.data(), projection_indices.data() + projection_indices.size(), compare); |
|
|
|
|
|
|
|
|
|
|
|
size_t firstVertexIdx = projection_indices[0]; |
|
|
VertexNodeIdx.resize(m_meshVertexNum); |
|
|
VertexNodeIdx.setConstant(-1); |
|
|
VertexNodeIdx[firstVertexIdx] = 0; |
|
|
size_t cur_node_idx = 0; |
|
|
|
|
|
m_vertexGraph.resize(m_meshVertexNum); |
|
|
VectorX weight_sum = VectorX::Zero(m_meshVertexNum); |
|
|
|
|
|
for(int idx = 0; idx < projection_indices.size(); idx++) |
|
|
{ |
|
|
int vertexIdx = projection_indices[idx]; |
|
|
if(VertexNodeIdx[vertexIdx] < 0 && m_vertexGraph.at(vertexIdx).empty()) |
|
|
{ |
|
|
m_nodeContainer.emplace_back(cur_node_idx, vertexIdx); |
|
|
VertexNodeIdx[vertexIdx] = cur_node_idx; |
|
|
|
|
|
std::vector<size_t> neighbor_verts; |
|
|
geodesic::GeodesicAlgorithmExact geoalg(&mesh, vertexIdx, m_sampleRadius); |
|
|
geoalg.propagate(vertexIdx, neighbor_verts); |
|
|
for(size_t i = 0; i < neighbor_verts.size(); i++) |
|
|
{ |
|
|
int neighIdx = neighbor_verts[i]; |
|
|
Scalar geodist = mesh.data(mesh.vertex_handle(neighIdx)).geodesic_distance; |
|
|
if(geodist < m_sampleRadius) |
|
|
{ |
|
|
Scalar weight = std::pow(1-std::pow(geodist/m_sampleRadius, 2), 3); |
|
|
m_vertexGraph.at(neighIdx).emplace(std::pair<int, Scalar>(cur_node_idx, weight)); |
|
|
weight_sum[neighIdx] += weight; |
|
|
} |
|
|
} |
|
|
cur_node_idx++; |
|
|
} |
|
|
} |
|
|
|
|
|
m_nodeGraph.resize(cur_node_idx); |
|
|
for(int idx = 0; idx < projection_indices.size(); idx++) |
|
|
{ |
|
|
int vertexIdx = projection_indices[idx]; |
|
|
for(auto &node: m_vertexGraph[vertexIdx]) |
|
|
{ |
|
|
size_t nodeIdx = node.first; |
|
|
for(auto &neighNode: m_vertexGraph[vertexIdx]) |
|
|
{ |
|
|
size_t neighNodeIdx = neighNode.first; |
|
|
if(nodeIdx != neighNodeIdx) |
|
|
{ |
|
|
m_nodeGraph.at(nodeIdx).emplace(std::pair<int, Scalar>(neighNodeIdx, 1.0)); |
|
|
} |
|
|
} |
|
|
m_vertexGraph.at(vertexIdx).at(nodeIdx) /= weight_sum[vertexIdx]; |
|
|
} |
|
|
} |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
return m_sampleRadius; |
|
|
} |
|
|
|
|
|
|
|
|
Scalar nodeSampler::SampleAndConstuctForSrcPoints(Mesh &mesh, Scalar sampleRadiusRatio, const Matrix3X & src_points, const Eigen::MatrixXi & src_knn_indices) |
|
|
{ |
|
|
|
|
|
m_meshVertexNum = src_points.cols(); |
|
|
int knn_num_neighbor = src_knn_indices.rows(); |
|
|
m_meshEdgeNum = m_meshVertexNum*knn_num_neighbor; |
|
|
m_mesh = & mesh; |
|
|
|
|
|
|
|
|
for (size_t i = 0; i < m_meshVertexNum; ++i) |
|
|
{ |
|
|
for(size_t j = 0; j < knn_num_neighbor; j++) |
|
|
{ |
|
|
Scalar edgeLen = (src_points.col(i) - src_points.col(src_knn_indices(j, i))).norm(); |
|
|
m_averageEdgeLen += edgeLen; |
|
|
} |
|
|
} |
|
|
m_averageEdgeLen /= m_meshEdgeNum; |
|
|
|
|
|
|
|
|
|
|
|
m_sampleRadius = sampleRadiusRatio * m_averageEdgeLen; |
|
|
|
|
|
|
|
|
Vector3 means = src_points.rowwise().mean(); |
|
|
Matrix3X demean_src_points = src_points.colwise() - means; |
|
|
Matrix33 covariance = demean_src_points * demean_src_points.transpose(); |
|
|
|
|
|
Eigen::SelfAdjointEigenSolver<Matrix33> eigensolver(covariance); |
|
|
Vector3 eigen_values = eigensolver.eigenvalues(); |
|
|
int max_idx = 2; |
|
|
if(eigen_values[0] > eigen_values[max_idx]) max_idx = 0; |
|
|
if(eigen_values[1] > eigen_values[max_idx]) max_idx = 1; |
|
|
|
|
|
Vector3 main_axis = eigensolver.eigenvectors().col(max_idx); |
|
|
|
|
|
VectorX projection = demean_src_points.transpose() * main_axis; |
|
|
|
|
|
projection_indices.resize(projection.size()); |
|
|
projection_indices.setLinSpaced(projection.size(), 0, projection.size() - 1); |
|
|
|
|
|
auto compare = [&projection](int i, int j) { |
|
|
return projection(i) < projection(j); |
|
|
}; |
|
|
|
|
|
std::sort(projection_indices.data(), projection_indices.data() + projection_indices.size(), compare); |
|
|
|
|
|
|
|
|
|
|
|
size_t firstVertexIdx = projection_indices[0]; |
|
|
VertexNodeIdx.resize(m_meshVertexNum); |
|
|
VertexNodeIdx.setConstant(-1); |
|
|
VertexNodeIdx[firstVertexIdx] = 0; |
|
|
size_t cur_node_idx = 0; |
|
|
|
|
|
m_vertexGraph.resize(m_meshVertexNum); |
|
|
VectorX weight_sum = VectorX::Zero(m_meshVertexNum); |
|
|
|
|
|
for(int idx = 0; idx < projection_indices.size(); idx++) |
|
|
{ |
|
|
int vertexIdx = projection_indices[idx]; |
|
|
if(VertexNodeIdx[vertexIdx] < 0 && m_vertexGraph.at(vertexIdx).empty()) |
|
|
{ |
|
|
m_nodeContainer.emplace_back(cur_node_idx, vertexIdx); |
|
|
VertexNodeIdx[vertexIdx] = cur_node_idx; |
|
|
|
|
|
std::queue<size_t> neighbor_verts; |
|
|
neighbor_verts.push(vertexIdx); |
|
|
Eigen::VectorXi visited(m_meshVertexNum); |
|
|
visited.setZero(); |
|
|
visited[vertexIdx] = 1; |
|
|
while(!neighbor_verts.empty()) |
|
|
{ |
|
|
size_t vidx = neighbor_verts.front(); |
|
|
neighbor_verts.pop(); |
|
|
Scalar dist = (src_points.col(vertexIdx) - src_points.col(vidx)).norm(); |
|
|
|
|
|
if(dist < m_sampleRadius) |
|
|
{ |
|
|
Scalar weight = std::pow(1-std::pow(dist/m_sampleRadius, 2), 3); |
|
|
m_vertexGraph.at(vidx).emplace(std::pair<int, Scalar>(cur_node_idx, weight)); |
|
|
weight_sum[vidx] += weight; |
|
|
for(int j = 0; j < knn_num_neighbor; j++) |
|
|
{ |
|
|
int neighbor_vidx = src_knn_indices(j, vidx); |
|
|
if(visited[neighbor_vidx]==0) |
|
|
{ |
|
|
neighbor_verts.push(neighbor_vidx); |
|
|
visited[neighbor_vidx]= 1; |
|
|
} |
|
|
} |
|
|
} |
|
|
} |
|
|
|
|
|
cur_node_idx++; |
|
|
} |
|
|
} |
|
|
|
|
|
m_nodeGraph.resize(cur_node_idx); |
|
|
for(int idx = 0; idx < projection_indices.size(); idx++) |
|
|
{ |
|
|
int vertexIdx = projection_indices[idx]; |
|
|
for(auto &node: m_vertexGraph[vertexIdx]) |
|
|
{ |
|
|
size_t nodeIdx = node.first; |
|
|
for(auto &neighNode: m_vertexGraph[vertexIdx]) |
|
|
{ |
|
|
size_t neighNodeIdx = neighNode.first; |
|
|
if(nodeIdx != neighNodeIdx) |
|
|
{ |
|
|
m_nodeGraph.at(nodeIdx).emplace(std::pair<int, Scalar>(neighNodeIdx, 1.0)); |
|
|
} |
|
|
} |
|
|
m_vertexGraph.at(vertexIdx).at(nodeIdx) /= weight_sum[vertexIdx]; |
|
|
} |
|
|
} |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
return m_sampleRadius; |
|
|
} |
|
|
|
|
|
|
|
|
Scalar nodeSampler::SampleAndConstuctFPS(Mesh &mesh, Scalar sampleRadiusRatio, const Matrix3X & src_points, const Eigen::MatrixXi & src_knn_indices, int num_vn, int num_nn) |
|
|
{ |
|
|
m_meshVertexNum = src_points.cols(); |
|
|
int knn_num_neighbor = src_knn_indices.rows(); |
|
|
m_meshEdgeNum = m_meshVertexNum*knn_num_neighbor; |
|
|
m_mesh = & mesh; |
|
|
|
|
|
|
|
|
for (size_t i = 0; i < m_meshVertexNum; ++i) |
|
|
{ |
|
|
for(size_t j = 0; j < knn_num_neighbor; j++) |
|
|
{ |
|
|
Scalar edgeLen = (src_points.col(i) - src_points.col(src_knn_indices(j, i))).norm(); |
|
|
m_averageEdgeLen += edgeLen; |
|
|
} |
|
|
} |
|
|
m_averageEdgeLen /= m_meshEdgeNum; |
|
|
|
|
|
|
|
|
m_sampleRadius = sampleRadiusRatio * m_averageEdgeLen; |
|
|
|
|
|
|
|
|
size_t startIndex = 0; |
|
|
|
|
|
VectorX minDistances(m_meshVertexNum); |
|
|
minDistances.setConstant(std::numeric_limits<Scalar>::max()); |
|
|
minDistances[startIndex] = 0; |
|
|
m_nodeContainer.emplace_back(startIndex, 0); |
|
|
|
|
|
Scalar minimal_dist = 1e10; |
|
|
int cur_node_idx = 1; |
|
|
|
|
|
while (minimal_dist > m_sampleRadius) { |
|
|
|
|
|
#pragma omp parallel for |
|
|
for (size_t i = 0; i < m_meshVertexNum; ++i) { |
|
|
if(i==startIndex) |
|
|
continue; |
|
|
|
|
|
Scalar dist = (src_points.col(startIndex) - src_points.col(i)).norm(); |
|
|
if(dist < minDistances[i]) |
|
|
minDistances[i] = dist; |
|
|
} |
|
|
|
|
|
|
|
|
int maxDistanceIndex; |
|
|
minimal_dist = minDistances.maxCoeff(&maxDistanceIndex); |
|
|
minDistances[maxDistanceIndex] = 0; |
|
|
|
|
|
|
|
|
startIndex = maxDistanceIndex; |
|
|
m_nodeContainer.emplace_back(cur_node_idx, startIndex); |
|
|
cur_node_idx++; |
|
|
} |
|
|
|
|
|
Matrix3X node_positions(3, cur_node_idx); |
|
|
for(int i = 0; i < cur_node_idx; i++) |
|
|
{ |
|
|
int vidx = getNodeVertexIdx(i); |
|
|
node_positions.col(i) = src_points.col(vidx); |
|
|
} |
|
|
KDtree node_tree(node_positions); |
|
|
|
|
|
|
|
|
m_vertexGraph.resize(m_meshVertexNum); |
|
|
VectorX weight_sum(m_meshVertexNum); |
|
|
weight_sum.setZero(); |
|
|
#pragma omp parallel for |
|
|
for(int vidx = 0; vidx < m_meshVertexNum; vidx++) |
|
|
{ |
|
|
std::vector<int> out_indices(num_vn); |
|
|
std::vector<Scalar> out_dists(num_vn); |
|
|
node_tree.query(src_points.col(vidx).data(), num_vn, out_indices.data(), out_dists.data()); |
|
|
for(int k = 0; k < num_vn; k++) |
|
|
{ |
|
|
Scalar weight = std::pow(1-std::pow(out_dists[k]/m_sampleRadius, 2), 3); |
|
|
m_vertexGraph.at(vidx).emplace(std::pair<int, Scalar>(out_indices[k], weight)); |
|
|
weight_sum[vidx] += weight; |
|
|
} |
|
|
} |
|
|
|
|
|
#pragma omp parallel for |
|
|
for(int vidx = 0; vidx < m_meshVertexNum; vidx++) |
|
|
{ |
|
|
for(auto &neighNode: m_vertexGraph[vidx]) |
|
|
{ |
|
|
size_t neighNodeIdx = neighNode.first; |
|
|
m_vertexGraph.at(vidx).at(neighNodeIdx) /= weight_sum[vidx]; |
|
|
} |
|
|
} |
|
|
|
|
|
|
|
|
|
|
|
m_nodeGraph.resize(cur_node_idx); |
|
|
#pragma omp parallel for |
|
|
for(int nidx = 0; nidx < cur_node_idx; nidx++) |
|
|
{ |
|
|
std::vector<int> out_indices(num_nn+1); |
|
|
std::vector<Scalar> out_dists(num_nn+1); |
|
|
int vidx = getNodeVertexIdx(nidx); |
|
|
node_tree.query(src_points.col(vidx).data(), num_nn+1, out_indices.data(), out_dists.data()); |
|
|
for(int k = 0; k < num_nn; k++) |
|
|
{ |
|
|
m_nodeGraph.at(nidx).emplace(std::pair<int, Scalar>(out_indices[k+1], 1.0)); |
|
|
} |
|
|
} |
|
|
return m_sampleRadius; |
|
|
} |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
void nodeSampler::initWeight(RowMajorSparseMatrix& matPV, VectorX & matP, RowMajorSparseMatrix& matB, VectorX& matD, VectorX& smoothw) |
|
|
{ |
|
|
Timer time; |
|
|
std::vector<Eigen::Triplet<Scalar>> coeff; |
|
|
matP.setZero(); |
|
|
Eigen::VectorXi nonzero_num = Eigen::VectorXi::Zero(m_mesh->n_vertices()); |
|
|
|
|
|
for (size_t vertexIdx = 0; vertexIdx < m_meshVertexNum; ++vertexIdx) |
|
|
{ |
|
|
Mesh::Point vi = m_mesh->point(m_mesh->vertex_handle(vertexIdx)); |
|
|
for (auto &eachNeighbor : m_vertexGraph[vertexIdx]) |
|
|
{ |
|
|
size_t nodeIdx = eachNeighbor.first; |
|
|
Scalar weight = m_vertexGraph.at(vertexIdx).at(nodeIdx); |
|
|
Mesh::Point pj = m_mesh->point(m_mesh->vertex_handle(getNodeVertexIdx(nodeIdx))); |
|
|
|
|
|
for (int k = 0; k < 3; k++) |
|
|
{ |
|
|
coeff.push_back(Eigen::Triplet<Scalar>(3 * vertexIdx + k, nodeIdx * 12 + k, weight * (vi[0] - pj[0]))); |
|
|
coeff.push_back(Eigen::Triplet<Scalar>(3 * vertexIdx + k, nodeIdx * 12 + k + 3, weight * (vi[1] - pj[1]))); |
|
|
coeff.push_back(Eigen::Triplet<Scalar>(3 * vertexIdx + k, nodeIdx * 12 + k + 6, weight * (vi[2] - pj[2]))); |
|
|
coeff.push_back(Eigen::Triplet<Scalar>(3 * vertexIdx + k, nodeIdx * 12 + k + 9, weight * 1.0)); |
|
|
} |
|
|
|
|
|
matP[vertexIdx * 3] += weight * pj[0]; |
|
|
matP[vertexIdx * 3 + 1] += weight * pj[1]; |
|
|
matP[vertexIdx * 3 + 2] += weight * pj[2]; |
|
|
} |
|
|
nonzero_num[vertexIdx] = m_vertexGraph[vertexIdx].size(); |
|
|
} |
|
|
matPV.setFromTriplets(coeff.begin(), coeff.end()); |
|
|
|
|
|
|
|
|
|
|
|
coeff.clear(); |
|
|
int max_edge_num = nodeSize() * (nodeSize()-1); |
|
|
matB.resize(max_edge_num * 3, 12 * nodeSize()); |
|
|
matD.resize(max_edge_num * 3); |
|
|
smoothw.resize(max_edge_num*3); |
|
|
smoothw.setZero(); |
|
|
matD.setZero(); |
|
|
int edge_id = 0; |
|
|
|
|
|
for (size_t nodeIdx = 0; nodeIdx < m_nodeContainer.size(); ++nodeIdx) |
|
|
{ |
|
|
size_t vIdx0 = getNodeVertexIdx(nodeIdx); |
|
|
Mesh::VertexHandle vh0 = m_mesh->vertex_handle(vIdx0); |
|
|
Mesh::Point v0 = m_mesh->point(vh0); |
|
|
|
|
|
for (auto &eachNeighbor : m_nodeGraph[nodeIdx]) |
|
|
{ |
|
|
size_t neighborIdx = eachNeighbor.first; |
|
|
size_t vIdx1 = getNodeVertexIdx(neighborIdx); |
|
|
Mesh::Point v1 = m_mesh->point(m_mesh->vertex_handle(vIdx1)); |
|
|
Mesh::Point dv = v0 - v1; |
|
|
int k = edge_id; |
|
|
|
|
|
for (int t = 0; t < 3; t++) |
|
|
{ |
|
|
coeff.push_back(Eigen::Triplet<Scalar>(k * 3 + t, neighborIdx * 12 + t, dv[0])); |
|
|
coeff.push_back(Eigen::Triplet<Scalar>(k * 3 + t, neighborIdx * 12 + t + 3, dv[1])); |
|
|
coeff.push_back(Eigen::Triplet<Scalar>(k * 3 + t, neighborIdx * 12 + t + 6, dv[2])); |
|
|
coeff.push_back(Eigen::Triplet<Scalar>(k * 3 + t, neighborIdx * 12 + t + 9, 1.0)); |
|
|
coeff.push_back(Eigen::Triplet<Scalar>(k * 3 + t, nodeIdx * 12 + t + 9, -1.0)); |
|
|
} |
|
|
|
|
|
Scalar dist = dv.norm(); |
|
|
if(dist > 0) |
|
|
{ |
|
|
smoothw[k * 3] = smoothw[k * 3 + 1] = smoothw[k * 3 + 2] = 1.0 / dist; |
|
|
} |
|
|
else |
|
|
{ |
|
|
|
|
|
std::cout << "node repeat"; |
|
|
exit(1); |
|
|
} |
|
|
matD[k * 3] = dv[0]; |
|
|
matD[k * 3 + 1] = dv[1]; |
|
|
matD[k * 3 + 2] = dv[2]; |
|
|
edge_id++; |
|
|
} |
|
|
} |
|
|
matB.setFromTriplets(coeff.begin(), coeff.end()); |
|
|
matD.conservativeResize(edge_id*3); |
|
|
matB.conservativeResize(edge_id*3, matPV.cols()); |
|
|
smoothw.conservativeResize(edge_id*3); |
|
|
smoothw *= edge_id/(smoothw.sum()/3.0); |
|
|
} |
|
|
|
|
|
|
|
|
void nodeSampler::print_nodes(Mesh & mesh, std::string file_path) |
|
|
{ |
|
|
std::string namev = file_path + "nodes.obj"; |
|
|
std::ofstream out1(namev); |
|
|
for (size_t i = 0; i < m_nodeContainer.size(); i++) |
|
|
{ |
|
|
int vexid = m_nodeContainer[i].second; |
|
|
out1 << "v " << mesh.point(mesh.vertex_handle(vexid))[0] << " " << mesh.point(mesh.vertex_handle(vexid))[1] |
|
|
<< " " << mesh.point(mesh.vertex_handle(vexid))[2] << std::endl; |
|
|
} |
|
|
Eigen::VectorXi nonzero_num = Eigen::VectorXi::Zero(m_nodeContainer.size()); |
|
|
for (size_t nodeIdx = 0; nodeIdx < m_nodeContainer.size(); ++nodeIdx) |
|
|
{ |
|
|
for (auto &eachNeighbor : m_nodeGraph[nodeIdx]) |
|
|
{ |
|
|
size_t neighborIdx = eachNeighbor.first; |
|
|
out1 << "l " << nodeIdx+1 << " " << neighborIdx+1 << std::endl; |
|
|
} |
|
|
nonzero_num[nodeIdx] = m_nodeGraph[nodeIdx].size(); |
|
|
} |
|
|
std::cout << "node neighbor min = " << nonzero_num.minCoeff() << " max = " |
|
|
<< nonzero_num.maxCoeff() << " average = " << nonzero_num.mean() << std::endl; |
|
|
out1.close(); |
|
|
} |
|
|
} |
|
|
|