origin version

.gitignore

@@ -32,3 +32,5 @@
 *.out
 *.app
 
+.idea/
+cmake-build-debug/

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)

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<vector<Range>> 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<float> srf1_, RationalSurface<float> srf2_);
+
+    const SrfMesh& mesh1;
+    const SrfMesh& mesh2;
+
+    // 在给定的mesh的最小网格上重新细分的采样数目
+    int reSampleCnt_u = 17;
+    int reSampleCnt_v = 17;
+
+    RationalSurface<float> srf1;
+    RationalSurface<float> srf2;
+
+    pair<bool, bool> c2OrC4(int patchIdx_u1, int patchIdx_v1, int patchIdx_u2, int patchIdx_v2);
+    vector<int> c4ExcludeCols;
+
+};
+
+
+#endif //C2C4_C2C4_H

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

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<float> &srf1_, RationalSurface<float> &srf2_, SrfMesh &mesh1_, SrfMesh &mesh2_);
+
+    SrfMesh &mesh1;
+    SrfMesh &mesh2;
+    // 调用各个工具直接传入mesh信息，避免nurbs曲面的重复采样
+    // 但曲面参数信息仍然需要传入，因为C2C4在最小的网格的基础上还需要进一步细分采样
+    RationalSurface<float> &srf1;
+    RationalSurface<float> &srf2;
+
+    // 关注mesh的采样网格中的哪个范围
+    // 假设mesh.sampleLevel = 5;
+    bool judge(pair<int, int> focusRange_u1, pair<int, int> focusRange_v1,
+               pair<int, int> focusRange_u2, pair<int, int> focusRange_v2);
+
+    vector<vector<char>> judgeRes;
+
+};
+
+
+#endif //SINGULARITYJUDGER_SINGULARITYJUDGER_H

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 <cmath>
+
+class AABB {
+public:
+    // 边界
+    glm::dvec3 pMin, pMax;
+
+public:
+    AABB() {
+        double minNum = std::numeric_limits<double>::lowest();
+        double maxNum = std::numeric_limits<double>::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

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<std::vector<glm::vec3>> & normals;
+    tinynurbs::RationalSurface<float> srf;
+    std::vector<GaussMapNode> tree;
+
+    void recursiveBuild(int level, int idx, int idx_u, int idx_v);
+
+    GaussMap(const std::vector<std::vector<glm::vec3>>& normals_);
+
+    void build();
+
+    void normalInit();
+
+    void printQuadTree();
+};
+
+std::vector<std::pair<int, int>> getOverlapLeafNodes(const GaussMap &gm1, const GaussMap &gm2);
+
+void recursiveGetOverlapLeafNodes(const GaussMap &gm1, const GaussMap &gm2, int idx1, int idx2,
+                                  std::vector<std::pair<int, int>> &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<int, int> idxRange_u1,
+                           std::pair<int, int> idxRange_v1, std::pair<int, int> idxRange_u2,
+                           std::pair<int, int> 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<float, float> range_u1,
+                           std::pair<float, float> range_v1, std::pair<float, float> range_u2,
+                           std::pair<float, float> range_v2, std::pair<float, float> paramRange_u1,
+                           std::pair<float, float> paramRange_v1,
+                           std::pair<float, float> paramRange_u2,
+                           std::pair<float, float> paramRange_v2);
+
+int getStartIdxOfLayerN(int layer);
+
+int getChildNodeIdx(int ix, int iy);
+
+#endif //GAUSSMAP_GAUSS_MAP_H

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<vector<glm::vec3>> &s_evaluation_, const vector<vector<glm::vec3>> &f_evaluation_,
+                 const vector<vector<glm::vec3>> &s_tangent_u_, const vector<vector<glm::vec3>> &f_tangent_u_,
+                 const vector<vector<glm::vec3>> &s_tangent_v_, const vector<vector<glm::vec3>> &f_tangent_v_,
+                 const vector<vector<glm::vec3>> &s_normal_, const vector<vector<glm::vec3>> &f_normal_);
+
+    // 两个曲面
+    tinynurbs::RationalSurface<float> s;
+    tinynurbs::RationalSurface<float> f;
+    // 细分采样的总层数（最高与BVH的总层数相等）。从第二层开始，每一层都表示把原来的一个小sub patch分成了更小的四份
+    int maxSplitLayer;
+
+    // 采样点的求值和求切向量先设为public，因为这些因该是算好的，LoopDetector不用重复计算
+    // 曲面s和f的采样点。
+    const vector<vector<glm::vec3>> &s_evaluation;
+    const vector<vector<glm::vec3>> &f_evaluation;
+    // 曲面s和f的切向量。
+    const vector<vector<glm::vec3>> &s_tangent_u;
+    const vector<vector<glm::vec3>> &s_tangent_v;
+    const vector<vector<glm::vec3>> &f_tangent_u;
+    const vector<vector<glm::vec3>> &f_tangent_v;
+    // 曲面s和f的法向量
+    const vector<vector<glm::vec3>> &s_normal;
+    const vector<vector<glm::vec3>> &f_normal;
+
+    // 有向距离计算结果。有向距离通过采样的方式计算，其结果与s曲面sub patch中的采样网格规模相同，即与s_evaluation大小一致
+//    vector<vector<glm::vec3>> distance;
+    // 确定有向距离后，对应的f曲面上的最短距离点的位置
+    vector<vector<pair<int, int>>> selectedPointsIdx;
+    // vector fields, 即有向距离关于u、v的导数
+    vector<vector<glm::vec2>> vectorFields;
+    vector<vector<int>> rotationNumbers;
+
+//    int subPatchEdgeSampleCnt;  // 在一个sub patch的采样网格中，边上的采样点个数
+//    void init(tinynurbs::RationalSurface<float> _s, tinynurbs::RationalSurface<float> _f, int _maxSplitLayer);
+    vector<pair<int, int>> detect(pair<int, int> _s_subPatchIdxRange_u, pair<int, int> _s_subPatchIdxRange_v,
+                                  pair<int, int> _f_subPatchIdxRange_u, pair<int, int> _f_subPatchIdxRange_v);
+
+
+    // 需要做Loop检测的sub patch在最后一层上的下标的范围(每个范围都真包含于[0, 2^(maxSplitLayer-1)-1])
+    pair<int, int> s_subPatchIdxRange_u;
+    pair<int, int> s_subPatchIdxRange_v;
+    pair<int, int> f_subPatchIdxRange_u;
+    pair<int, int> 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

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<vector<glm::vec3>> evaluation;
+    vector<vector<glm::vec3>> tangent_u;
+    vector<vector<glm::vec3>> tangent_v;
+    vector<vector<glm::vec3>> normal;
+};
+
+
+#endif //SINGULARITYJUDGER_SRF_MESH_H

