From 9b4474612ca7d05cc50e8662623620b30eae719f Mon Sep 17 00:00:00 2001 From: Dtouch Date: Sun, 12 Feb 2023 22:24:10 +0800 Subject: [PATCH] origin version --- .gitignore | 2 + CMakeLists.txt | 16 +++ include/C2C4.h | 50 +++++++++ include/Range.h | 21 ++++ include/SingularityJudger.h | 34 +++++++ include/aabb.h | 35 +++++++ include/gauss_map.h | 80 +++++++++++++++ include/loop_detector.h | 74 ++++++++++++++ include/srf_mesh.h | 19 ++++ main.cpp | 100 ++++++++++++++++++ src/C2C4.cpp | 114 +++++++++++++++++++++ src/Range.cpp | 20 ++++ src/SingularityJudger.cpp | 53 ++++++++++ src/aabb.cpp | 45 ++++++++ src/gauss_map.cpp | 197 ++++++++++++++++++++++++++++++++++++ src/loop_detector.cpp | 145 ++++++++++++++++++++++++++ src/srf_mesh.cpp | 1 + 17 files changed, 1006 insertions(+) create mode 100644 CMakeLists.txt create mode 100644 include/C2C4.h create mode 100644 include/Range.h create mode 100644 include/SingularityJudger.h create mode 100644 include/aabb.h create mode 100644 include/gauss_map.h create mode 100644 include/loop_detector.h create mode 100644 include/srf_mesh.h create mode 100644 main.cpp create mode 100644 src/C2C4.cpp create mode 100644 src/Range.cpp create mode 100644 src/SingularityJudger.cpp create mode 100644 src/aabb.cpp create mode 100644 src/gauss_map.cpp create mode 100644 src/loop_detector.cpp create mode 100644 src/srf_mesh.cpp diff --git a/.gitignore b/.gitignore index e257658..2f2ce88 100644 --- a/.gitignore +++ b/.gitignore @@ -32,3 +32,5 @@ *.out *.app +.idea/ +cmake-build-debug/ \ No newline at end of file diff --git a/CMakeLists.txt b/CMakeLists.txt new file mode 100644 index 0000000..8c3a7fe --- /dev/null +++ b/CMakeLists.txt @@ -0,0 +1,16 @@ +cmake_minimum_required(VERSION 3.21) +project(SingularityJudger) + +set(CMAKE_CXX_STANDARD 14) + +include_directories(include) +# 引入glm tinynurbs依赖 +include_directories(E:/CLib/tinynurbs/include E:/CLib/glm) + +add_executable(SingularityJudger main.cpp + include/loop_detector.h src/loop_detector.cpp + include/gauss_map.h src/gauss_map.cpp + include/aabb.h src/aabb.cpp + include/C2C4.h src/C2C4.cpp + include/Range.h src/Range.cpp + src/SingularityJudger.cpp include/SingularityJudger.h src/srf_mesh.cpp include/srf_mesh.h) diff --git a/include/C2C4.h b/include/C2C4.h new file mode 100644 index 0000000..310a839 --- /dev/null +++ b/include/C2C4.h @@ -0,0 +1,50 @@ +// +// Created by 14727 on 2022/12/7. +// + +#ifndef C2C4_C2C4_H +#define C2C4_C2C4_H + +#include "glm/glm.hpp" +#include "../include/Range.h" +#include "vector" +#include "srf_mesh.h" +#include "tinynurbs/tinynurbs.h" + +using namespace tinynurbs; +using namespace std; + + +class C2C4 { +private: + // 雅可比矩阵【行】【列】<左界，右界> + vector> jacobian; + /** + * 求range去掉某一列的三阶子行列式是否含零 + * @return 若去掉某些列，三阶子行列式含零，输出去掉的列号 + */ + Range determinant3(int c0, int c1, int c2); + Range determinant2(int l0, int l1, int c0, int c1); + + +public: + + C2C4(const SrfMesh& mesh1_, const SrfMesh& mesh2_, RationalSurface srf1_, RationalSurface srf2_); + + const SrfMesh& mesh1; + const SrfMesh& mesh2; + + // 在给定的mesh的最小网格上重新细分的采样数目 + int reSampleCnt_u = 17; + int reSampleCnt_v = 17; + + RationalSurface srf1; + RationalSurface srf2; + + pair c2OrC4(int patchIdx_u1, int patchIdx_v1, int patchIdx_u2, int patchIdx_v2); + vector c4ExcludeCols; + +}; + + +#endif //C2C4_C2C4_H diff --git a/include/Range.h b/include/Range.h new file mode 100644 index 0000000..0fc0046 --- /dev/null +++ b/include/Range.h @@ -0,0 +1,21 @@ +// +// Created by 14727 on 2022/12/7. +// + +#ifndef C2C4_RANGE_H +#define C2C4_RANGE_H + + +class Range { +public: + float a; + float b; + Range(float _a, float _b); + Range()= default; + Range operator + (const Range & r2) const; + Range operator - (const Range & r2) const; + Range operator * (const Range & r2) const; + bool hasZero() const; +}; + +#endif //C2C4_RANGE_H diff --git a/include/SingularityJudger.h b/include/SingularityJudger.h new file mode 100644 index 0000000..aba42f6 --- /dev/null +++ b/include/SingularityJudger.h @@ -0,0 +1,34 @@ +#ifndef SINGULARITYJUDGER_SINGULARITYJUDGER_H +#define SINGULARITYJUDGER_SINGULARITYJUDGER_H + +#include "tinynurbs/tinynurbs.h" +#include "vector" +#include "srf_mesh.h" +#include "tinynurbs/tinynurbs.h" + +using namespace std; +using namespace tinynurbs; + +class SingularityJudger { + +public: + SingularityJudger(RationalSurface &srf1_, RationalSurface &srf2_, SrfMesh &mesh1_, SrfMesh &mesh2_); + + SrfMesh &mesh1; + SrfMesh &mesh2; + // 调用各个工具直接传入mesh信息，避免nurbs曲面的重复采样 + // 但曲面参数信息仍然需要传入，因为C2C4在最小的网格的基础上还需要进一步细分采样 + RationalSurface &srf1; + RationalSurface &srf2; + + // 关注mesh的采样网格中的哪个范围 + // 假设mesh.sampleLevel = 5; + bool judge(pair focusRange_u1, pair focusRange_v1, + pair focusRange_u2, pair focusRange_v2); + + vector> judgeRes; + +}; + + +#endif //SINGULARITYJUDGER_SINGULARITYJUDGER_H diff --git a/include/aabb.h b/include/aabb.h new file mode 100644 index 0000000..58c66f8 --- /dev/null +++ b/include/aabb.h @@ -0,0 +1,35 @@ +// AABB部分的代码直接原封不动用的gitea/GeometryMain/surface-surface-intersection.git中bvh的代码，之后可以合并过去 +#ifndef AABB_HPP +#define AABB_HPP + +#include "glm/glm.hpp" +#include + +class AABB { +public: + // 边界 + glm::dvec3 pMin, pMax; + +public: + AABB() { + double minNum = std::numeric_limits::lowest(); + double maxNum = std::numeric_limits::max(); + pMax = glm::dvec3(minNum, minNum, minNum); + pMin = glm::dvec3(maxNum, maxNum, maxNum); + } + AABB(const glm::dvec3 p) : pMin(p), pMax(p) {} + AABB(const glm::dvec3 p1, const glm::dvec3 p2) { + pMin = glm::dvec3(fmin(p1.x, p2.x), fmin(p1.y, p2.y), fmin(p1.z, p2.z)); + pMax = glm::dvec3(fmax(p1.x, p2.x), fmax(p1.y, p2.y), fmax(p1.z, p2.z)); + } +}; +// 判断点是否在包围盒内，是返回true，否返回false +bool isPointInBox(const glm::dvec3 &pt, const AABB& box); +// aabb包围盒合并操作，包围盒和点 +AABB Union(const AABB& b1, const glm::dvec3& p); +// aabb包围盒合并操作，两个包围盒 +AABB Union(const AABB& b1, const AABB& b2); +// 判断两个aabb包围盒是否重叠 +bool IsOverlap(AABB b1, AABB b2); + +#endif diff --git a/include/gauss_map.h b/include/gauss_map.h new file mode 100644 index 0000000..ab634d4 --- /dev/null +++ b/include/gauss_map.h @@ -0,0 +1,80 @@ +// +// Created by 14727 on 2022/12/24. +// + +#ifndef GAUSSMAP_GAUSS_MAP_H +#define GAUSSMAP_GAUSS_MAP_H + +#include "vector" +#include "tinynurbs/tinynurbs.h" +#include "glm/glm.hpp" +#include "aabb.h" +#include "map" + +// 一个节点就表示一个球上矩形面片 +struct GaussMapNode { + // 四个顶点的法向量值确定一个平面的法向量值 + AABB nBound; + int level, firstChild; +}; + +class GaussMap { +public: + int maxLevel; + int leafSampleCnt; + const std::vector> & normals; + tinynurbs::RationalSurface srf; + std::vector tree; + + void recursiveBuild(int level, int idx, int idx_u, int idx_v); + + GaussMap(const std::vector>& normals_); + + void build(); + + void normalInit(); + + void printQuadTree(); +}; + +std::vector> getOverlapLeafNodes(const GaussMap &gm1, const GaussMap &gm2); + +void recursiveGetOverlapLeafNodes(const GaussMap &gm1, const GaussMap &gm2, int idx1, int idx2, + std::vector> &pairs); + +/** + * 判断两个曲面的Gauss Map在指定的、各自的参数范围内有没有交 + * @param idxRange_u1 第一个gauss map的参数u范围对应的u方向上的采样网格的格子下标范围 + * @param idxRange_v1 第一个gauss map的参数v范围对应的u方向上的采样网格的格子下标范围 + * @param idxRange_u2 第而个gauss map的参数u范围对应的u方向上的采样网格的格子下标范围 + * @param idxRange_v2 第二个gauss map的参数v范围对应的u方向上的采样网格的格子下标范围 + * @return true:gauss map的包围盒有交集，说明gauss map<可能>有重合；false: gauss map一定没有重合 + */ +bool isGaussMapsOverlapped(const GaussMap &gm1, const GaussMap &gm2, std::pair idxRange_u1, + std::pair idxRange_v1, std::pair idxRange_u2, + std::pair idxRange_v2); + +/** + * 判断两个曲面的Gauss Map在指定的、各自的参数范围内有没有交 + * @param range_u1 第一个gauss map的参数u范围 + * @param range_v1 第一个gauss map的参数v范围 + * @param range_u2 第二个gauss map的参数u范围 + * @param range_v2 第二个gauss map的参数v范围 + * @param paramRange_u1 第一个gauss map的参数u定义域 + * @param paramRange_v1 第一个gauss map的参数v定义域 + * @param paramRange_u2 第二个gauss map的参数u定义域 + * @param paramRange_v2 第二个gauss map的参数v定义域 + * @return true:gauss map的包围盒有交集，说明gauss map<可能>有重合；false: gauss map一定没有重合 + */ +bool isGaussMapsOverlapped(const GaussMap &gm1, const GaussMap &gm2, std::pair range_u1, + std::pair range_v1, std::pair range_u2, + std::pair range_v2, std::pair paramRange_u1, + std::pair paramRange_v1, + std::pair paramRange_u2, + std::pair paramRange_v2); + +int getStartIdxOfLayerN(int layer); + +int getChildNodeIdx(int ix, int iy); + +#endif //GAUSSMAP_GAUSS_MAP_H diff --git a/include/loop_detector.h b/include/loop_detector.h new file mode 100644 index 0000000..749c835 --- /dev/null +++ b/include/loop_detector.h @@ -0,0 +1,74 @@ +#ifndef LOOPDETECTOR_LOOP_DETECTOR_H +#define LOOPDETECTOR_LOOP_DETECTOR_H + +#include "tinynurbs/tinynurbs.h" +#include "glm/glm.hpp" + +using namespace std; + +class LoopDetector { +public: + LoopDetector(const vector> &s_evaluation_, const vector> &f_evaluation_, + const vector> &s_tangent_u_, const vector> &f_tangent_u_, + const vector> &s_tangent_v_, const vector> &f_tangent_v_, + const vector> &s_normal_, const vector> &f_normal_); + + // 两个曲面 + tinynurbs::RationalSurface s; + tinynurbs::RationalSurface f; + // 细分采样的总层数（最高与BVH的总层数相等）。从第二层开始，每一层都表示把原来的一个小sub patch分成了更小的四份 + int maxSplitLayer; + + // 采样点的求值和求切向量先设为public，因为这些因该是算好的，LoopDetector不用重复计算 + // 曲面s和f的采样点。 + const vector> &s_evaluation; + const vector> &f_evaluation; + // 曲面s和f的切向量。 + const vector> &s_tangent_u; + const vector> &s_tangent_v; + const vector> &f_tangent_u; + const vector> &f_tangent_v; + // 曲面s和f的法向量 + const vector> &s_normal; + const vector> &f_normal; + + // 有向距离计算结果。有向距离通过采样的方式计算，其结果与s曲面sub patch中的采样网格规模相同，即与s_evaluation大小一致 +// vector> distance; + // 确定有向距离后，对应的f曲面上的最短距离点的位置 + vector>> selectedPointsIdx; + // vector fields, 即有向距离关于u、v的导数 + vector> vectorFields; + vector> rotationNumbers; + +// 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); + + + // 需要做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; + +private: + // 整个曲面一条边上的采样点个数 + int edgeSampleCnt; + + + void initOrientedDistance(); + + void initVectorField(); + + void getRotationNumber(); +}; + + +#endif //LOOPDETECTOR_LOOP_DETECTOR_H diff --git a/include/srf_mesh.h b/include/srf_mesh.h new file mode 100644 index 0000000..f1f0b5f --- /dev/null +++ b/include/srf_mesh.h @@ -0,0 +1,19 @@ +#ifndef SINGULARITYJUDGER_SRF_MESH_H +#define SINGULARITYJUDGER_SRF_MESH_H + +#include "vector" +#include "glm/glm.hpp" +using namespace std; +// 用采样网格表示的一个曲面 +class SrfMesh { +public: + int sampleLevel; // 采样层级；每多一层，uv方向的网格数都分别x2，整个mesh的网格数x4；sampleLevel为1时，采样整个面，仅有一个网格 + int sampleCntOnSingleDir; + vector> evaluation; + vector> tangent_u; + vector> tangent_v; + vector> normal; +}; + + +#endif //SINGULARITYJUDGER_SRF_MESH_H diff --git a/main.cpp b/main.cpp new file mode 100644 index 0000000..2a72e40 --- /dev/null +++ b/main.cpp @@ -0,0 +1,100 @@ +#include +#include "SingularityJudger.h" + +SrfMesh getMesh(const RationalSurface& s, int sampleLevel){ + SrfMesh res; + res.sampleLevel = sampleLevel; + res.sampleCntOnSingleDir = int(pow(2, sampleLevel - 1) + 1); + res.evaluation = vector>(res.sampleCntOnSingleDir, vector(res.sampleCntOnSingleDir)); + res.tangent_u = vector>(res.sampleCntOnSingleDir, vector(res.sampleCntOnSingleDir)); + res.tangent_v = vector>(res.sampleCntOnSingleDir, vector(res.sampleCntOnSingleDir)); + res.normal = vector>(res.sampleCntOnSingleDir, vector(res.sampleCntOnSingleDir)); + + 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(res.sampleCntOnSingleDir - 1); + auto s_step_v = (*(s.knots_v.end() - 1) - s_first_v) / float(res.sampleCntOnSingleDir - 1); + + for(int i = 0; i < res.sampleCntOnSingleDir; i++) { + auto u = s_first_u + s_step_u * float(i); + for(int j = 0; j < res.sampleCntOnSingleDir; j++) { + auto v = s_first_v + s_step_v * float(j); + res.evaluation[i][j] = tinynurbs::surfacePoint(s, u, v); + auto der = tinynurbs::surfaceDerivatives(s, 1, u, v); + res.tangent_u[i][j] = der(1, 0); + res.tangent_v[i][j] = der(0, 1); + res.normal[i][j] = glm::cross(res.tangent_u[i][j], res.tangent_v[i][j]); + } + } +} + + +int main() { + RationalSurface s; + 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, + } + }; + + vector> s_evaluation; + vector> f_evaluation; + // 曲面s和f的切向量。 + vector> s_tangent_u; + vector> s_tangent_v; + vector> f_tangent_u; + const vector> f_tangent_v; + // 曲面s和f的法向量 + const vector> s_normal; + const vector> f_normal; + + auto mesh1 = getMesh(s, 7); + auto mesh2 = getMesh(f, 7); + + SingularityJudger singularityJudger(s, f, mesh1, mesh2); + auto hasSingularity = singularityJudger.judge({2,7}, {4,12}, {3, 9}, {10, 16}); + for(auto line: singularityJudger.judgeRes) { + for(auto el: line) { + cout< + +pair C2C4::c2OrC4(int patchIdx_u1, int patchIdx_v1, int patchIdx_u2, int patchIdx_v2) { + jacobian = vector>(3, vector(4, Range(999999.f, -999999.f))); + float knotsRange_u1 = srf1.knots_u[srf1.knots_u.size() - 1] - srf1.knots_u[0]; + float knotsRange_u2 = srf2.knots_u[srf2.knots_u.size() - 1] - srf2.knots_u[0]; + float knotsRange_v1 = srf1.knots_v[srf1.knots_v.size() - 1] - srf1.knots_v[0]; + float knotsRange_v2 = srf2.knots_v[srf2.knots_v.size() - 1] - srf2.knots_v[0]; + auto patchEdge_u1 = knotsRange_u1 / (mesh1.sampleCntOnSingleDir - 1); + auto patchEdge_u2 = knotsRange_u2 / (mesh2.sampleCntOnSingleDir - 1); + auto patchEdge_v1 = knotsRange_v1 / (mesh1.sampleCntOnSingleDir - 1); + auto patchEdge_v2 = knotsRange_v2 / (mesh2.sampleCntOnSingleDir - 1); + float u1Base = srf1.knots_u[0] + patchIdx_u1 * patchEdge_u1; + float u2Base = srf2.knots_u[0] + patchIdx_u2 * patchEdge_u2; + float v1Base = srf1.knots_v[0] + patchIdx_v1 * patchEdge_v1; + float v2Base = srf2.knots_v[0] + patchIdx_v2 * patchEdge_v2; + // 采样 + for (int i = 0; i < reSampleCnt_u; ++i) { + auto u1 = u1Base + i * patchEdge_u1 / (reSampleCnt_u - 1); + auto u2 = u2Base + i * patchEdge_u2 / (reSampleCnt_u - 1); + for (int j = 0; j < reSampleCnt_v; ++j) { + glm::vec<3, float> der1_u{}, der1_v{}; + glm::vec<3, float> der2_u{}, der2_v{}; + // 尽量减少重复采样：四个边界点的值是已经采样过的 + if (i == 0 && j == 0) { + der1_u = mesh1.tangent_u[patchIdx_u1][patchIdx_v1]; + der1_v = mesh1.tangent_v[patchIdx_u1][patchIdx_v1]; + der2_u = mesh2.tangent_u[patchIdx_u2][patchIdx_v2]; + der2_v = mesh2.tangent_v[patchIdx_u2][patchIdx_v2]; + } else if (i == 0 && j == reSampleCnt_v - 1) { + der1_u = mesh1.tangent_u[patchIdx_u1][patchIdx_v1 + 1]; + der1_v = mesh1.tangent_v[patchIdx_u1][patchIdx_v1 + 1]; + der2_u = mesh2.tangent_u[patchIdx_u2][patchIdx_v2 + 1]; + der2_v = mesh2.tangent_v[patchIdx_u2][patchIdx_v2 + 1]; + } else if (i == reSampleCnt_u - 1 && j == 0) { + der1_u = mesh1.tangent_u[patchIdx_u1 + 1][patchIdx_v1]; + der1_v = mesh1.tangent_v[patchIdx_u1 + 1][patchIdx_v1]; + der2_u = mesh2.tangent_u[patchIdx_u2 + 1][patchIdx_v2]; + der2_v = mesh2.tangent_v[patchIdx_u2 + 1][patchIdx_v2]; + } else if (i == reSampleCnt_u - 1 || j == reSampleCnt_v - 1) { + der1_u = mesh1.tangent_u[patchIdx_u1 + 1][patchIdx_v1 + 1]; + der1_v = mesh1.tangent_v[patchIdx_u1 + 1][patchIdx_v1 + 1]; + der2_u = mesh2.tangent_u[patchIdx_u2 + 1][patchIdx_v2 + 1]; + der2_v = mesh2.tangent_v[patchIdx_u2 + 1][patchIdx_v2 + 1]; + } else { + float v1 = v1Base + j * patchEdge_v1 / (reSampleCnt_v - 1); + float v2 = v2Base + j * patchEdge_v2 / (reSampleCnt_v - 1); + auto der1 = tinynurbs::surfaceDerivatives(srf1, 1, u1, v1); + auto der2 = tinynurbs::surfaceDerivatives(srf2, 1, u2, v2); + der1_u = der1[2], der1_v = der1[1]; + der2_u = -der2[2], der2_v = -der2[1]; + } + jacobian[0][0].a = min(jacobian[0][0].a, der1_u.x), jacobian[0][0].b = max(jacobian[0][0].b, der1_u.x); + jacobian[1][0].a = min(jacobian[1][0].a, der1_u.y), jacobian[1][0].b = max(jacobian[1][0].b, der1_u.y); + jacobian[2][0].a = min(jacobian[2][0].a, der1_u.z), jacobian[2][0].b = max(jacobian[2][0].b, der1_u.z); + + jacobian[0][1].a = min(jacobian[0][1].a, der1_v.x), jacobian[0][1].b = max(jacobian[0][1].b, der1_v.x); + jacobian[1][1].a = min(jacobian[1][1].a, der1_v.y), jacobian[1][1].b = max(jacobian[1][1].b, der1_v.y); + jacobian[2][1].a = min(jacobian[2][1].a, der1_v.z), jacobian[2][1].b = max(jacobian[2][1].b, der1_v.z); + + jacobian[0][2].a = min(jacobian[0][2].a, der2_u.x), jacobian[0][2].b = max(jacobian[0][2].b, der2_u.x); + jacobian[1][2].a = min(jacobian[1][2].a, der2_u.y), jacobian[1][2].b = max(jacobian[1][2].b, der2_u.y); + jacobian[2][2].a = min(jacobian[2][2].a, der2_u.z), jacobian[2][2].b = max(jacobian[2][2].b, der2_u.z); + + jacobian[0][3].a = min(jacobian[0][3].a, der2_v.x), jacobian[0][3].b = max(jacobian[0][3].b, der2_v.x); + jacobian[1][3].a = min(jacobian[1][3].a, der2_v.y), jacobian[1][3].b = max(jacobian[1][3].b, der2_v.y); + jacobian[2][3].a = min(jacobian[2][3].a, der2_v.z), jacobian[2][3].b = max(jacobian[2][3].b, der2_v.z); + } + } + +// 判断存在一个2阶子式行列式非0 + bool hasNotZero = false; // true: 存在一个二阶子式非0 + for (int l0 = 0; l0 < 2 && !hasNotZero; l0++) { + for (int l1 = 1; l1 < 3 && !hasNotZero; l1++) { + for (int c0 = 0; c0 < 3 && !hasNotZero; c0++) { + for (int c1 = 1; c1 < 4 && !hasNotZero; c1++) { + if (!determinant2(l0, l1, c0, c1).hasZero()) hasNotZero = true; + } + } + } + } + +// 判断3阶子式行列式含0情况 + if (!determinant3(0, 1, 2).hasZero()) c4ExcludeCols.push_back(3); + if (!determinant3(0, 1, 3).hasZero()) c4ExcludeCols.push_back(2); + if (!determinant3(0, 2, 3).hasZero()) c4ExcludeCols.push_back(1); + if (!determinant3(1, 2, 3).hasZero()) c4ExcludeCols.push_back(0); + pair res(false, false); + if (hasNotZero && c4ExcludeCols.empty()) res.first = true; // C2 + else if (!c4ExcludeCols.empty()) res.second = true; // C4 + return res; +} + +Range C2C4::determinant3(int c0, int c1, int c2) { + return jacobian[0][c0] * jacobian[1][c1] * jacobian[2][c2] + + jacobian[0][c1] * jacobian[1][c2] * jacobian[2][c0] + + jacobian[0][c2] * jacobian[1][c0] * jacobian[2][c1] - + jacobian[0][c2] * jacobian[1][c1] * jacobian[2][c0] - + jacobian[0][c0] * jacobian[1][c2] * jacobian[2][c1] - + jacobian[0][c1] * jacobian[1][c0] * jacobian[2][c2]; +} + +Range C2C4::determinant2(int l0, int l1, int c0, int c1) { + return jacobian[l0][c0] * jacobian[l1][c1] - jacobian[l0][c1] * jacobian[l1][c0]; +} + +C2C4::C2C4(const SrfMesh &mesh1_, const SrfMesh &mesh2_, RationalSurface srf1_, RationalSurface srf2_) : + mesh1(mesh1_), mesh2(mesh2_), srf1(std::move(srf1_)), srf2(std::move(srf2_)) {} diff --git a/src/Range.cpp b/src/Range.cpp new file mode 100644 index 0000000..1ca44d9 --- /dev/null +++ b/src/Range.cpp @@ -0,0 +1,20 @@ +#include "../include/Range.h" +#include "cmath" + +Range::Range(float _a, float _b): a(_a), b(_b) {} + +Range Range::operator-(const Range &r2) const { + return {a + r2.a, b + r2.b}; +} + +Range Range::operator+(const Range &r2) const { + return {a-r2.b, b-r2.a}; +} + +Range Range::operator*(const Range &r2) const { + return {fmin(fmin(fmin(a * r2.a, a * r2.b), b * r2.a), b * r2.b), fmin(fmin(fmin(a * r2.a, a * r2.b), b * r2.a), b * r2.b)}; +} + +bool Range::hasZero() const { + return a < 0 && b > 0; +} \ No newline at end of file diff --git a/src/SingularityJudger.cpp b/src/SingularityJudger.cpp new file mode 100644 index 0000000..0fcaf70 --- /dev/null +++ b/src/SingularityJudger.cpp @@ -0,0 +1,53 @@ +#include "SingularityJudger.h" + +#include +#include "gauss_map.h" +#include "loop_detector.h" +#include "C2C4.h" + + +bool SingularityJudger::judge(pair focusRange_u1, pair focusRange_v1, pair focusRange_u2, + pair focusRange_v2) { + // TODO gauss map to exclude patch pair that must not make loop or singular intersection + // gauss map 只需要法向量信息，不需要evaluations + GaussMap gaussMap1(mesh1.normal); + GaussMap gaussMap2(mesh2.normal); + if (!isGaussMapsOverlapped(gaussMap1, gaussMap2, focusRange_u1, focusRange_u2, focusRange_v1, focusRange_v2)) { + // 一定没有交 + return false; + } + // TODO loop detection to retain patch pair that must have loop or singular intersection + LoopDetector loopDetector(mesh1.evaluation, mesh2.evaluation, mesh1.tangent_u, mesh2.tangent_u, mesh1.tangent_v, + mesh2.tangent_v, mesh1.normal, mesh2.normal); + loopDetector.detect(focusRange_u1, focusRange_v1, focusRange_u2, focusRange_v2); + judgeRes = vector>(loopDetector.s_subPatchEdgeSampleCnt_u - 1, + vector(loopDetector.s_subPatchEdgeSampleCnt_v - 1, 0)); + // TODO c2 determination for tangency case + bool hasLoop = false; + C2C4 c2C4(mesh1, mesh2, srf1, srf2); + for (int i = 0; i < judgeRes.size(); i++) { + for (int j = 0; j < judgeRes[0].size(); j++) { + if (loopDetector.rotationNumbers[i][j] != 0) { + + hasLoop = true; + + // 非零，有loop + // 依次测试C2 + auto accordPointOnF = loopDetector.selectedPointsIdx[i][j]; // 有向距离最短的，f曲面上的对应面片的坐标 + if (c2C4.c2OrC4(focusRange_u1.first + i, focusRange_v1.first + j, + focusRange_u2.first + accordPointOnF.first, + focusRange_v2.first + accordPointOnF.second).first) { + //C2 + judgeRes[i][j] = -1; + } else judgeRes[i][j] = 1; // 圆环 + } else { + judgeRes[i][j] = 0; // 0表示啥也没有 + } + } + } + return hasLoop; +} + +SingularityJudger:: +SingularityJudger(RationalSurface &srf1_, RationalSurface &srf2_, SrfMesh &mesh1_, SrfMesh &mesh2_) +:srf1(srf1_), srf2(srf2_), mesh1(mesh1_), mesh2(mesh2_){} diff --git a/src/aabb.cpp b/src/aabb.cpp new file mode 100644 index 0000000..12d7cdf --- /dev/null +++ b/src/aabb.cpp @@ -0,0 +1,45 @@ +// AABB部分的代码直接原封不动用的gitea/GeometryMain/surface-surface-intersection.git中bvh的代码，之后可以合并过去 +#include "aabb.h" + +bool isPointInBox(const glm::dvec3 &pt, const AABB& box) { + double x = pt.x; + double y = pt.y; + double z = pt.z; + if (x < box.pMin.x || x > box.pMax.x) return false; + if (y < box.pMin.y || y > box.pMax.y) return false; + if (z < box.pMin.z || z > box.pMax.z) return false; + return true; +} + +AABB Union(const AABB& b1, const glm::dvec3& p) { + AABB ret; + ret.pMin = glm::dvec3(fmin(b1.pMin.x, p.x), + fmin(b1.pMin.y, p.y), + fmin(b1.pMin.z, p.z)); + ret.pMax = glm::dvec3(fmax(b1.pMax.x, p.x), + fmax(b1.pMax.y, p.y), + fmax(b1.pMax.z, p.z)); + return ret; +} + +AABB Union(const AABB& b1, const AABB& b2) { + AABB ret; + ret.pMin = glm::dvec3(fmin(b1.pMin.x, b2.pMin.x), + fmin(b1.pMin.y, b2.pMin.y), + fmin(b1.pMin.z, b2.pMin.z)); + ret.pMax = glm::dvec3(fmax(b1.pMax.x, b2.pMax.x), + fmax(b1.pMax.y, b2.pMax.y), + fmax(b1.pMax.z, b2.pMax.z)); + return ret; +} + +bool IsOverlap(AABB b1, AABB b2) { + if (b1.pMin.x > b2.pMax.x) return false; + if (b2.pMin.x > b1.pMax.x) return false; + if (b1.pMin.y > b2.pMax.y) return false; + if (b2.pMin.y > b1.pMax.y) return false; + if (b1.pMin.z > b2.pMax.z) return false; + if (b2.pMin.z > b1.pMax.z) return false; + return true; +} + diff --git a/src/gauss_map.cpp b/src/gauss_map.cpp new file mode 100644 index 0000000..0cfe199 --- /dev/null +++ b/src/gauss_map.cpp @@ -0,0 +1,197 @@ +// +// Created by 14727 on 2022/12/24. +// + +#include "gauss_map.h" + +#include + + +GaussMap::GaussMap(const std::vector> &normals_): normals(normals_) { + leafSampleCnt = int(normals_.size()); + maxLevel = int(log2(leafSampleCnt - 1) + 1); + tree.resize(int((pow(4, maxLevel) - 1)) / 3); +} + +void GaussMap::recursiveBuild(int level, int idx, int idx_u, int idx_v) { + AABB bound; + if (level == maxLevel) { + // 叶子节点 + for (int i = 0; i < 2; i++) + for (int j = 0; j < 2; j++) + bound = Union(bound, normals[idx_u + i][idx_v + j]); + tree[idx].nBound = bound; + tree[idx].level = level; + tree[idx].firstChild = -1; + return; + } else { + int firstChild = 4 * idx + 1; + int halfRange = int(std::pow(2, maxLevel - level - 1)); // 当前层的曲面片的边长采样宽度的一半 + recursiveBuild(level + 1, firstChild, idx_u, idx_v); + recursiveBuild(level + 1, firstChild + 1, idx_u, idx_v + halfRange); + recursiveBuild(level + 1, firstChild + 2, idx_u + halfRange, idx_v); + recursiveBuild(level + 1, firstChild + 3, idx_u + halfRange, idx_v + halfRange); + for (int i = 0; i < 4; i++) + bound = Union(bound, tree[firstChild + i].nBound); + tree[idx].level = level; + tree[idx].nBound = bound; + tree[idx].firstChild = firstChild; + } +} + +//void GaussMap::normalInit() { +// auto first_u = *(srf.knots_u.begin()); +// auto first_v = *(srf.knots_v.begin()); +// auto step_u = (*(srf.knots_u.end() - 1) - first_u) / (leafSampleCnt - 1); +// auto step_v = (*(srf.knots_v.end() - 1) - first_v) / (leafSampleCnt - 1); +// normals = std::vector>(leafSampleCnt, std::vector(leafSampleCnt, {0, 0, 0})); +// for (int i = 0; i < leafSampleCnt; i++) { +// auto u = first_u + i * step_u; +// for (int j = 0; j < leafSampleCnt; j++) { +// auto v = first_v + j * step_v; +// auto res = tinynurbs::surfaceDerivatives(srf, 1, u, v); +// auto &du = res(1, 0); +// auto &dv = res(0, 1); +// normals[i][j] = glm::normalize(glm::cross(du, dv)); +//// if(i == 1 && j == 1){ +//// if(du == res[2])printf("Woooowl\n"); +//// printf("x(%f, %f, %f), y(%f, %f, %f)\n", du.x, du.y, du.z, dv.x, dv.y, dv.z); +//// printf("normal(%f, %f, %f)\n", normals[i][j].x, normals[i][j].y, normals[i][j].z); +//// } +// } +// } +//} + +void GaussMap::build() { + recursiveBuild(1, 0, 0, 0); +} + +void GaussMap::printQuadTree() { + for (int l = 1, levelNodesCnt = 1, baseIdx = 0; l <= maxLevel; l++, levelNodesCnt *= 4, baseIdx = baseIdx * 4 + 1) { + for (int biasIdx = 0; biasIdx < levelNodesCnt; biasIdx++) { + int idx = baseIdx + biasIdx; + auto pMax = tree[idx].nBound.pMax; + auto pMin = tree[idx].nBound.pMin; + printf("<<%g, %g, %g>,<%g, %g, %g>> ", pMin.x, pMin.y, pMin.z, pMax.x, pMax.y, pMax.z); + } + printf("\n =============$$$$$$$============= \n"); + } +} + + +std::vector> getOverlapLeafNodes(const GaussMap &gm1, const GaussMap &gm2) { + std::vector> resPairs; + recursiveGetOverlapLeafNodes(gm1, gm2, 0, 0, resPairs); + return resPairs; +} + +void recursiveGetOverlapLeafNodes(const GaussMap &gm1, const GaussMap &gm2, int idx1, int idx2, + std::vector> &pairs) { + GaussMapNode A = gm1.tree[idx1]; + GaussMapNode B = gm2.tree[idx2]; + auto AABBSize = [](AABB aabb) { + return (aabb.pMax.z - aabb.pMin.z) * + (aabb.pMax.y - aabb.pMin.y) * + (aabb.pMax.x - aabb.pMin.x); + }; + + // 两个包围盒不相交，返回 + if (!IsOverlap(A.nBound, B.nBound)) return; + // 相交 + if (A.firstChild == -1 && B.firstChild == -1) { + // 两者都是叶子节点 + pairs.emplace_back(idx1, idx2); + } else if (A.firstChild != -1 && B.firstChild == -1) { + // A是中间结点，B是叶子结点 + for (int i = 0; i < 4; i++) recursiveGetOverlapLeafNodes(gm1, gm2, A.firstChild + i, idx2, pairs); + } else if (A.firstChild == -1 && B.firstChild != -1) { + // A是叶子结点，B是中间结点 + for (int i = 0; i < 4; i++) recursiveGetOverlapLeafNodes(gm1, gm2, idx1, B.firstChild + i, pairs); + } else { + // 都是中间结点 + if (AABBSize(A.nBound) > AABBSize(B.nBound)) { + // A的包围盒更大 + for (int i = 0; i < 4; i++) recursiveGetOverlapLeafNodes(gm1, gm2, A.firstChild + i, idx2, pairs); + } else { + // B的包围盒更大 + for (int i = 0; i < 4; i++) recursiveGetOverlapLeafNodes(gm1, gm2, idx1, B.firstChild + i, pairs); + } + } +} + +bool isGaussMapsOverlapped(const GaussMap &gm1, const GaussMap &gm2, std::pair range_u1, + std::pair range_v1, std::pair range_u2, + std::pair range_v2, std::pair paramRange_u1, + std::pair paramRange_v1, std::pair paramRange_u2, + std::pair paramRange_v2) { + if (gm1.maxLevel != gm2.maxLevel || gm1.maxLevel <= 0) { + printf("BVH Layer error!\n"); + return false; + } + int edgeCellCnt = pow(2, gm1.maxLevel - 1); + // 根据所给参数的范围和参数的定义域范围，获得对应的采样网格中的范围 + auto getIdxRange = [](std::pair range, std::pair paramRange, int edgeCellCnt) { + float paramStep = (paramRange.second - paramRange.first) / edgeCellCnt; + return std::pair({int((range.first - paramRange.first) / paramStep), + int((range.second - paramRange.first) / paramStep)}); + }; + auto idxRange_u1 = getIdxRange(range_u1, paramRange_u1, edgeCellCnt); + auto idxRange_v1 = getIdxRange(range_v1, paramRange_v1, edgeCellCnt); + auto idxRange_u2 = getIdxRange(range_u2, paramRange_u2, edgeCellCnt); + auto idxRange_v2 = getIdxRange(range_v2, paramRange_v2, edgeCellCnt); + return isGaussMapsOverlapped(gm1, gm2, idxRange_u1, idxRange_v1, idxRange_u2, idxRange_v2); +} + +bool isGaussMapsOverlapped(const GaussMap &gm1, const GaussMap &gm2, std::pair idxRange_u1, + std::pair idxRange_v1, std::pair idxRange_u2, + std::pair idxRange_v2) { + if (gm1.maxLevel != gm2.maxLevel || gm1.maxLevel <= 0) { + printf("BVH Layer error!\n"); + return false; + } + int commonMaxLayer = gm1.maxLevel; + int edgeCellCnt = pow(2, commonMaxLayer - 1); + if (idxRange_u1.first < 0 || idxRange_u2.first < 0 || idxRange_v1.first < 0 || idxRange_v2.first < 0 || + idxRange_u1.second >= edgeCellCnt || idxRange_u2.second >= edgeCellCnt || + idxRange_v1.second >= edgeCellCnt || idxRange_v2.second >= edgeCellCnt) { + printf("Error when detecting overlapping: idx range invalid!\n"); + return false; + } + auto getRangedBox = [&commonMaxLayer](const GaussMap &gm, const std::pair &idxRange_u, + const std::pair idxRange_v) { + AABB bounding; + for (int i = idxRange_u.first; i <= idxRange_u.second; ++i) { + for (int j = idxRange_v.first; j <= idxRange_v.second; ++j) { + bounding = Union(bounding, + gm.tree[getStartIdxOfLayerN(commonMaxLayer) + getChildNodeIdx(i, j)].nBound); + } + } + return bounding; + }; + return IsOverlap(getRangedBox(gm1, idxRange_u1, idxRange_v1), getRangedBox(gm2, idxRange_u2, idxRange_v2)); +} + +int getStartIdxOfLayerN(int layer) { + if (layer == 1) return 0; + int num = 1; + // 找规律得出 + for (int i = 2; i < layer; i++) { + num = num * 4 + 1; + } + return num; +} + +int getChildNodeIdx(int ix, int iy) { + int idx = 0; + while (ix > 0) { + int log2XCeil = int(log2f(ix)); + idx += int(pow(4, log2XCeil)); + ix -= int(pow(2, log2XCeil)); + } + while (iy > 0) { + int log2XCeil = int(log2f(iy)); + idx += 2 * int(pow(4, log2XCeil)); + iy -= int(pow(2, log2XCeil)); + } + return idx; +} \ No newline at end of file diff --git a/src/loop_detector.cpp b/src/loop_detector.cpp new file mode 100644 index 0000000..d9a95aa --- /dev/null +++ b/src/loop_detector.cpp @@ -0,0 +1,145 @@ +#include "loop_detector.h" +#include + +const float PI = 3.1415927; + +LoopDetector:: +LoopDetector(const vector> &s_evaluation_, const vector> &f_evaluation_, + const vector> &s_tangent_u_, const vector> &f_tangent_u_, + const vector> &s_tangent_v_, const vector> &f_tangent_v_, + const vector> &s_normal_, const vector> &f_normal_) + : s_evaluation(s_evaluation_), f_evaluation(f_evaluation_), + s_tangent_u(s_tangent_u_), s_tangent_v(s_tangent_v_), + f_tangent_u(f_tangent_u_), f_tangent_v(f_tangent_v_), + s_normal(s_normal_), f_normal(f_normal_) { +} + +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[s_subPatchIdxRange_u.first + i][s_subPatchIdxRange_v.first + j], + f_evaluation[f_subPatchIdxRange_u.first + k][f_subPatchIdxRange_v.first + 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[s_subPatchIdxRange_u.first + i][s_subPatchIdxRange_v.first + j] - + f_evaluation[f_subPatchIdxRange_u.first + fPointIdx.first] + [f_subPatchIdxRange_v.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(s_normal[s_subPatchIdxRange_u.first + i][s_subPatchIdxRange_v.first + j] + + f_normal[f_subPatchIdxRange_u.first + fPointIdx.first] + [f_subPatchIdxRange_v.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[s_subPatchIdxRange_u.first + i][s_subPatchIdxRange_v.first + j] + , ps_pv = s_tangent_v[s_subPatchIdxRange_u.first + i][s_subPatchIdxRange_v.first + j]; + auto pf_pp = f_tangent_u[f_subPatchIdxRange_u.first + fPointIdx.first] + [f_subPatchIdxRange_v.first + fPointIdx.second], pf_pq = f_tangent_v[f_subPatchIdxRange_u.first + fPointIdx.first] + [f_subPatchIdxRange_v.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); + + rotationNumbers[i][j] = int(round(rotationNumber)); + + 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; + + // 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; + + initOrientedDistance(); + initVectorField(); + getRotationNumber(); + return {}; +} + diff --git a/src/srf_mesh.cpp b/src/srf_mesh.cpp new file mode 100644 index 0000000..45f5aa8 --- /dev/null +++ b/src/srf_mesh.cpp @@ -0,0 +1 @@ +#include "srf_mesh.h"