[mod]: quadric smooth

.gitmodules

@@ -1,6 +1,6 @@
-[submodule "external/probabilistic-quadrics"]
- path = external/probabilistic-quadrics
- url = [email protected]:anapupa/probabilistic-quadrics.git
 [submodule "external/happly"]
  path = external/happly
  url = [email protected]:nmwsharp/happly.git
+[submodule "external/probabilistic-quadrics"]
+ path = external/probabilistic-quadrics
+ url = https://github.com/Philip-Trettner/probabilistic-quadrics

CMakeLists.txt

@@ -84,11 +84,11 @@ FILE(GLOB CU_UTILS_SRC
  cuMesh/**/*.cpp)
 
 set(CUDA_LINK_LIBRARIES_KEYWORD PRIVATE)
+set( CUDA_NVCC_FLAGS "${CUDA_NVCC_FLAGS}" "--extended-lambda" )
 
 include_directories(${PROJECT_SOURCE_DIR})
-
 cuda_add_library(cuMesh_Impl ${CU_UTILS_SRC_IMPL} )
-TARGET_LINK_LIBRARIES(cuMesh_Impl
+target_link_libraries(cuMesh_Impl
  PRIVATE ${CUDA_LIBRARIES}
  PRIVATE ${CUDA_CUDART_LIBRARY}
  PRIVATE ${CUDA_cublas_LIBRARY}

cuMesh/HalfEdgeNavigators.cuh

@@ -93,7 +93,7 @@ namespace cuMesh {
 
 
  // Geometry Query
- __host__ __device__ __forceinline__ float3 triangle_centroid(Index fid) {
+ __host__ __device__ __forceinline__ float3 triangle_centroid(Index fid) const {
  return (_vert[_face[fid].x]+_vert[_face[fid].y]+_vert[_face[fid].z])/3.0f;
  }
 

cuMesh/algorithm/SurfaceMeshOperation.h

@@ -4,13 +4,17 @@
 #pragma once
 #include <cuMesh/HalfEdgeNavigators.cuh>
 
+#define Q_API __device__ __host__
+#include <probabilistic-quadrics.hh>
+#include <minimal-math.hh>
 // Helper struct for thrust
 
-struct Quadric;
-
 namespace cuMesh {
 
+using Quadric = pq::quadric< pq::math<double, pq::vec<3, double>, pq::vec<3, double>, pq::mat<3, 3, double>> >;
+float3 Q_API from_pq(const pq::dvec3& v) ;
 
+pq::dvec3 Q_API to_pq(const float3& v) ;
 /**
  * @brief function struct for thrust
  */
@@ -26,11 +30,16 @@ struct ComputeFaceNormal: MeshNavigators {
  __device__ float3 operator()(Index fid) ;
 };
 
+struct ComputeVertexNormal: MeshNavigators {
+ using MeshNavigators::MeshNavigators;
+ __device__ float3 operator()(Index vid) ;
+};
+
 struct ComputeFaceQuadric: MeshNavigators {
  using MeshNavigators::MeshNavigators;
  float _stddev{0};
  __host__ ComputeFaceQuadric(VFMeshData& meshData, float stddev): MeshNavigators(meshData), _stddev(stddev) {}
- __device__ Quadric operator()(Index& fid,float3& normal) ;
+ __device__ Quadric operator()(Index& fid) ;
 };
 
 struct CheckTriangleIsDeleted: MeshNavigators {

cuMesh/algorithm/SurfaceSmoothing.h

@@ -3,11 +3,15 @@
 //
 #pragma once
 #include <cuMesh/HalfEdgeElementType.h>
-#include <probabilistic-quadrics.cuh>
+
+#define Q_API __device__ __host__
+#include <probabilistic-quadrics.hh>
+#include <minimal-math.hh>
 
 namespace cuMesh {
 
 class MeshNavigators;
+using Quadric = pq::quadric< pq::math<double, pq::vec<3, double>, pq::vec<3, double>, pq::mat<3, 3, double>> >;
 
 /**
  * @brief Mesh Smoothing/Denoising
@@ -31,6 +35,7 @@ void PQGFSmoothing(VFMeshData& mesh_data, double sigma_s, double sigma_r, double
 namespace Kernel{
  using cuMesh::MeshNavigators;
  using cuMesh::Index;
+ using cuMesh::Quadric;
  /**
  * Facet's Quadric Diffusion
  * @param navigators
@@ -40,10 +45,12 @@ namespace Kernel{
  * @param sigma_r
  */
  __global__ void DiffuseTriQuadric(MeshNavigators navigators, Quadric* tri_quadric, Quadric* tri_quadric_out,
- float sigma_s, float sigma_r);
+ double sigma_s, double  sigma_r);
 
  __global__ void UpdateVertexPosition(MeshNavigators navigators, Quadric* tri_quadric, double sigma_v);
 
+ __global__ void UpdateVertexQuadric(MeshNavigators navi, Quadric* tri_quadric, Quadric* vert_quadric, double sigma_v_sq);
+ __global__ void UpdateVertexPositionConstraint(MeshNavigators navi, Quadric* vert_quadric, float3* f_n);
 }
 
 

cuMesh/cuda/SurfaceMeshOperation.cu

@@ -2,8 +2,17 @@
 // Created by pupa on 2021/11/13.
 //
 #include <cuMesh/HalfEdgeElementType.h>
+
 #include "cuMesh/algorithm/SurfaceMeshOperation.h"
-#include <probabilistic-quadrics.cuh>
+#include "cuMesh/algorithm/SurfaceSmoothing.h"
+
+float3 Q_API cuMesh::from_pq(const pq::dvec3 & v) {
+ return make_float3(v.x, v.y, v.z);
+}
+
+pq::dvec3 Q_API cuMesh::to_pq(const float3& v) {
+ return pq::dvec3 {v.x, v.y, v.z};
+}
 
 namespace cuMesh::Operation{
 
@@ -19,9 +28,22 @@ __device__ float3 ComputeFaceNormal::operator()(Index fid) {
  return normalize(cross(v1-v0, v2-v0));
 }
 
-__device__ Quadric ComputeFaceQuadric::operator()(Index& fid, float3& normal) {
+__device__ float3 ComputeVertexNormal::operator()(Index vid) {
+ float3 normal_sum{0 ,0,0};
+ for(auto it = outgoing_hedges(vid); !it.is_end(); ++it) {
+ Index ffid = hedge2face(*it);
+ normal_sum += this->cross_dot(ffid);
+ }
+}
+
+__device__ Quadric ComputeFaceQuadric::operator()(Index& fid) {
  float3 v0 = position(vertex(fid, 0));
- return Quadric::PlaneQuadric(v0, normal, _stddev, 0.1);
+ float3 v1 = position(vertex(fid, 1));
+ float3 v2 = position(vertex(fid, 2));
+
+// return cuMesh::Quadric::plane_quadric(to_pq(v0), to_pq(normal));
+ return Quadric::probabilistic_plane_quadric(to_pq((v0+v1+v2)/3.0f), to_pq(normal(fid)), _stddev, 0.1);
+// return Quadric::probabilistic_triangle_quadric(to_pq(v0), to_pq(v1), to_pq(v2), _stddev);
 }
 
 __host__ __device__ bool CheckHedgeIsBoundary::operator()(Index hid) {

cuMesh/cuda/SurfaceSmoothing.cu

@@ -18,7 +18,7 @@ void Smoothing::PQGFSmoothing(VFMeshData &mesh_data, double sigma_s, double sigm
  MeshNavigators navigators(mesh_data);
  thrust::device_vector<Index> face_index(navigators.n_f);
  thrust::sequence(face_index.begin(), face_index.end());
-
+ make_double3(1, 2, 3);
  float mean_e_len = 0;
  {
  thrust::device_vector<float> hedge_length(navigators.n_f);
@@ -29,28 +29,35 @@ void Smoothing::PQGFSmoothing(VFMeshData &mesh_data, double sigma_s, double sigm
  }
 
  sigma_r *= mean_e_len;
- sigma_s *= mean_e_len;
  sigma_v *= mean_e_len;
 
-
  //step1. compute triangle quadric
- thrust::device_vector<float3> tri_normal(navigators.n_f);
- thrust::transform(face_index.begin(), face_index.end(), tri_normal.begin(),
- Operation::ComputeFaceNormal(mesh_data));
  thrust::device_vector<Quadric> tri_quadric(navigators.n_f);
- thrust::transform(face_index.begin(), face_index.end(), tri_normal.begin(), tri_quadric.begin(),
+ thrust::transform(face_index.begin(), face_index.end(), tri_quadric.begin(),
  Operation::ComputeFaceQuadric(mesh_data, mean_e_len*0.05));
 
 
  //step2. tri_quadric diffusion
  thrust::device_vector<Quadric> tri_quadric2(navigators.n_f);
- for(int i = 0; i < 3; i++) {
+ for(int i = 0; i < 5; i++) {
  Kernel::DiffuseTriQuadric<<<GetGridDim(navigators.n_f), blk_size, 0>>>(
  navigators, tri_quadric.data().get(), tri_quadric2.data().get(),
  sigma_s * sigma_s, sigma_r*sigma_r);
  thrust::swap(tri_quadric, tri_quadric2);
  }
 
+ thrust::device_vector<Quadric> vert_quadric(navigators.n_v);
+ Kernel::UpdateVertexQuadric<<<GetGridDim(navigators.n_v), blk_size, 0>>>(
+ navigators, tri_quadric.data().get(), vert_quadric.data().get(), sigma_v * sigma_v);
+
+// thrust::device_vector<float3> tri_nromal(navigators.n_f);
+// for(int i = 0; i < 10; i++) {
+// thrust::transform(face_index.begin(), face_index.end(), tri_nromal.begin(),
+// Operation::ComputeFaceNormal(mesh_data));
+// Kernel::UpdateVertexPositionConstraint<<<GetGridDim(navigators.n_v), blk_size, 0>>>(
+// navigators, vert_quadric.data().get(), tri_nromal.data().get() );
+// }
+
  Kernel::UpdateVertexPosition<<<GetGridDim(navigators.n_v), blk_size, 0>>>(
  navigators, tri_quadric.data().get(), sigma_v * sigma_v);
 }
@@ -60,26 +67,40 @@ void Smoothing::PQGFSmoothing(VFMeshData &mesh_data, double sigma_s, double sigm
 
 
 __global__ void Kernel::DiffuseTriQuadric(MeshNavigators navi, Quadric* tri_quadric, Quadric* tri_quadric_out,
- float sigma_s_sq, float sigma_r_sq) {
+ double sigma_s_sq, double sigma_r_sq) {
  uint32_t fid = blockIdx.x * blockDim.x + threadIdx.x;
- if(fid > navi.n_f ) return ;
- Quadric Qf_new, Qf = tri_quadric[fid];
- float ww = 0, w, d_s, d_r;
+ if(fid >= navi.n_f ) return ;
+ Quadric Qf_new, Qf = tri_quadric[fid], Qff;
+ double ww = 0, w, d_s, d_r, w_s, w_r, w_angle;
+
+ double qf_norm = pq::minimal_math<double>::frobenius_inner_product(Qf.A()), qff_norm;
  float3 fid_center = navi.triangle_centroid(fid);
 #pragma unroll 3
  for(int i = 0; i < 3; i++) {
  for(auto it = navi.outgoing_hedges(navi.vertex(fid, i)); !it.is_end(); ++it) {
  Index ffid = navi.hedge2face(*it);
  if( ffid == NonIndex) continue;
+ Qff = tri_quadric[ffid] ;
+ qff_norm = pq::minimal_math<double>::frobenius_inner_product(Qff.A());
+
  float3 _v0 = navi.position( navi.vertex(ffid, 0) );
  float3 _v1 = navi.position( navi.vertex(ffid, 1) );
  float3 _v2 = navi.position( navi.vertex(ffid, 2) );
 
  d_s = length( fid_center - (_v0+_v1+_v2)/3.0f );
- d_r = (Qf(_v0)+Qf(_v1)+Qf(_v2)) * length( navi.cross_dot(ffid) ) / 6.0f;
+ d_s = ::exp(-0.5 * d_s *d_s / sigma_s_sq )* length(navi.cross_dot(ffid));
 
- w = ::exp(-0.5 * d_s *d_s / sigma_s_sq ) * ::exp(-0.5 * d_r / sigma_r_sq );
- Qf_new += tri_quadric[ffid] * w;
+ d_r = (Qf(_v0.x, _v0.y, _v0.z)+Qf(_v1.x, _v1.y, _v1.z)+Qf(_v2.x, _v2.y, _v2.z))
+ * length( navi.cross_dot(ffid) ) / 6.0f;
+ d_r = exp(-0.5*d_r/sigma_r_sq);
+ w_angle = pq::minimal_math<double>::trace_of_product(Qff.A(), Qf.A())/sqrt(qf_norm*qff_norm);
+
+//// printf("[face %d]: |A|: %lf, |B|: %lf > <A,B>/|A||B|: %lf\n", fid, qf_norm, qff_norm, d_r);
+//
+ w_angle = (w_angle > 0.11)?w_angle:0; //
+// d_r = w_angle = 1;
+ w = max(d_r * d_s * w_angle, 0.0001);
+ Qf_new += Qff * w;
  ww += w;
  }
  }
@@ -92,20 +113,24 @@ __global__ void Kernel::DiffuseTriQuadric(MeshNavigators navi, Quadric* tri_quad
  for(auto it = navi.outgoing_hedges(vid); !it.is_end(); ++it) {
  Index ffid = navi.hedge2face(*it);
  if( ffid == NonIndex) continue;
+ Qff = tri_quadric[ffid] ;
+ qff_norm = pq::minimal_math<double>::frobenius_inner_product(Qff.A());
+
  float3 _v0 = navi.position( navi.vertex(ffid, 0) );
  float3 _v1 = navi.position( navi.vertex(ffid, 1) );
  float3 _v2 = navi.position( navi.vertex(ffid, 2) );
 
  d_s = length( fid_center - (_v0+_v1+_v2)/3.0f );
- d_r = (Qf(_v0)+Qf(_v1)+Qf(_v2)) * length( navi.cross_dot(ffid) ) / 6.0f;
+ d_s = ::exp(-0.5 * d_s *d_s / sigma_s_sq ) * length(navi.cross_dot(ffid));
 
  w = ::exp(-0.5 * d_s *d_s / sigma_s_sq ) * ::exp(-0.5 * d_r / sigma_r_sq );
  Qf_new += tri_quadric[ffid] * w;
  ww += w;
  }
  }
 
- tri_quadric_out[fid] = Qf_new / ww;
+ assert(abs(ww) > 0.0001 );
+ tri_quadric_out[fid] = Qf_new / ww;
 }
 
 
@@ -114,29 +139,71 @@ __global__ void Kernel::UpdateVertexPosition(MeshNavigators navi, Quadric* tri_q
  if(vid >= navi.n_v || navi.is_v_boundary(vid) ) return ;
  float3 p = navi.position(vid);
  Quadric Qv;
- float w, ww = 0;
+ double w, ww = 0;
  for(auto it = navi.outgoing_hedges(vid); !it.is_end(); ++it) {
  Index ffid = navi.hedge2face(*it);
  if( ffid == NonIndex ) continue;
  float3 _v0 = navi.position( navi.vertex(ffid, 0) );
  float3 _v1 = navi.position( navi.vertex(ffid, 1) );
  float3 _v2 = navi.position( navi.vertex(ffid, 2) );
 
- float3 diff = p - (_v0+_v1+_v2)/3.0f ;
- w = ::exp(-0.5 * dot(diff, diff) / sigma_v_sq ) ;
+ Quadric q = Quadric::point_quadric( cuMesh::to_pq((_v0+_v1+_v2)/3.0f  ) );
+ w = length( navi.cross_dot(ffid) );
 
- Qv += tri_quadric[ffid] * w;
+ Qv += (tri_quadric[ffid]) * w;
  ww += w;
- Qv += Quadric::PointQuadric(p) * w * 0.05f;
+
  }
- Qv /= ww;
+ Qv += Quadric::point_quadric(cuMesh::to_pq(p)) * 0.005;
 
- float3 min_p = Qv.minimizer();
+ auto min_p = Qv.minimizer();
 // printf("%f %f %f\n", min_p.x, min_p.y, min_p.z);
 
  if(min_p.x == min_p.x && min_p.y == min_p.y && min_p.z == min_p.z) {
 // printf("%f %f %f --> %f %f %f\n", p.x, p.y, p.z, min_p.x, min_p.y, min_p.z);
- navi._vert[vid] = min_p;
+ navi._vert[vid] = make_float3(min_p.x, min_p.y, min_p.z);
+ }
+ else printf("failed\n");
+}
+
+
+
+__global__ void Kernel::UpdateVertexQuadric(MeshNavigators navi, Quadric* tri_quadric, Quadric* vert_quadric, double sigma_v_sq) {
+ uint32_t vid = blockIdx.x * blockDim.x + threadIdx.x;
+ if(vid >= navi.n_v || navi.is_v_boundary(vid) ) return ;
+ float3 p = navi.position(vid);
+ Quadric Qv;
+ double w, ww = 0;
+ for(auto it = navi.outgoing_hedges(vid); !it.is_end(); ++it) {
+ Index ffid = navi.hedge2face(*it);
+ if( ffid == NonIndex ) continue;
+// float3 _v0 = navi.position( navi.vertex(ffid, 0) );
+// float3 _v1 = navi.position( navi.vertex(ffid, 1) );
+// float3 _v2 = navi.position( navi.vertex(ffid, 2) );
+//
+// Quadric q = Quadric::point_quadric( cuMesh::to_pq((_v0+_v1+_v2)/3.0f ) );
+ w = length( navi.cross_dot(ffid) );
+
+ Qv += (tri_quadric[ffid]) * w;
+ ww += w;
+
+ }
+ vert_quadric[vid] = Qv/ww + Quadric::point_quadric(cuMesh::to_pq(p)) * 0.005;
+}
+
+__global__ void Kernel::UpdateVertexPositionConstraint(MeshNavigators navi, Quadric* vert_quadric, float3* f_n){
+ uint32_t vid = blockIdx.x * blockDim.x + threadIdx.x;
+ if(vid >= navi.n_v || navi.is_v_boundary(vid) ) return ;
+
+ double ww = 0;
+ float3 v_motion{0 ,0 ,0};
+ for(auto it = navi.outgoing_hedges(vid); !it.is_end(); ++it) {
+ Index ffid = navi.hedge2face(*it);
+ if( ffid == NonIndex ) continue;
+ double tri_area = length(navi.cross_dot(ffid));
+ double lambda = vert_quadric[vid].constraint_minimizer(cuMesh::to_pq(navi.position(vid)), cuMesh::to_pq(f_n[ffid]));
+ v_motion += lambda * f_n[ffid] * tri_area;
+ ww += tri_area;
  }
-// else printf("failed\n");
+ navi._vert[vid] += v_motion/ww;
 }

tests/TestCudaConvertHedge.cpp

@@ -11,7 +11,9 @@ int main (int argc, char **argv) {
  cuMesh::readPLY(argv[1], mesh_data);
 
  cuMesh::Topology::UpdateHalfEdge(mesh_data);
- cuMesh::Smoothing::PQGFSmoothing(mesh_data, 2, 5, 1);
-
- cuMesh::writePLY(argv[2], mesh_data);
+ cuMesh::writePLY(argv[2]+std::string("0.ply"), mesh_data);
+ for(int i = 0; i < 5; i++) {
+ cuMesh::Smoothing::PQGFSmoothing(mesh_data, 2, 1, 1);
+ cuMesh::writePLY(argv[2]+std::to_string(i+1)+".ply", mesh_data);
+ }
 }