main.cpp

@@ -0,0 +1,100 @@
+#include <iostream>
+#include "SingularityJudger.h"
+
+SrfMesh getMesh(const RationalSurface<float>& s, int sampleLevel){
+    SrfMesh res;
+    res.sampleLevel = sampleLevel;
+    res.sampleCntOnSingleDir = int(pow(2, sampleLevel - 1) + 1);
+    res.evaluation = vector<vector<glm::vec3>>(res.sampleCntOnSingleDir, vector<glm::vec3>(res.sampleCntOnSingleDir));
+    res.tangent_u = vector<vector<glm::vec3>>(res.sampleCntOnSingleDir, vector<glm::vec3>(res.sampleCntOnSingleDir));
+    res.tangent_v = vector<vector<glm::vec3>>(res.sampleCntOnSingleDir, vector<glm::vec3>(res.sampleCntOnSingleDir));
+    res.normal = vector<vector<glm::vec3>>(res.sampleCntOnSingleDir, vector<glm::vec3>(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<float> s;
+    RationalSurface<float> 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<vector<glm::vec3>> s_evaluation;
+    vector<vector<glm::vec3>> f_evaluation;
+    // 曲面s和f的切向量。
+    vector<vector<glm::vec3>> s_tangent_u;
+    vector<vector<glm::vec3>> s_tangent_v;
+    vector<vector<glm::vec3>> f_tangent_u;
+    const vector<vector<glm::vec3>> f_tangent_v;
+    // 曲面s和f的法向量
+    const vector<vector<glm::vec3>> s_normal;
+    const vector<vector<glm::vec3>> 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<<el<<' ';
+        }
+        cout<<endl;
+    }
+    return 0;
+}

src/C2C4.cpp

