commit a1efe768ea3cf593c5cdffe4477bacd46f675898 Author: Dtouch Date: Sun Jan 8 20:08:29 2023 +0800 origin diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..8617437 --- /dev/null +++ b/.gitignore @@ -0,0 +1,36 @@ +# ---> C++ +# Prerequisites +*.d + +# Compiled Object files +*.slo +*.lo +*.o +*.obj + +# Precompiled Headers +*.gch +*.pch + +# Compiled Dynamic libraries +*.so +*.dylib +*.dll + +# Fortran module files +*.mod +*.smod + +# Compiled Static libraries +*.lai +*.la +*.a +*.lib + +# Executables +*.exe +*.out +*.app + +cmake-build-debug/ +.idea/ \ No newline at end of file diff --git a/CMakeLists.txt b/CMakeLists.txt new file mode 100644 index 0000000..0cb8ae3 --- /dev/null +++ b/CMakeLists.txt @@ -0,0 +1,9 @@ +cmake_minimum_required(VERSION 3.21) +project(LoopDetector) + +set(CMAKE_CXX_STANDARD 14) + +add_executable(LoopDetector main.cpp src/loop_detector.cpp include/loop_detector.h) + +include_directories(E:/CLib/tinynurbs/include E:/CLib/glm) +include_directories(include) \ No newline at end of file diff --git a/README.md b/README.md new file mode 100644 index 0000000..b5d28da --- /dev/null +++ b/README.md @@ -0,0 +1,3 @@ +# GaussMap + +A Gauss mapping tool for Nurbs surfaces on CPUs. \ No newline at end of file diff --git a/include/loop_detector.h b/include/loop_detector.h new file mode 100644 index 0000000..4e92ff8 --- /dev/null +++ b/include/loop_detector.h @@ -0,0 +1,64 @@ +#ifndef LOOPDETECTOR_LOOP_DETECTOR_H +#define LOOPDETECTOR_LOOP_DETECTOR_H + +#include "tinynurbs/tinynurbs.h" +#include "glm/glm.hpp" + +using namespace std; + +class LoopDetector { +public: + // 两个曲面 + tinynurbs::RationalSurface s; + tinynurbs::RationalSurface f; + // 细分采样的总层数（最高与BVH的总层数相等）。从第二层开始，每一层都表示把原来的一个小sub patch分成了更小的四份 + int maxSplitLayer; + + // 采样点的求值和求切向量先设为public，因为这些因该是算好的，LoopDetector不用重复计算 + // 曲面s和f的采样点。 + vector> s_evaluation; + vector> f_evaluation; + // 曲面s和f的切向量。 + vector> s_tangent_u; + vector> s_tangent_v; + vector> f_tangent_u; + vector> f_tangent_v; + +// int subPatchEdgeSampleCnt; // 在一个sub patch的采样网格中，边上的采样点个数 +// void init(tinynurbs::RationalSurface _s, tinynurbs::RationalSurface _f, int _maxSplitLayer); + vector> detect(pair _s_subPatchIdxRange_u, pair _s_subPatchIdxRange_v, + pair _f_subPatchIdxRange_u, pair _f_subPatchIdxRange_v); + +private: + // 需要做Loop检测的sub patch在最后一层上的下标的范围(每个范围都真包含于[0, 2^(maxSplitLayer-1)-1]) + pair s_subPatchIdxRange_u; + pair s_subPatchIdxRange_v; + pair f_subPatchIdxRange_u; + pair f_subPatchIdxRange_v; + // subPatch每条边上的采样点个数 + // 这四个值根据subPatchRange的范围得到，不可被外部修改 + int s_subPatchEdgeSampleCnt_u; + int s_subPatchEdgeSampleCnt_v; + int f_subPatchEdgeSampleCnt_u; + int f_subPatchEdgeSampleCnt_v; + // 整个曲面一条边上的采样点个数 + int edgeSampleCnt; + + // 有向距离计算结果。有向距离通过采样的方式计算，其结果与s曲面sub patch中的采样网格规模相同，即与s_evaluation大小一致 +// vector> distance; + // 确定有向距离后，对应的f曲面上的最短距离点的位置 + vector>> selectedPointsIdx; + // vector fields, 即有向距离关于u、v的导数 + vector> vectorFields; + + void initEvaluation(); + + void initOrientedDistance(); + + void initVectorField(); + + void getRotationNumber(); +}; + + +#endif //LOOPDETECTOR_LOOP_DETECTOR_H diff --git a/main.cpp b/main.cpp new file mode 100644 index 0000000..8c80ed9 --- /dev/null +++ b/main.cpp @@ -0,0 +1,77 @@ +#include +#include "loop_detector.h" + +int main() { + LoopDetector loopDetector; + tinynurbs::RationalSurface s; + tinynurbs::RationalSurface f; + s.degree_u = 3; + s.degree_v = 3; + s.knots_u = {0, 0, 0, 0, 1, 1, 1, 1}; + s.knots_v = {0, 0, 0, 0, 1, 1, 1, 1}; + s.control_points = {4, 4, { + glm::vec3(0, 0.3, 0.9), glm::vec3(0, 0.6, 1), glm::vec3(0, 0.9, 1.1), glm::vec3(0, 1.2, 1), + glm::vec3(0.33, 0.3, 0.12), glm::vec3(0.33, 0.6, 0.12), glm::vec3(0.33, 0.9, 0.12), + glm::vec3(0.33, 1.2, 0.12), + glm::vec3(0.66, 0.3, 0.12), glm::vec3(0.66, 0.6, 0.12), glm::vec3(0.66, 0.9, 0.12), + glm::vec3(0.66, 1.2, 0.12), + glm::vec3(1, 0.3, 0.8), glm::vec3(1, 0.6, 1), glm::vec3(1, 0.9, 1.1), glm::vec3(1, 1.2, 1) + }}; + s.weights = {4, 4, + { + 1, 1, 1, 1, + 1, 1, 1, 1, + 1, 1, 1, 1, + 1, 1, 1, 1, + } + }; + + f.degree_u = 3; + f.degree_v = 3; + f.knots_u = {0, 0, 0, 0, 1, 1, 1, 1}; + f.knots_v = {0, 0, 0, 0, 1, 1, 1, 1}; + f.control_points = {4, 4, { + glm::vec3(0, 0.2, 0.9), glm::vec3(0, 0.5, 1.8), glm::vec3(0, 0.8, 1.1), glm::vec3(0, 1.2, 1), + glm::vec3(0.33, 0.2, 0.12), glm::vec3(0.33, 0.5, 0.42), glm::vec3(0.33, 0.9, -0.62), + glm::vec3(0.33, 1.1, -1.756), + glm::vec3(0.66, 0.2, 0.12), glm::vec3(0.66, 0.5, 0.42), glm::vec3(0.66, 0.9, -0.62), + glm::vec3(0.66, 1.0, -1.756), + glm::vec3(1, 0.2, 0.8), glm::vec3(1, 0.5, 1), glm::vec3(1, 0.9, 1.1), glm::vec3(1, 1.2, 1) + }}; + f.weights = {4, 4, + { + 1, 1, 1, 1, + 1, 1, 1, 1, + 1, 1, 1, 1, + 1, 1, 1, 1, + } + }; + + loopDetector.s = s; + loopDetector.f = f; + loopDetector.maxSplitLayer = 6; + // 需要做Loop检测的sub patch在最后一层上的下标的范围(每个范围都真包含于[0, 2^(maxSplitLayer-1)-1]) + // 这里范围真包含于[0, 31] + loopDetector.detect({3, 11}, {4, 11}, {2, 7}, {6, 15}); + + glm::vec3 a(2, 3, 4); + glm::vec3 b(3, 7, 1); + auto nab = glm::normalize(a - b); + auto res = nab * (a - b); + auto ab = a * b; + cout << res.x << ", " << res.y << ", " << res.z << endl; + cout << nab.x << ", " << nab.y << ", " << nab.z << endl; + cout << ab.x << ", " << ab.y << ", " << ab.z << endl; + std::cout << "Hello, World!" << std::endl; + auto m = glm::mat2x2(1, 2, 3, 4); + cout << m[0][0] << " " << m[0][1] << endl; + auto dotRes = glm::dot(a, b); + cout << dotRes << endl; + + cout< + +const float PI = 3.1415927; + +void LoopDetector::initEvaluation() { + // 初始化subPatch上所有采样点的值，以及它们的切向量，并存储 + // 当该工具集成到系统中后，这一部分应该是可以直接全部去掉的 + // 因为在GPU上已经做过对每个点的最细粒度的evaluation了，只需要从整个曲面把需要的subPatch找出来即可 + + // subPatch每条边上的采样点个数。边上的格子个数=range.second-range.first+1，采样点个数=格子个数+1 + s_subPatchEdgeSampleCnt_u = s_subPatchIdxRange_u.second - s_subPatchIdxRange_u.first + 2; + s_subPatchEdgeSampleCnt_v = s_subPatchIdxRange_v.second - s_subPatchIdxRange_v.first + 2; + f_subPatchEdgeSampleCnt_u = f_subPatchIdxRange_u.second - f_subPatchIdxRange_u.first + 2; + f_subPatchEdgeSampleCnt_v = f_subPatchIdxRange_v.second - f_subPatchIdxRange_v.first + 2; + + edgeSampleCnt = int(pow(2, maxSplitLayer - 1)) + 1; + + s_evaluation = vector>(s_subPatchEdgeSampleCnt_u, vector(s_subPatchEdgeSampleCnt_v)); + f_evaluation = vector>(f_subPatchEdgeSampleCnt_u, vector(f_subPatchEdgeSampleCnt_v)); + // 曲面s和f的切向量。 + s_tangent_u = vector>(s_subPatchEdgeSampleCnt_u, vector(s_subPatchEdgeSampleCnt_v)); + s_tangent_v = vector>(s_subPatchEdgeSampleCnt_u, vector(s_subPatchEdgeSampleCnt_v)); + f_tangent_u = vector>(f_subPatchEdgeSampleCnt_u, vector(f_subPatchEdgeSampleCnt_v)); + f_tangent_v = vector>(f_subPatchEdgeSampleCnt_u, vector(f_subPatchEdgeSampleCnt_v)); + + auto s_first_u = *(s.knots_u.begin()); + auto s_first_v = *(s.knots_v.begin()); + auto s_step_u = (*(s.knots_u.end() - 1) - s_first_u) / float(edgeSampleCnt - 1); + auto s_step_v = (*(s.knots_v.end() - 1) - s_first_v) / float(edgeSampleCnt - 1); + + auto f_first_u = *(f.knots_u.begin()); + auto f_first_v = *(f.knots_v.begin()); + auto f_step_u = (*(f.knots_u.end() - 1) - f_first_u) / float(edgeSampleCnt - 1); + auto f_step_v = (*(f.knots_v.end() - 1) - f_first_v) / float(edgeSampleCnt - 1); + + for (int i = 0; i < s_subPatchEdgeSampleCnt_u; i++) { + auto u = s_first_u + s_step_u * float(s_subPatchIdxRange_u.first + i); + for (int j = 0; j < s_subPatchEdgeSampleCnt_v; j++) { + auto v = s_first_v + s_step_v * float(s_subPatchIdxRange_v.first + j); + s_evaluation[i][j] = tinynurbs::surfacePoint(s, u, v); + auto der = tinynurbs::surfaceDerivatives(s, 1, u, v); + s_tangent_u[i][j] = der(1, 0); + s_tangent_v[i][j] = der(0, 1); + } + } + for (int i = 0; i < f_subPatchEdgeSampleCnt_u; i++) { + auto u = f_first_u + f_step_u * float(f_subPatchIdxRange_u.first + i); + for (int j = 0; j < f_subPatchEdgeSampleCnt_v; j++) { + auto v = f_first_v + f_step_v * float(f_subPatchIdxRange_v.first + j); + f_evaluation[i][j] = tinynurbs::surfacePoint(f, u, v); + auto der = tinynurbs::surfaceDerivatives(f, 1, u, v); + f_tangent_u[i][j] = der(1, 0); + f_tangent_v[i][j] = der(0, 1); + } + } +} + +void LoopDetector::initOrientedDistance() { + + selectedPointsIdx = vector>>(s_subPatchEdgeSampleCnt_u, + vector>(s_subPatchEdgeSampleCnt_v)); +// distance = vector>(s_subPatchEdgeSampleCnt_u, vector(s_subPatchEdgeSampleCnt_v)); + + for (int i = 0; i < s_subPatchEdgeSampleCnt_u; i++) { + for (int j = 0; j < s_subPatchEdgeSampleCnt_v; j++) { + float minDis = FLT_MAX; + int minDisFIdx_u = -1, minDisFIdx_v = -1; + for (int k = 0; k < f_subPatchEdgeSampleCnt_u; k++) { + for (int l = 0; l < f_subPatchEdgeSampleCnt_v; l++) { + auto dis = glm::distance(s_evaluation[i][j], f_evaluation[k][l]); + // 确定f上的对应点时，依据的最小欧氏距离，而不是有向距离 + if (dis < minDis) { + minDis = dis; + minDisFIdx_u = k; + minDisFIdx_v = l; + } + } + } +// auto N = glm::normalize(s_evaluation[i][j] - f_evaluation[minDisFIdx_u][minDisFIdx_v]); +// // 直接normalize得到的单位向量总是与SF同向，但实际上N应当总是与f的一侧法向量成锐角 +// auto testAngle = N * glm::cross(f_tangent_u[minDisFIdx_u][minDisFIdx_v], f_tangent_v[minDisFIdx_u][minDisFIdx_v]); +// N = testAngle.x + testAngle.y + testAngle.z > 0 ? N : -N; +// distance[i][j] = N * (s_evaluation[i][j] - f_evaluation[minDisFIdx_u][minDisFIdx_v]); + selectedPointsIdx[i][j] = {minDisFIdx_u, minDisFIdx_v}; + } + } +} + +void LoopDetector::initVectorField() { + vectorFields = vector>(s_subPatchEdgeSampleCnt_u, vector(s_subPatchEdgeSampleCnt_v)); + for (int i = 0; i < s_subPatchEdgeSampleCnt_u; i++) { + for (int j = 0; j < s_subPatchEdgeSampleCnt_v; j++) { + auto fPointIdx = selectedPointsIdx[i][j]; + auto N = glm::normalize(s_evaluation[i][j] - f_evaluation[fPointIdx.first][fPointIdx.second]); + if (isnan(N.x) || isnan(N.y) || isnan(N.z)) { + // 绝了，曲面s和f的交线正好经过采样点，此时就无法定义向量N + // 这里取两个曲面在交点的法向量的平均做近似 + auto tmp1 = glm::cross(s_tangent_u[i][j], s_tangent_v[i][j]); + auto tmp2 = glm::cross(f_tangent_u[fPointIdx.first][fPointIdx.second], + f_tangent_v[fPointIdx.first][fPointIdx.second]); + N = glm::normalize(glm::cross(s_tangent_u[i][j], s_tangent_v[i][j]) + + glm::cross(f_tangent_u[fPointIdx.first][fPointIdx.second], + f_tangent_v[fPointIdx.first][fPointIdx.second])); + } + // u1、u2两个向量构成和N垂直的曲面 + glm::vec3 u1(1, 1, 1); + if (N.z != 0)u1.z = (-N.x - N.y) / N.z; + else if (N.y != 0)u1.y = (-N.x - N.z) / N.y; + else u1.x = (-N.y - N.z) / N.x; + u1 = glm::normalize(u1); + auto u2 = glm::normalize(glm::cross(N, u1)); + // s，f曲面在两个方向上的偏导 + auto ps_pu = s_tangent_u[i][j], ps_pv = s_tangent_v[i][j]; + auto pf_pp = f_tangent_u[fPointIdx.first][fPointIdx.second], pf_pq = f_tangent_v[fPointIdx.first][fPointIdx.second]; + // 构造Aij矩阵，见 APPENDIX (A7) + auto inverseAij = glm::inverse( + glm::mat2x2(glm::dot(u1, pf_pp), glm::dot(u1, pf_pq), glm::dot(u2, pf_pp), glm::dot(u2, pf_pq)) + ); + auto mBmn = glm::mat2x2(glm::dot(u1, ps_pu), glm::dot(u1, ps_pv), glm::dot(u2, ps_pu), glm::dot(u2, ps_pv)); + auto tmp = glm::vec2(glm::dot(N, pf_pp), glm::dot(N, pf_pq)); + auto vNfpNfq_inverseAij = glm::vec2(glm::dot(N, pf_pp), glm::dot(N, pf_pq)) * inverseAij; + auto pd_pu = glm::dot(N, ps_pu) - + glm::dot(vNfpNfq_inverseAij, glm::vec2(glm::dot(u1, ps_pu), glm::dot(u2, ps_pu))); + auto pd_pv = glm::dot(N, ps_pv) - + glm::dot(vNfpNfq_inverseAij, glm::vec2(glm::dot(u1, ps_pv), glm::dot(u2, ps_pv))); + vectorFields[i][j] = glm::vec2(pd_pu, pd_pv); + } + } + int a = 1; + int b = 2; +} + +void LoopDetector::getRotationNumber() { + // 以小格子为单位遍历 + for (int i = 0; i < s_subPatchEdgeSampleCnt_u - 1; i++) { + for (int j = 0; j < s_subPatchEdgeSampleCnt_v - 1; j++) { + auto rotationNumber = 0.; + + for (int biasI = 0; biasI < 2; biasI++) { + for (int biasJ = 0; biasJ < 2; biasJ++) { + auto v = vectorFields[i + biasI][j + biasJ]; + rotationNumber += atan(double(v.x / v.y)); + } + } + auto v00 = vectorFields[i][j]; + auto v01 = vectorFields[i][j + 1]; + auto v11 = vectorFields[i + 1][j + 1]; + auto v10 = vectorFields[i + 1][j]; + + rotationNumber = 2 * atan(double (v11.x/v11.y)) - 2 * atan(double(v00.x/v00.y)); + rotationNumber /= (2 * PI); + printf("rotationNumber: %lf\n", rotationNumber); + } + } +} + +vector> LoopDetector::detect(pair _s_subPatchIdxRange_u, pair _s_subPatchIdxRange_v, + pair _f_subPatchIdxRange_u, + pair _f_subPatchIdxRange_v) { + s_subPatchIdxRange_u = _s_subPatchIdxRange_u; + s_subPatchIdxRange_v = _s_subPatchIdxRange_v; + f_subPatchIdxRange_u = _f_subPatchIdxRange_u; + f_subPatchIdxRange_v = _f_subPatchIdxRange_v; + initEvaluation(); // 当evaluation已经完成，这一步就不用了，直接传入两个evaluation矩阵和四个tangent矩阵的值即可 + initOrientedDistance(); + initVectorField(); + getRotationNumber(); + return {}; +} \ No newline at end of file