@@ -0,0 +1,114 @@
+//
+// Created by 14727 on 2022/12/7.
+//
+
+#include "../include/C2C4.h"
+
+#include <utility>
+
+pair<bool, bool> C2C4::c2OrC4(int patchIdx_u1, int patchIdx_v1, int patchIdx_u2, int patchIdx_v2) {
+    jacobian = vector<vector<Range>>(3, vector<Range>(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<bool, bool> 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<float> srf1_, RationalSurface<float> srf2_) :
+        mesh1(mesh1_), mesh2(mesh2_), srf1(std::move(srf1_)), srf2(std::move(srf2_)) {}

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;
+}

src/SingularityJudger.cpp

@@ -0,0 +1,53 @@
+#include "SingularityJudger.h"
+
+#include <utility>
+#include "gauss_map.h"
+#include "loop_detector.h"
+#include "C2C4.h"
+
+
+bool SingularityJudger::judge(pair<int, int> focusRange_u1, pair<int, int> focusRange_v1, pair<int, int> focusRange_u2,
+                              pair<int, int> 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<vector<char>>(loopDetector.s_subPatchEdgeSampleCnt_u - 1,
+                                    vector<char>(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<float> &srf1_, RationalSurface<float> &srf2_, SrfMesh &mesh1_, SrfMesh &mesh2_)
+:srf1(srf1_), srf2(srf2_), mesh1(mesh1_), mesh2(mesh2_){}

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;
+}
+

src/gauss_map.cpp

@@ -0,0 +1,197 @@
+//
+// Created by 14727 on 2022/12/24.
+//
+
+#include "gauss_map.h"
+
+#include <utility>
+
+
+GaussMap::GaussMap(const std::vector<std::vector<glm::vec3>> &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<std::vector<glm::vec3>>(leafSampleCnt, std::vector<glm::vec3>(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<std::pair<int, int>> getOverlapLeafNodes(const GaussMap &gm1, const GaussMap &gm2) {
+    std::vector<std::pair<int, int>> resPairs;
+    recursiveGetOverlapLeafNodes(gm1, gm2, 0, 0, resPairs);
+    return resPairs;
+}
+
+void recursiveGetOverlapLeafNodes(const GaussMap &gm1, const GaussMap &gm2, int idx1, int idx2,
+                                  std::vector<std::pair<int, int>> &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<float, float> range_u1,
+                           std::pair<float, float> range_v1, std::pair<float, float> range_u2,
+                           std::pair<float, float> range_v2, std::pair<float, float> paramRange_u1,
+                           std::pair<float, float> paramRange_v1, std::pair<float, float> paramRange_u2,
+                           std::pair<float, float> 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<float, float> range, std::pair<float, float> paramRange, int edgeCellCnt) {
+        float paramStep = (paramRange.second - paramRange.first) / edgeCellCnt;
+        return std::pair<int, int>({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<int, int> idxRange_u1,
+                           std::pair<int, int> idxRange_v1, std::pair<int, int> idxRange_u2,
+                           std::pair<int, int> 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<int, int> &idxRange_u,
+                                          const std::pair<int, int> 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;
+}

src/loop_detector.cpp

@@ -0,0 +1,145 @@
+#include "loop_detector.h"
+#include <utility>
+
+const float PI = 3.1415927;
+
+LoopDetector::
+LoopDetector(const vector<vector<glm::vec3>> &s_evaluation_, const vector<vector<glm::vec3>> &f_evaluation_,
+             const vector<vector<glm::vec3>> &s_tangent_u_, const vector<vector<glm::vec3>> &f_tangent_u_,
+             const vector<vector<glm::vec3>> &s_tangent_v_, const vector<vector<glm::vec3>> &f_tangent_v_,
+             const vector<vector<glm::vec3>> &s_normal_, const vector<vector<glm::vec3>> &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<vector<pair<int, int>>>(s_subPatchEdgeSampleCnt_u,
+                                                       vector<pair<int, int>>(s_subPatchEdgeSampleCnt_v));
+//    distance = vector<vector<glm::vec3>>(s_subPatchEdgeSampleCnt_u, vector<glm::vec3>(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<vector<glm::vec2>>(s_subPatchEdgeSampleCnt_u, vector<glm::vec2>(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矩阵，见<Detection of loops and singularities of surface intersections> 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<pair<int, int>> LoopDetector::detect(pair<int, int> _s_subPatchIdxRange_u, pair<int, int> _s_subPatchIdxRange_v,
+                                            pair<int, int> _f_subPatchIdxRange_u,
+                                            pair<int, int> _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 {};
+}
+

src/srf_mesh.cpp

@@ -0,0 +1 @@
+#include "srf_mesh.h"