优化&测试

CMakeLists.txt

@@ -17,6 +17,14 @@ add_executable(SingularityJudger main.cpp
         include/aabb.h src/aabb.cpp
         include/C2C4.h src/C2C4.cpp
         include/Range.h src/Range.cpp
+        include/bvh.h src/bvh.cpp
+        include/utils.h src/utils.cpp
         src/SingularityJudger.cpp include/SingularityJudger.h src/srf_mesh.cpp include/srf_mesh.h)
 
+# 将intersectTest的内容拷贝到build文件夹下
+# 参考https://stackoverflow.com/questions/13429656/how-to-copy-contents-of-a-directory-into-build-directory-after-make-with-cmake
+add_custom_command(TARGET SingularityJudger PRE_BUILD
+        COMMAND ${CMAKE_COMMAND} -E copy_directory
+        ${CMAKE_SOURCE_DIR}/intersectTest/ $<TARGET_FILE_DIR:SingularityJudger>/intersectTest)
+
 target_link_libraries(SingularityJudger PRIVATE Qt${QT_VERSION_MAJOR}::Widgets)

include/SingularityJudger.h

@@ -4,6 +4,7 @@
 #include "tinynurbs/tinynurbs.h"
 #include "vector"
 #include "srf_mesh.h"
+#include "gauss_map.h"
 #include "tinynurbs/tinynurbs.h"
 
 using namespace std;
@@ -14,6 +15,9 @@ class SingularityJudger {
 public:
     SingularityJudger(RationalSurface<float> &srf1_, RationalSurface<float> &srf2_, SrfMesh &mesh1_, SrfMesh &mesh2_);
 
+    double gaussMapTimeCost;
+
+
     SrfMesh &mesh1;
     SrfMesh &mesh2;
     // 调用各个工具直接传入mesh信息，避免nurbs曲面的重复采样
@@ -21,16 +25,25 @@ public:
     RationalSurface<float> &srf1;
     RationalSurface<float> &srf2;
 
+    void init();
+
     // 关注mesh的采样网格中的哪个范围
     // 假设mesh.sampleLevel = 5, 则sampleRange 属于[0, 15];
     bool judge(pair<int, int> focusRange_u1, pair<int, int> focusRange_v1,
                pair<int, int> focusRange_u2, pair<int, int> focusRange_v2);
+
     // 没有loop detection的版本
     void judge2(pair<int, int> focusRange_u1, pair<int, int> focusRange_v1,
+                pair<int, int> focusRange_u2, pair<int, int> focusRange_v2);
+
+    void judge(const vector<pair<pair<int, int>, pair<int, int>>> &intersectBoxPairs,
+               pair<int, int> focusRange_u1, pair<int, int> focusRange_v1,
                pair<int, int> focusRange_u2, pair<int, int> focusRange_v2);
 
     vector<vector<char>> judgeRes;
-
+private:
+    GaussMap gaussMap1;
+    GaussMap gaussMap2;
 };
 
 

include/bvh.h

@@ -0,0 +1,27 @@
+#ifndef BVH_HPP
+#define BVH_HPP
+
+#include "aabb.h"
+#include <tinynurbs/tinynurbs.h>
+#include <vector>
+
+struct BVHNode {
+    AABB bound;    // AABB包围盒
+    int firstChild, level;  // 第一个孩子节点的下标，该结点所在层次
+    int idx_u;
+    int idx_v;
+};
+
+class BVH {
+public:
+    int maxLevel;         // 树高和结点数量
+    const std::vector<std::vector<glm::vec3>> & evaluations;
+    std::vector<BVHNode> tree;
+    explicit BVH (const std::vector<std::vector<glm::vec3>>& evaluations_);
+    void recursiveBuild(int level, int idx, int idx_u, int idx_v);
+    void build();
+};
+void recursiveGetOverlapLeafNodes(const BVH &bvh1, const BVH &bvh2, int idx1, int idx2,
+                                  std::vector<std::pair<int, int>> &pairs);
+std::vector<std::pair<std::pair<int, int>, std::pair<int,int>>> getOverlapLeafNodes(const BVH &bvh1, const BVH &bvh2);
+#endif

include/gauss_map.h

@@ -14,7 +14,7 @@
 // 一个节点就表示一个球上矩形面片
 struct GaussMapNode {
     // 四个顶点的法向量值确定一个平面的法向量值
-    AABB nBound;
+    AABB bound;
     int level, firstChild;
 };
 
@@ -22,17 +22,14 @@ class GaussMap {
 public:
     int maxLevel;
     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_);
+    explicit GaussMap(const std::vector<std::vector<glm::vec3>>& normals_);
 
     void build();
 
-    void normalInit();
-
     void printQuadTree();
 };
 
@@ -72,8 +69,14 @@ bool isGaussMapsOverlapped(const GaussMap &gm1, const GaussMap &gm2, std::pair<f
                            std::pair<float, float> paramRange_u2,
                            std::pair<float, float> paramRange_v2);
 
+/**
+ * 当前层的第一个节点的坐标
+ * @param layer 层号。根节点为第一层
+ */
 int getStartIdxOfLayerN(int layer);
-
+/**
+ * 根据u、v参数域中的小矩形的位置，判断它在BVH叶节点层（也就是最后一层）中的位置
+ */
 int getChildNodeIdx(int ix, int iy);
 
 #endif //GAUSSMAP_GAUSS_MAP_H

include/loop_detector.h

@@ -44,6 +44,9 @@ public:
 //    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);
+    vector<pair<int, int>> detect(const vector<pair<pair<int, int>, pair<int, int>>> &intersectBoxPairs,
+                                  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])
@@ -64,6 +67,8 @@ private:
 
     void initOrientedDistance();
 
+    void initOrientedDisAndVecFields(const vector<pair<pair<int, int>, pair<int, int>>> &intersectBoxPairs);
+
     // Yawei Ma的文章中的做法
     void initVectorField0();
 
@@ -71,6 +76,7 @@ private:
     void initVectorField();
 
     void getRotationNumber();
+    void getRotationNumber(const vector<pair<pair<int, int>, pair<int, int>>> &intersectBoxPairs);
 };
 
 

include/utils.h

@@ -0,0 +1,23 @@
+#ifndef UNTITLED1_UTILS_H
+#define UNTITLED1_UTILS_H
+
+
+#define IN_UNIX 0  // 确定当前运行的操作系统（需要通过系统调用获得时间）
+namespace utils {
+#if IN_UNIX
+
+#include <sys/time.h>
+#include <ctime>
+
+double get_time();
+
+#else
+#include <windows.h>
+
+    double get_time();
+
+#endif
+}
+
+
+#endif //UNTITLED1_UTILS_H

intersectTest/case10/surfaces.txt

@@ -0,0 +1,3 @@
+W:=Matrix([[1.01,0.99,1,1.1],[1,1,1,0.98],[1,1,1,0.9],[1,1,1.1,1.1]]);
+> P1:=Matrix(4, 4, {(1, 1) = [1, -1, -1], (1, 2) = [1, -1, -1/3], (1, 3) = [1, -1, 1/3], (1, 4) = [1, -1, 1], (2, 1) = [-1/3, 1, -1], (2, 2) = [-1/3, 1, -1/3], (2, 3) = [-1/3, 1, 1/3], (2, 4) = [-1/3, 1, 1], (3, 1) = [-1/3, -1, -1], (3, 2) = [-1/3, -1, -1/3], (3, 3) = [-1/3, -1, 1/3], (3, 4) = [-1/3, -1, 1], (4, 1) = [1, 1, -1], (4, 2) = [1, 1, -1/3], (4, 3) = [1, 1, 1/3], (4, 4) = [1, 1, 1]});
+> P2:=Matrix(4, 4, {(1, 1) = [-1, 0, -1], (1, 2) = [-1, 0, -1/3], (1, 3) = [-1, 0, 1/3], (1, 4) = [-1, 0, 1], (2, 1) = [-1/3, 0, -1], (2, 2) = [-1/3, 0, -1/3], (2, 3) = [-1/3, 0, 1/3], (2, 4) = [-1/3, 0, 1], (3, 1) = [1/3, 0, -1], (3, 2) = [1/3, 0, -1/3], (3, 3) = [1/3, 0, 1/3], (3, 4) = [1/3, 0, 1], (4, 1) = [1, 0, -1], (4, 2) = [1, 0, -1/3], (4, 3) = [1, 0, 1/3], (4, 4) = [1, 0, 1]});

intersectTest/case12/surfaces.txt

@@ -0,0 +1,3 @@
+>W:=Matrix([[1.01,0.99,1,1.1],[1,1,1,0.98],[1,1,1,0.9],[1,1,1.1,1.1]]);
+> P1:=Matrix(4, 4, {(1, 1) = [0.5625, -0.703125, -1.25], (1, 2) = [0.5625, -0.703125, -0.4167], (1, 3) = [0.5625, -0.703125, 0.4167], (1, 4) = [0.5625, -0.703125, 1.25], (2, 1) = [-1.52083, 2.36979167, -1.25], (2, 2) = [-1.52083, 2.36979167, -0.4167], (2, 3) = [-1.52083, 2.36979167, 0.4167], (2, 4) = [-1.52083, 2.36979167, 1.25], (3, 1) = [-1.52083, -2.36979167, -1.25], (3, 2) = [-1.52083, -2.36979167, -0.4167], (3, 3) = [-1.52083, -2.36979167, 0.4167], (3, 4) = [-1.52083, -2.36979167, 1.25], (4, 1) = [0.5625, 0.703125, -1.25], (4, 2) = [0.5625, 0.703125, -0.4167], (4, 3) = [0.5625, 0.703125, 0.4167], (4, 4) = [0.5625, 0.703125, 1.25]});
+> P2:=Matrix(4, 4, {(1, 1) = [-1.5625, -1.25, -1.953125], (1, 2) = [-1.5625, -0.4167, -1.953125], (1, 3) = [-1.5625, 0.4167, -1.953125], (1, 4) = [-1.5625, 1.25, -1.953125], (2, 1) = [0.52083, -1.25, 1.953125], (2, 2) = [0.52083, -0.4167, 1.953125], (2, 3) = [0.52083, 0.4167, 1.953125], (2, 4) = [0.52083, 1.25, 1.953125], (3, 1) = [0.52083, -1.25, -1.953125], (3, 2) = [0.52083, -0.4167, -1.953125], (3, 3) = [0.52083, 0.4167, -1.953125], (3, 4) = [0.52083, 1.25, -1.953125], (4, 1) = [-1.5625, -1.25, 1.953125], (4, 2) = [-1.5625, -0.4167, 1.953125], (4, 3) = [-1.5625, 0.4167, 1.953125], (4, 4) = [-1.5625, 1.25, 1.953125]});

intersectTest/case14/surfaces.txt

@@ -0,0 +1,3 @@
+W:=Matrix([[1.01,0.99,1,1.1],[1,1,1,0.98],[1,1,1,0.9],[1,1,1.1,1.1]]);
+> P1:=Matrix(4, 4, {(1, 1) = [0., .3000000119, .8999999762], (1, 2) = [0., .6000000239, 1.000000000], (1, 3) = [0., .9000000358, 1.100000024], (1, 4) = [0., 1.200000048, 1.], (2, 1) = [.3300000131, .3000000119, .1199999973], (2, 2) = [.3300000131, .6000000239, .1199999973], (2, 3) = [.3300000131, .9000000358, .1199999973], (2, 4) = [.3300000131, 1.200000048, .1199999973], (3, 1) = [.6600000262, .3000000119, .1199999973], (3, 2) = [.6600000262, .6000000239, .1199999973], (3, 3) = [.6600000262, .9000000358, .1199999973], (3, 4) = [.6600000262, 1.200000048, .1199999973], (4, 1) = [1.000000000, .3000000119, .8000000119], (4, 2) = [1.000000000, .6000000239, 1.000000000], (4, 3) = [1.000000000, .9000000358, 1.100000024], (4, 4) = [1.000000000, 1.200000048, 1.]});
+> P2:=Matrix(4, 4, {(1, 1) = [0., .3000000119, .8999999762], (1, 2) = [0., .6000000239, 1.000000000], (1, 3) = [0., .9000000358, 1.100000024], (1, 4) = [0., 1.200000048, 1.], (2, 1) = [.3300000131, .3000000119, .1199999973], (2, 2) = [.3300000131, .6000000239, .4199999869], (2, 3) = [.3300000131, .9000000358, -.6200000048], (2, 4) = [.3300000131, 1.200000048, -1.756000042], (3, 1) = [.6600000262, .3000000119, .1199999973], (3, 2) = [.6600000262, .6000000239, .4199999869], (3, 3) = [.6600000262, .9000000358, -.6200000048], (3, 4) = [.6600000262, 1.200000048, -1.756000042], (4, 1) = [1.000000000, .3000000119, .8000000119], (4, 2) = [1.000000000, .6000000239, 1.000000000], (4, 3) = [1.000000000, .9000000358, 1.100000024], (4, 4) = [1.000000000, 1.200000048, 1.]});

intersectTest/case17/surfaces.txt

@@ -0,0 +1,3 @@
+W:=Matrix([[1,1,1,1],[1,1,1,1],[1,1,1,1],[1,1,1,1]]);
+P1:=Matrix(4, 4, {(1, 1) = [-3., -6., 0.], (1, 2) = [-1., -6., 4.], (1, 3) = [1., -6., 4.], (1, 4) = [3., -6., 0.], (2, 1) = [-3., -2., 0.], (2, 2) = [-1., -2., 4.], (2, 3) = [1., -2., 4.], (2, 4) = [3., -2., 0.], (3, 1) = [-3., 2., 0.], (3, 2) = [-1., 2., 4.], (3, 3) = [1., 2., 4.], (3, 4) = [3., 2., 0.], (4, 1) = [-3., 6., 0.], (4, 2) = [-1., 6., 4.], (4, 3) = [1., 6., 4.], (4, 4) = [3., 6., 0.]});
+P2:=Matrix(4, 4, {(1, 1) = [-6., -3., 0.], (1, 2) = [-2., -3., 0.], (1, 3) = [2., -3., 0.], (1, 4) = [6., -3., 0.], (2, 1) = [-6., -1., 4.], (2, 2) = [-2., -1., 4.], (2, 3) = [2., -1., 4.], (2, 4) = [6., -1., 4.], (3, 1) = [-6., 1., 4.], (3, 2) = [-2., 1., 4.], (3, 3) = [2., 1., 4.], (3, 4) = [6., 1., 4.], (4, 1) = [-6., 3., 0.], (4, 2) = [-2., 3., 0.], (4, 3) = [2., 3., 0.], (4, 4) = [6., 3., 0.]});

intersectTest/case2/surfaces.txt

@@ -0,0 +1,3 @@
+W:=Matrix([[1.1,0.9,1.1,0.9,1,1],[1,1,1,1.1,1,0.9],[0.9,1,1.1,1,0.9,1],[0.9,0.9,1,0.9,1,1],[1.1,1,1.1,1,1,0.9],[1.1,1,1,1,1.2,1.1]]);
+P1 := Matrix(6, 6, {(1, 1) = [0, 0, .2], (1, 2) = [0, .2, .3], (1, 3) = [0, .4, .25], (1, 4) = [0, .6, .12], (1, 5) = [0, .8, .22], (1, 6) = [0, 1, .12], (2, 1) = [.2, 0, .22], (2, 2) = [.2, .2, -.3], (2, 3) = [.2, .4, .4], (2, 4) = [.2, .6, -.1], (2, 5) = [.2, .8, -.1], (2, 6) = [.2, 1, .2], (3, 1) = [.4, 0, .18], (3, 2) = [.4, .2, -.3], (3, 3) = [.4, .4, .4], (3, 4) = [.4, .6, -.1], (3, 5) = [.4, .8, -.1], (3, 6) = [.4, 1, .23], (4, 1) = [.6, 0, .2], (4, 2) = [.6, .2, .4], (4, 3) = [.6, .4, .4], (4, 4) = [.6, .6, .4], (4, 5) = [.6, .8, .4], (4, 6) = [.6, 1, .2], (5, 1) = [.8, 0, .19], (5, 2) = [.8, .2, -.4], (5, 3) = [.8, .4, .4], (5, 4) = [.8, .6, -.27], (5, 5) = [.8, .8, -.27], (5, 6) = [.8, 1, .2], (6, 1) = [1, 0, .2], (6, 2) = [1, .2, .12], (6, 3) = [1, .4, .22], (6, 4) = [1, .6, .32], (6, 5) = [1, .8, .22], (6, 6) = [1, 1, .12]});
+P2 := Matrix(6, 6, {(1, 1) = [0, 0, 0.1], (1, 2) = [0, .2, -0.1],(1, 3) = [0, .4, 0],  (1, 4) = [0, .6, 0.1], (1, 5) = [0, .8, 0], (1, 6) = [0, 0.9, 0], (2, 1) = [.2, 0, 0], (2, 2) = [.2, .2, .42], (2, 3) = [.2, .4, .42], (2, 4) = [.2, .6, -.2], (2, 5) = [.2, .8, .6], (2, 6) = [.2, 1, 0], (3, 1) = [.4, 0, 0], (3, 2) = [.4, .2, -.2], (3, 3) = [.4, .4, -.2], (3, 4) = [.4, .6, .4], (3, 5) = [.4, .8, -.2], (3, 6) = [.4, 1, 0], (4, 1) = [.6, 0, 0], (4, 2) = [.6, .2, .4], (4, 3) = [.6, .4, .4], (4, 4) = [.6, .6, -.2], (4, 5) = [.6, .8, .41], (4, 6) = [.6, 1, 0], (5, 1) = [.8, 0, 0], (5, 2) = [.8, .2, .4], (5, 3) = [.8, .4, .4], (5, 4) = [.8, .6, -.2], (5, 5) = [.8, .8, .41], (5, 6) = [.8, 1, 0], (6, 1) = [1.1, 0, 0], (6, 2) = [1, .2, 0], (6, 3) = [1, .4, 0], (6, 4) = [1, .6, 0], (6, 5) = [1, .8, 0], (6, 6) = [1, 1, 0.075]});

intersectTest/case20/surfaces.txt

@@ -0,0 +1,3 @@
+> W:=Matrix([[1.01,0.99,1,1.1],[1,1,1,0.98],[1,1,1,0.9],[1,1,1.1,1.1]]);
+> P1:=Matrix(4, 4, {(1, 1) = [.2000000030, 0., 0.], (1, 2) = [.2000000030, .3300000131, 0.], (1, 3) = [.2000000030, .6600000262, 0.], (1, 4) = [.2000000030, 1.000000000, 0.], (2, 1) = [1.200000048, 0., 1.200000048], (2, 2) = [1.200000048, .3300000131, 1.200000048], (2, 3) = [1.200000048, .6600000262, 1.200000048], (2, 4) = [1.200000048, 1.000000000, 1.200000048], (3, 1) = [-.2000000030, 0., 1.200000048], (3, 2) = [-.2000000030, .3300000131, 1.200000048], (3, 3) = [-.2000000030, .6600000262, 1.200000048], (3, 4) = [-.2000000030, 1.000000000, 1.200000048], (4, 1) = [.8000000119, 0., 0.], (4, 2) = [.8000000119, .3300000131, 0.], (4, 3) = [.8000000119, .6600000262, 0.], (4, 4) = [.8000000119, 1.000000000, 0.]});
+> P2:=Matrix(4, 4, {(1, 1) = [.2000000030, .2000000030, .8999999762], (1, 2) = [.2000000030, .3300000131, 1.000000000], (1, 3) = [.2000000030, .6600000262, 1.100000024], (1, 4) = [.2000000030, .8000000119, 1.], (2, 1) = [.3300000131, .2000000030, 0.1999999955e-1], (2, 2) = [.3300000131, .3300000131, 0.1999999955e-1], (2, 3) = [.3300000131, .6600000262, 0.1999999955e-1], (2, 4) = [.3300000131, .8000000119, 0.1999999955e-1], (3, 1) = [.6600000262, .2000000030, 0.1999999955e-1], (3, 2) = [.6600000262, .3300000131, 0.1999999955e-1], (3, 3) = [.6600000262, .6600000262, 0.1999999955e-1], (3, 4) = [.6600000262, .8000000119, 0.1999999955e-1], (4, 1) = [.8000000119, .2000000030, .8000000119], (4, 2) = [.8000000119, .3300000131, 1.000000000], (4, 3) = [.8000000119, .6600000262, 1.100000024], (4, 4) = [.8000000119, .8000000119, 1.]});

intersectTest/case22/surfaces.txt

@@ -0,0 +1,3 @@
+> W:=Matrix([[1,1,1,1,0.9,1,1,1],[1,1,1,1,1,1,0.9,1],[1,1,1,1,1,1,1,1],[1,1.1,1,1,1.1,1,1,1],[1,1,1,1,1,1,1,1],[1,1,1,1,0.9,1,1,1],[1,1,0.9,1,1,1,1,1],[1,0.9,1,1,1.1,1,1,1]]);
+> P1:=Matrix(8, 8, {(1, 1) = [0.,0.,.3], (1, 2) = [0.,.1428571429,.3], (1, 3) = [0.,.2857142857,.3], (1, 4) = [0.,.4285714286,.3], (1, 5) = [0.,.5714285714,.3], (1, 6) = [0.,.7142857143,.3], (1, 7) = [0.,.8571428571,.3], (1, 8) = [0.,1.,.3], (2, 1) = [.1428571429,0.,.3], (2, 2) = [.1428571429,.1428571429,.8], (2, 3) = [.1428571429,.2857142857,.1], (2, 4) = [.1428571429,.4285714286,.6], (2, 5) = [.1428571429,.5714285714,-2.5], (2, 6) = [.1428571429,.7142857143,.6], (2, 7) = [.1428571429,.8571428571,.5], (2, 8) = [.1428571429,1.,.3], (3, 1) = [3.2857142857,0.,.3], (3, 2) = [.2857142857,.1428571429,.8], (3, 3) = [.2857142857,.2857142857,.1], (3, 4) = [.2857142857,.4285714286,.6], (3, 5) = [.2857142857,.5714285714,-2.5], (3, 6) = [.2857142857,.7142857143,.6], (3, 7) = [.2857142857,.8571428571,.5], (3, 8) = [.2857142857,1.,.3], (4, 1) = [.4285714286,0.,.3], (4, 2) = [.4285714286,.1428571429,.1], (4, 3) = [.4285714286,.2857142857,.1], (4, 4) = [.4285714286,.4285714286,.1], (4, 5) = [.4285714286,.5714285714,.1], (4, 6) = [.4285714286,.7142857143,.1], (4, 7) = [.4285714286,.8571428571,.1], (4, 8) = [.4285714286,1.,.3], (5, 1) = [.5714285714,0.,.3], (5, 2) = [.5714285714,.1428571429,.1], (5, 3) = [.5714285714,.2857142857,.1], (5, 4) = [.5714285714,.4285714286,.1], (5, 5) = [.5714285714,.5714285714,.1], (5, 6) = [.5714285714,.7142857143,.1], (5, 7) = [.5714285714,.8571428571,.1], (5, 8) = [.5714285714,1.,.3], (6, 1) = [.7142857143,0.,.3], (6, 2) = [.7142857143,.1428571429,.9], (6, 3) = [.7142857143,.2857142857,.1], (6, 4) = [.7142857143,.4285714286,.77], (6, 5) = [.7142857143,.5714285714,-2.5], (6, 6) = [.7142857143,.7142857143,.77], (6, 7) = [.7142857143,.8571428571,.5], (6, 8) = [.7142857143,1.,.3], (7, 1) = [.8571428571,0.,.3], (7, 2) = [.8571428571,.1428571429,.9], (7, 3) = [.8571428571,.2857142857,.1], (7, 4) = [.8571428571,.4285714286,.77], (7, 5) = [.8571428571,.5714285714,-2.5], (7, 6) = [.8571428571,.7142857143,.77], (7, 7) = [.8571428571,.8571428571,.5], (7, 8) = [.8571428571,1.,.3], (8, 1) = [1.,0.,.3], (8, 2) = [1.,.1428571429,.3], (8, 3) = [1.,.2857142857,.3], (8, 4) = [1.,.4285714286,.3], (8, 5) = [1.,.5714285714,.3], (8, 6) = [1.,.7142857143,.3], (8, 7) = [1.,.8571428571,.3], (8, 8) = [1.,1.,.3]});
+> P2:=Matrix(8, 8, {(1, 1) = [0.,0.,0], (1, 2) = [0.,.1428571429,0.], (1, 3) = [0.,.2857142857,0.], (1, 4) = [0.,.4285714286,0.], (1, 5) = [0.,.5714285714,0.], (1, 6) = [0.,.7142857143,0.], (1, 7) = [0.,.8571428571,0.], (1, 8) = [0.,1.,0.], (2, 1) = [.1428571429,0.,0.], (2, 2) = [.1428571429,.1428571429,.42], (2, 3) = [.1428571429,.2857142857,-.3], (2, 4) = [.1428571429,.4285714286,.42], (2, 5) = [.1428571429,.5714285714,-.2], (2, 6) = [.1428571429,.7142857143,.6], (2, 7) = [.1428571429,.8571428571,.2], (2, 8) = [.1428571429,1.,0.], (3, 1) = [.2857142857,0.,0.], (3, 2) = [.2857142857,.1428571429,.42], (3, 3) = [.2857142857,.2857142857,-.3], (3, 4) = [.2857142857,.4285714286,.42], (3, 5) = [.2857142857,.5714285714,-.2], (3, 6) = [.2857142857,.7142857143,.6], (3, 7) = [.2857142857,.8571428571,.2], (3, 8) = [.2857142857,1.,0.], (4, 1) = [.4285714286,0.,0.], (4, 2) = [.4285714286,.1428571429,.42], (4, 3) = [.4285714286,.2857142857,-.3], (4, 4) = [.4285714286,.4285714286,.42], (4, 5) = [.4285714286,.5714285714,-.2], (4, 6) = [.4285714286,.7142857143,.6], (4, 7) = [.4285714286,.8571428571,.2], (4, 8) = [.4285714286,1.,0.], (5, 1) = [.5714285714,0.,0.], (5, 2) = [.5714285714,.1428571429,.4], (5, 3) = [.5714285714,.2857142857,-.3], (5, 4) = [.5714285714,.4285714286,.4], (5, 5) = [.5714285714,.5714285714,-.2], (5, 6) = [.5714285714,.7142857143,.41], (5, 7) = [.5714285714,.8571428571,.3], (5, 8) = [.5714285714,1.,0.], (6, 1) = [.7142857143,0.,0.], (6, 2) = [.7142857143,.1428571429,.4], (6, 3) = [.7142857143,.2857142857,-.3], (6, 4) = [.7142857143,.4285714286,.4], (6, 5) = [.7142857143,.5714285714,-.2], (6, 6) = [.7142857143,.7142857143,.41], (6, 7) = [.7142857143,.8571428571,.3], (6, 8) = [.7142857143,1.,0.], (7, 1) = [.8571428571,0.,0.], (7, 2) = [.8571428571,.1428571429,.42], (7, 3) = [.8571428571,.2857142857,-.3], (7, 4) = [.8571428571,.4285714286,.42], (7, 5) = [.8571428571,.5714285714,-.2], (7, 6) = [.8571428571,.7142857143,.6], (7, 7) = [.8571428571,.8571428571,.2], (7, 8) = [.8571428571,1.,0.], (8, 1) = [1.,0.,0.], (8, 2) = [1.,.1428571429,0.], (8, 3) = [1.,.2857142857,0.], (8, 4) = [1.,.4285714286,0.], (8, 5) = [1.,.5714285714,0.], (8, 6) = [1.,.7142857143,0.], (8, 7) = [1.,.8571428571,0.], (8, 8) = [1.,1.,0.]});

intersectTest/case3/surfaces.txt

@@ -0,0 +1,3 @@
+W:=Matrix([[1,1,1,1],[1,1,1,1],[1,1,1,1],[1,1,1,1]]);
+P1 := Matrix(4, 4, {(1, 1) = [0, 0, 0], (1, 2) = [0, .33, .25], (1, 3) = [0, .66, 0], (1, 4) = [0, 1, -.25], (2, 1) = [.6, 0, .67], (2, 2) = [.6, .33, .67], (2, 3) = [.6, .66, .67], (2, 4) = [.6, 1, .67], (3, 1) = [.3, 0, .67], (3, 2) = [.3, .33, .67], (3, 3) = [.3, .66, .67], (3, 4) = [.3, 1, .67], (4, 1) = [1, 0, .31], (4, 2) = [1, .33, 0], (4, 3) = [1, .66, .31], (4, 4) = [1, 1, 0]});
+P2 := Matrix(4, 4, {(1, 1) = [0, 0, 1], (1, 2) = [0, .33, 1.25], (1, 3) = [0, .66, 1], (1, 4) = [0, 1, .75], (2, 1) = [.6, 0, .33], (2, 2) = [.6, .33, .33], (2, 3) = [.6, .66, .33], (2, 4) = [.6, 1, .33], (3, 1) = [.3, 0, .33], (3, 2) = [.3, .33, .33], (3, 3) = [.3, .66, .33], (3, 4) = [.3, 1, .33], (4, 1) = [1, 0, 1.31], (4, 2) = [1, .33, 1], (4, 3) = [1, .66, 1.31], (4, 4) = [1, 1, 1]});

intersectTest/case4/surfaces.txt

@@ -0,0 +1,3 @@
+W:=Matrix([[1.01,0.98,1,1],[1,1.02,1,1.03],[1,0.97,1,0.9],[1,1,1,1]]);
+P1 := Matrix(4, 4, {(1, 1) = [0, 0, 0], (1, 2) = [0, .33, .25], (1, 3) = [0, .66, 0], (1, 4) = [0, 1, -.25],(2, 1) = [.6, 0, .67], (2, 2) = [.6, .33, .67], (2, 3) = [.6, .66, .67], (2, 4) = [.6, 1, .67],(3, 1) = [.3, 0, .67], (3, 2) = [.3, .33, .67], (3, 3) = [.3, .66, .67], (3, 4) = [.3, 1, .67],(4, 1) = [1, 0, .31], (4, 2) = [1, .33, 0], (4, 3) = [1, .66, .31], (4, 4) = [1, 1, 0]});
+P2 := Matrix(4, 4, {(1, 1) = [0, 0, 1], (1, 2) = [0, .33, 1.25], (1, 3) = [0, .66, 1], (1, 4) = [0, 1, .75],(2, 1) = [.6, 0, .33], (2, 2) = [.6, .33, .33], (2, 3) = [.6, .66, .33], (2, 4) = [.6, 1, .33],(3, 1) = [.3, 0, .33], (3, 2) = [.3, .33, .33], (3, 3) = [.3, .66, .33], (3, 4) = [.3, 1, .33],(4, 1) = [1, 0, 1.31], (4, 2) = [1, .33, 1], (4, 3) = [1, .66, 1.31], (4, 4) = [1, 1, 1]});

intersectTest/case6/surfaces.txt

@@ -0,0 +1,3 @@
+W:=Matrix([[1.01,0.99,1,1.1],[1,1,1,0.98],[1,1,1,1],[1,1,1,1]]);
+> P1:=Matrix(4, 4, {(1, 1) = [9/16, -45/64, -5/4], (1, 2) = [9/16, -45/64, -5/12], (1, 3) = [9/16, -45/64, 5/12], (1, 4) = [9/16, -45/64, 5/4], (2, 1) = [-73/48, 455/192, -5/4], (2, 2) = [-73/48, 455/192, -5/12], (2, 3) = [-73/48, 455/192, 5/12], (2, 4) = [-73/48, 455/192, 5/4], (3, 1) = [-73/48, -455/192, -5/4], (3, 2) = [-73/48, -455/192, -5/12], (3, 3) = [-73/48, -455/192, 5/12], (3, 4) = [-73/48, -455/192, 5/4], (4, 1) = [9/16, 45/64, -5/4], (4, 2) = [9/16, 45/64, -5/12], (4, 3) = [9/16, 45/64, 5/12], (4, 4) = [9/16, 45/64, 5/4]});
+> P2:=Matrix(4, 4, {(1, 1) = [-1, 2, -5/4], (1, 2) = [-1, 2, -5/12], (1, 3) = [-1, 2, 5/12], (1, 4) = [-1, 2, 5/4], (2, 1) = [-1/3, 0, -5/4], (2, 2) = [-1/3, 0, -5/12], (2, 3) = [-1/3, 0, 5/12], (2, 4) = [-1/3, 0, 5/4], (3, 1) = [1/3, -2/3, -5/4], (3, 2) = [1/3, -2/3, -5/12], (3, 3) = [1/3, -2/3, 5/12], (3, 4) = [1/3, -2/3, 5/4], (4, 1) = [1, 0, -5/4], (4, 2) = [1, 0, -5/12], (4, 3) = [1, 0, 5/12], (4, 4) = [1, 0, 5/4]});

intersectTest/case8/surfaces.txt

@@ -0,0 +1,3 @@
+W:=Matrix([[1.01,0.99,1,1.1],[1,1,1,0.98],[1,1,1,1],[1,1,1,1]]);
+P1:=Matrix(4, 4, {(1, 1) = [1, -1, -1], (1, 2) = [1, -1, -1/3], (1, 3) = [1, -1, 1/3], (1, 4) = [1, -1, 1], (2, 1) = [-1/3, 1, -1], (2, 2) = [-1/3, 1, -1/3], (2, 3) = [-1/3, 1, 1/3], (2, 4) = [-1/3, 1, 1], (3, 1) = [-1/3, -1, -1], (3, 2) = [-1/3, -1, -1/3], (3, 3) = [-1/3, -1, 1/3], (3, 4) = [-1/3, -1, 1], (4, 1) = [1, 1, -1], (4, 2) = [1, 1, -1/3], (4, 3) = [1, 1, 1/3], (4, 4) = [1, 1, 1]});
+P2:=Matrix(4, 4, {(1, 1) = [1, -1, -1], (1, 2) = [1, -1/3, -1], (1, 3) = [1, 1/3, -1], (1, 4) = [1, 1, -1], (2, 1) = [-1/3, -1, 1], (2, 2) = [-1/3, -1/3, 1], (2, 3) = [-1/3, 1/3, 1], (2, 4) = [-1/3, 1, 1], (3, 1) = [-1/3, -1, -1], (3, 2) = [-1/3, -1/3, -1], (3, 3) = [-1/3, 1/3, -1], (3, 4) = [-1/3, 1, -1], (4, 1) = [1, -1, 1], (4, 2) = [1, -1/3, 1], (4, 3) = [1, 1/3, 1], (4, 4) = [1, 1, 1]});

intersectTest/casea1/surfaces.txt

@@ -0,0 +1,3 @@
+W:= Matrix([[1,1,1,1,1],[1,1,1,1,1],[1,1,1,1,1],[1,1,1,1,1],[1,1,1,1,1]]);
+P1 := Matrix(5, 5, {(1, 1) = [0, 0, 1], (1, 2) = [0, 1/4, 1], (1, 3) = [0, 1/2, 1], (1, 4) = [0, 3/4, 1], (1, 5) = [0, 1, 1],(2, 1) = [1/4, 0, 1], (2, 2) = [1/4, 1/4, -730861249/243360000], (2, 3) = [1/4, 1/2, 1], (2, 4) = [1/4, 3/4, -730861249/243360000], (2, 5) = [1/4, 1, 1],(3, 1) = [1/2, 0, 1], (3, 2) = [1/2, 1/4, 1], (3, 3) = [1/2, 1/2, 137671249/136890000], (3, 4) = [1/2, 3/4, 1], (3, 5) = [1/2, 1, 1],(4, 1) = [3/4, 0, 1], (4, 2) = [3/4, 1/4, -730861249/243360000], (4, 3) = [3/4, 1/2, 1], (4, 4) = [3/4, 3/4, -730861249/243360000], (4, 5) = [3/4, 1, 1],(5, 1) = [1, 0, 1], (5, 2) = [1, 1/4, 1], (5, 3) = [1, 1/2, 1], (5, 4) = [1, 3/4, 1], (5, 5) = [1, 1, 1]});
+P2 := Matrix(5, 5, {(1, 1) = [0, 0, -1], (1, 2) = [0, 1/4, -1], (1, 3) = [0, 1/2, -1], (1, 4) = [0, 3/4, -1], (1, 5) = [0, 1, -1],(2, 1) = [1/4, 0, -1], (2, 2) = [1/4, 1/4, 730861249/243360000], (2, 3) = [1/4, 1/2, -1], (2, 4) = [1/4, 3/4, 730861249/243360000], (2, 5) = [1/4, 1, -1],(3, 1) = [1/2, 0, -1], (3, 2) = [1/2, 1/4, -1], (3, 3) = [1/2, 1/2, -137671249/136890000], (3, 4) = [1/2, 3/4, -1], (3, 5) = [1/2, 1, -1],(4, 1) = [3/4, 0, -1], (4, 2) = [3/4, 1/4, 730861249/243360000], (4, 3) = [3/4, 1/2, -1], (4, 4) = [3/4, 3/4, 730861249/243360000], (4, 5) = [3/4, 1, -1],(5, 1) = [1, 0, -1], (5, 2) = [1, 1/4, -1], (5, 3) = [1, 1/2, -1], (5, 4) = [1, 3/4, -1], (5, 5) = [1, 1, -1]});

intersectTest/casea2/curve.txt

@@ -0,0 +1,8 @@
+n_vectors=1
+0.215542 0.215541 -0.000000
+n_vectors=1
+0.215541 0.784458 -0.000000
+n_vectors=1
+0.784458 0.215541 0.000000
+n_vectors=1
+0.784459 0.784458 -0.000000

intersectTest/casea2/surfaces.txt

@@ -0,0 +1,3 @@
+W:= Matrix([[1,1,1,1,1],[1,1,1,1,1],[1,1,1,1,1],[1,1,1,1,1],[1,1,1,1,1]]);
+P1 := Matrix(5, 5, {(1, 1) = [0., 0., 1.], (1, 2) = [0., .2500000000, 1.], (1, 3) = [0., .5000000000, 1.], (1, 4) = [0., .7500000000, 1.], (1, 5) = [0., 1., 1.],(2, 1) = [.2500000000, 0., 1.], (2, 2) = [.2500000000, .2500000000, -5.606324097], (2, 3) = [.2500000000, .5000000000, 1.], (2, 4) = [.2500000000, .7500000000, -5.606324097], (2, 5) = [.2500000000, 1., 1.],(3, 1) = [.5000000000, 0., 1.], (3, 2) = [.5000000000, .2500000000, 1.], (3, 3) = [.5000000000, .5000000000, 12.], (3, 4) = [.5000000000, .7500000000, 1.], (3, 5) = [.5000000000, 1., 1.],(4, 1) = [.7500000000, 0., 1.], (4, 2) = [.7500000000, .2500000000, -5.606324097], (4, 3) = [.7500000000, .5000000000, 1.], (4, 4) = [.7500000000, .7500000000, -5.606324097], (4, 5) = [.7500000000, 1., 1.],(5, 1) = [1., 0., 1.], (5, 2) = [1., .2500000000, 1.], (5, 3) = [1., .5000000000, 1.], (5, 4) = [1., .7500000000, 1.], (5, 5) = [1., 1., 1.]});
+P2 := Matrix(5, 5, {(1, 1) = [0., 0., -1.], (1, 2) = [0., .2500000000, -1.], (1, 3) = [0., .5000000000, -1.], (1, 4) = [0., .7500000000, -1.], (1, 5) = [0., 1., -1.],(2, 1) = [.2500000000, 0., -1.], (2, 2) = [.2500000000, .2500000000, 5.606324097], (2, 3) = [.2500000000, .5000000000, -1.], (2, 4) = [.2500000000, .7500000000, 5.606324097], (2, 5) = [.2500000000, 1., -1.],(3, 1) = [.5000000000, 0., -1.], (3, 2) = [.5000000000, .2500000000, -1.], (3, 3) = [.5000000000, .5000000000, -12.], (3, 4) = [.5000000000, .7500000000, -1.], (3, 5) = [.5000000000, 1., -1.],(4, 1) = [.7500000000, 0., -1.], (4, 2) = [.7500000000, .2500000000, 5.606324097], (4, 3) = [.7500000000, .5000000000, -1.], (4, 4) = [.7500000000, .7500000000, 5.606324097], (4, 5) = [.7500000000, 1., -1.],(5, 1) = [1., 0., -1.], (5, 2) = [1., .2500000000, -1.], (5, 3) = [1., .5000000000, -1.], (5, 4) = [1., .7500000000, -1.], (5, 5) = [1., 1., -1.]});

intersectTest/casea3/surfaces.txt

@@ -0,0 +1,3 @@
+W:= Matrix([[1,1,1,1,1],[1,1,1,1,1],[1,1,1,1,1],[1,1,1,1,1],[1,1,1,1,1]]);
+P1 := Matrix(5, 5, {(1, 1) = [0., 0., 1.], (1, 2) = [0., .2500000000, 1.], (1, 3) = [0., .5000000000, 1.], (1, 4) = [0., .7500000000, 1.], (1, 5) = [0., 1., 1.],(2, 1) = [.2500000000, 0., 1.], (2, 2) = [.2500000000, .2500000000, -5.706324097], (2, 3) = [.2500000000, .5000000000, 1.], (2, 4) = [.2500000000, .7500000000, -5.706324097], (2, 5) = [.2500000000, 1., 1.],(3, 1) = [.5000000000, 0., 1.], (3, 2) = [.5000000000, .2500000000, 1.], (3, 3) = [.5000000000, .5000000000, 12.], (3, 4) = [.5000000000, .7500000000, 1.], (3, 5) = [.5000000000, 1., 1.],(4, 1) = [.7500000000, 0., 1.], (4, 2) = [.7500000000, .2500000000, -5.706324097], (4, 3) = [.7500000000, .5000000000, 1.], (4, 4) = [.7500000000, .7500000000, -5.706324097], (4, 5) = [.7500000000, 1., 1.],(5, 1) = [1., 0., 1.], (5, 2) = [1., .2500000000, 1.], (5, 3) = [1., .5000000000, 1.], (5, 4) = [1., .7500000000, 1.], (5, 5) = [1., 1., 1.]});
+P2 := Matrix(5, 5, {(1, 1) = [0., 0., -1.], (1, 2) = [0., .2500000000, -1.], (1, 3) = [0., .5000000000, -1.], (1, 4) = [0., .7500000000, -1.], (1, 5) = [0., 1., -1.],(2, 1) = [.2500000000, 0., -1.], (2, 2) = [.2500000000, .2500000000, 5.706324097], (2, 3) = [.2500000000, .5000000000, -1.], (2, 4) = [.2500000000, .7500000000, 5.706324097], (2, 5) = [.2500000000, 1., -1.],(3, 1) = [.5000000000, 0., -1.], (3, 2) = [.5000000000, .2500000000, -1.], (3, 3) = [.5000000000, .5000000000, -12.], (3, 4) = [.5000000000, .7500000000, -1.], (3, 5) = [.5000000000, 1., -1.],(4, 1) = [.7500000000, 0., -1.], (4, 2) = [.7500000000, .2500000000, 5.706324097], (4, 3) = [.7500000000, .5000000000, -1.], (4, 4) = [.7500000000, .7500000000, 5.706324097], (4, 5) = [.7500000000, 1., -1.],(5, 1) = [1., 0., -1.], (5, 2) = [1., .2500000000, -1.], (5, 3) = [1., .5000000000, -1.], (5, 4) = [1., .7500000000, -1.], (5, 5) = [1., 1., -1.]});

intersectTest/casea4/surfaces.txt

@@ -0,0 +1,3 @@
+W:= Matrix([[1,1,1,1,1],[1,1,1,1,1],[1,1,1,1,1],[1,1,1,1,1],[1,1,1,1,1]]);
+P1 := Matrix(5, 5, {(1, 1) = [0., 0., 1.], (1, 2) = [0., .2500000000, 1.], (1, 3) = [0., .5000000000, 1.], (1, 4) = [0., .7500000000, 1.], (1, 5) = [0., 1., 1.],(2, 1) = [.2500000000, 0., 1.], (2, 2) = [.2500000000, .2500000000, -5.607324097], (2, 3) = [.2500000000, .5000000000, 1.], (2, 4) = [.2500000000, .7500000000, -5.607324097], (2, 5) = [.2500000000, 1., 1.],(3, 1) = [.5000000000, 0., 1.], (3, 2) = [.5000000000, .2500000000, 1.], (3, 3) = [.5000000000, .5000000000, 12.], (3, 4) = [.5000000000, .7500000000, 1.], (3, 5) = [.5000000000, 1., 1.],(4, 1) = [.7500000000, 0., 1.], (4, 2) = [.7500000000, .2500000000, -5.607324097], (4, 3) = [.7500000000, .5000000000, 1.], (4, 4) = [.7500000000, .7500000000, -5.607324097], (4, 5) = [.7500000000, 1., 1.],(5, 1) = [1., 0., 1.], (5, 2) = [1., .2500000000, 1.], (5, 3) = [1., .5000000000, 1.], (5, 4) = [1., .7500000000, 1.], (5, 5) = [1., 1., 1.]});
+P2 := Matrix(5, 5, {(1, 1) = [0., 0., -1.], (1, 2) = [0., .2500000000, -1.], (1, 3) = [0., .5000000000, -1.], (1, 4) = [0., .7500000000, -1.], (1, 5) = [0., 1., -1.],(2, 1) = [.2500000000, 0., -1.], (2, 2) = [.2500000000, .2500000000, 5.607324097], (2, 3) = [.2500000000, .5000000000, -1.], (2, 4) = [.2500000000, .7500000000, 5.607324097], (2, 5) = [.2500000000, 1., -1.],(3, 1) = [.5000000000, 0., -1.], (3, 2) = [.5000000000, .2500000000, -1.], (3, 3) = [.5000000000, .5000000000, -12.], (3, 4) = [.5000000000, .7500000000, -1.], (3, 5) = [.5000000000, 1., -1.],(4, 1) = [.7500000000, 0., -1.], (4, 2) = [.7500000000, .2500000000, 5.607324097], (4, 3) = [.7500000000, .5000000000, -1.], (4, 4) = [.7500000000, .7500000000, 5.607324097], (4, 5) = [.7500000000, 1., -1.],(5, 1) = [1., 0., -1.], (5, 2) = [1., .2500000000, -1.], (5, 3) = [1., .5000000000, -1.], (5, 4) = [1., .7500000000, -1.], (5, 5) = [1., 1., -1.]});

intersectTest/casea5/curve.txt

@@ -0,0 +1,275 @@
+n_vectors=70
+0.215592 0.215391 0.000000
+0.215606 0.215389 0.000000
+0.215619 0.215389 0.000000
+0.215623 0.215389 -0.000000
+0.215635 0.215391 -0.000000
+0.215646 0.215393 -0.000000
+0.215654 0.215397 -0.000000
+0.215663 0.215402 0.000000
+0.215671 0.215408 0.000000
+0.215672 0.215409 0.000000
+0.215679 0.215417 0.000000
+0.215685 0.215426 0.000000
+0.215685 0.215427 0.000000
+0.215690 0.215437 0.000000
+0.215692 0.215448 0.000000
+0.215694 0.215461 0.000000
+0.215694 0.215475 -0.000000
+0.215693 0.215484 -0.000000
+0.215690 0.215501 -0.000000
+0.215685 0.215518 0.000000
+0.215678 0.215535 -0.000000
+0.215675 0.215539 0.000000
+0.215665 0.215558 0.000000
+0.215654 0.215576 -0.000000
+0.215640 0.215594 -0.000000
+0.215626 0.215611 -0.000000
+0.215622 0.215615 -0.000000
+0.215607 0.215630 0.000000
+0.215590 0.215643 0.000000
+0.215574 0.215656 -0.000000
+0.215556 0.215667 -0.000000
+0.215545 0.215673 0.000000
+0.215529 0.215680 -0.000000
+0.215513 0.215686 0.000000
+0.215497 0.215690 -0.000000
+0.215483 0.215693 -0.000000
+0.215470 0.215694 0.000000
+0.215457 0.215694 0.000000
+0.215446 0.215692 -0.000000
+0.215435 0.215689 -0.000000
+0.215431 0.215687 -0.000000
+0.215422 0.215683 0.000000
+0.215414 0.215677 -0.000000
+0.215412 0.215675 -0.000000
+0.215405 0.215668 0.000000
+0.215399 0.215659 -0.000000
+0.215398 0.215658 -0.000000
+0.215394 0.215647 0.000000
+0.215391 0.215637 0.000000
+0.215389 0.215624 0.000000
+0.215389 0.215610 -0.000000
+0.215390 0.215602 0.000000
+0.215393 0.215586 -0.000000
+0.215397 0.215569 -0.000000
+0.215404 0.215552 -0.000000
+0.215405 0.215548 0.000000
+0.215415 0.215530 0.000000
+0.215426 0.215511 -0.000000
+0.215439 0.215493 -0.000000
+0.215453 0.215476 0.000000
+0.215455 0.215474 0.000000
+0.215471 0.215458 -0.000000
+0.215487 0.215444 -0.000000
+0.215504 0.215431 0.000000
+0.215521 0.215419 0.000000
+0.215529 0.215415 0.000000
+0.215544 0.215407 0.000000
+0.215560 0.215400 -0.000000
+0.215576 0.215395 -0.000000
+0.215592 0.215391 0.000000
+n_vectors=67
+0.215416 0.784321 0.000000
+0.215425 0.784316 0.000000
+0.215435 0.784311 -0.000000
+0.215438 0.784310 0.000000
+0.215450 0.784307 -0.000000
+0.215463 0.784306 0.000000
+0.215469 0.784306 -0.000000
+0.215485 0.784307 0.000000
+0.215500 0.784310 0.000000
+0.215517 0.784315 -0.000000
+0.215535 0.784323 0.000000
+0.215553 0.784332 0.000000
+0.215571 0.784343 0.000000
+0.215586 0.784353 0.000000
+0.215604 0.784368 0.000000
+0.215620 0.784383 0.000000
+0.215635 0.784399 -0.000000
+0.215648 0.784416 0.000000
+0.215658 0.784430 -0.000000
+0.215668 0.784446 -0.000000
+0.215677 0.784463 -0.000000
+0.215683 0.784480 -0.000000
+0.215687 0.784489 0.000000
+0.215691 0.784504 0.000000
+0.215693 0.784518 0.000000
+0.215694 0.784532 0.000000
+0.215693 0.784544 0.000000
+0.215692 0.784555 -0.000000
+0.215690 0.784561 -0.000000
+0.215686 0.784571 0.000000
+0.215681 0.784580 0.000000
+0.215680 0.784582 -0.000000
+0.215673 0.784590 0.000000
+0.215665 0.784597 -0.000000
+0.215663 0.784598 0.000000
+0.215653 0.784604 0.000000
+0.215643 0.784607 0.000000
+0.215631 0.784610 0.000000
+0.215618 0.784611 -0.000000
+0.215611 0.784611 0.000000
+0.215596 0.784609 0.000000
+0.215580 0.784606 0.000000
+0.215563 0.784601 0.000000
+0.215544 0.784593 -0.000000
+0.215526 0.784583 0.000000
+0.215508 0.784572 0.000000
+0.215492 0.784560 -0.000000
+0.215475 0.784545 0.000000
+0.215459 0.784530 0.000000
+0.215444 0.784513 0.000000
+0.215432 0.784498 -0.000000
+0.215421 0.784482 -0.000000
+0.215412 0.784465 -0.000000
+0.215404 0.784448 0.000000
+0.215398 0.784433 0.000000
+0.215393 0.784418 0.000000
+0.215390 0.784404 -0.000000
+0.215389 0.784391 -0.000000
+0.215389 0.784378 -0.000000
+0.215390 0.784366 0.000000
+0.215392 0.784359 0.000000
+0.215395 0.784349 -0.000000
+0.215400 0.784340 -0.000000
+0.215402 0.784337 -0.000000
+0.215408 0.784329 0.000000
+0.215415 0.784322 0.000000
+0.215416 0.784321 0.000000
+n_vectors=66
+0.784321 0.215416 0.000000
+0.784328 0.215409 0.000000
+0.784336 0.215402 0.000000
+0.784337 0.215401 -0.000000
+0.784347 0.215396 -0.000000
+0.784357 0.215392 0.000000
+0.784369 0.215390 0.000000
+0.784382 0.215389 -0.000000
+0.784389 0.215389 -0.000000
+0.784405 0.215391 0.000000
+0.784421 0.215394 0.000000
+0.784438 0.215400 -0.000000
+0.784457 0.215408 0.000000
+0.784475 0.215417 -0.000000
+0.784493 0.215429 -0.000000
+0.784509 0.215441 0.000000
+0.784526 0.215455 -0.000000
+0.784542 0.215471 0.000000
+0.784557 0.215488 0.000000
+0.784569 0.215504 0.000000
+0.784580 0.215520 -0.000000
+0.784589 0.215537 -0.000000
+0.784597 0.215553 0.000000
+0.784602 0.215567 -0.000000
+0.784607 0.215582 0.000000
+0.784610 0.215597 0.000000
+0.784611 0.215610 -0.000000
+0.784611 0.215623 -0.000000
+0.784609 0.215635 -0.000000
+0.784608 0.215641 0.000000
+0.784604 0.215651 -0.000000
+0.784600 0.215661 0.000000
+0.784598 0.215663 0.000000
+0.784592 0.215671 0.000000
+0.784584 0.215678 0.000000
+0.784581 0.215680 -0.000000
+0.784572 0.215686 0.000000
+0.784562 0.215690 0.000000
+0.784550 0.215693 -0.000000
+0.784537 0.215694 -0.000000
+0.784531 0.215694 -0.000000
+0.784516 0.215693 -0.000000
+0.784500 0.215690 -0.000000
+0.784484 0.215685 0.000000
+0.784466 0.215678 0.000000
+0.784448 0.215669 -0.000000
+0.784430 0.215658 -0.000000
+0.784415 0.215647 -0.000000
+0.784397 0.215633 -0.000000
+0.784381 0.215618 0.000000
+0.784366 0.215602 -0.000000
+0.784352 0.215585 -0.000000
+0.784343 0.215571 -0.000000
+0.784333 0.215555 -0.000000
+0.784324 0.215538 -0.000000
+0.784317 0.215521 0.000000
+0.784313 0.215511 0.000000
+0.784309 0.215496 0.000000
+0.784307 0.215482 0.000000
+0.784306 0.215469 -0.000000
+0.784306 0.215457 -0.000000
+0.784308 0.215445 0.000000
+0.784310 0.215439 0.000000
+0.784314 0.215429 -0.000000
+0.784319 0.215420 0.000000
+0.784321 0.215416 0.000000
+n_vectors=68
+0.784610 0.784371 0.000000
+0.784611 0.784385 -0.000000
+0.784610 0.784399 -0.000000
+0.784607 0.784415 -0.000000
+0.784607 0.784418 -0.000000
+0.784601 0.784436 0.000000
+0.784594 0.784453 0.000000
+0.784585 0.784471 -0.000000
+0.784579 0.784482 -0.000000
+0.784566 0.784500 0.000000
+0.784552 0.784518 -0.000000
+0.784537 0.784534 0.000000
+0.784521 0.784550 -0.000000
+0.784511 0.784558 0.000000
+0.784495 0.784570 0.000000
+0.784478 0.784581 0.000000
+0.784461 0.784590 0.000000
+0.784444 0.784598 -0.000000
+0.784440 0.784600 0.000000
+0.784425 0.784605 0.000000
+0.784410 0.784608 -0.000000
+0.784397 0.784610 0.000000
+0.784384 0.784611 -0.000000
+0.784372 0.784610 0.000000
+0.784363 0.784609 -0.000000
+0.784353 0.784606 -0.000000
+0.784343 0.784602 0.000000
+0.784340 0.784600 0.000000
+0.784332 0.784594 0.000000
+0.784324 0.784587 0.000000
+0.784321 0.784584 0.000000
+0.784316 0.784575 0.000000
+0.784311 0.784565 -0.000000
+0.784308 0.784553 0.000000
+0.784306 0.784541 -0.000000
+0.784306 0.784535 0.000000
+0.784307 0.784521 -0.000000
+0.784309 0.784505 -0.000000
+0.784313 0.784491 -0.000000
+0.784319 0.784473 0.000000
+0.784328 0.784455 0.000000
+0.784338 0.784437 0.000000
+0.784346 0.784424 -0.000000
+0.784360 0.784406 0.000000
+0.784374 0.784389 0.000000
+0.784390 0.784373 0.000000
+0.784407 0.784359 -0.000000
+0.784419 0.784349 0.000000
+0.784436 0.784338 -0.000000
+0.784452 0.784329 0.000000
+0.784469 0.784321 0.000000
+0.784484 0.784315 0.000000
+0.784499 0.784311 -0.000000
+0.784514 0.784308 0.000000
+0.784526 0.784306 -0.000000
+0.784539 0.784306 0.000000
+0.784551 0.784307 0.000000
+0.784558 0.784309 0.000000
+0.784568 0.784312 0.000000
+0.784577 0.784317 0.000000
+0.784580 0.784319 -0.000000
+0.784588 0.784325 -0.000000
+0.784595 0.784332 -0.000000
+0.784597 0.784335 0.000000
+0.784603 0.784344 0.000000
+0.784607 0.784355 0.000000
+0.784610 0.784366 0.000000
+0.784610 0.784371 0.000000

intersectTest/casea5/surfaces.txt

@@ -0,0 +1,3 @@
+W:= Matrix([[1,1,1,1,1],[1,1,1,1,1],[1,1,1,1,1],[1,1,1,1,1],[1,1,1,1,1]]);
+P1 := Matrix(5, 5, {(1, 1) = [0., 0., 1.], (1, 2) = [0., .2500000000, 1.], (1, 3) = [0., .5000000000, 1.], (1, 4) = [0., .7500000000, 1.], (1, 5) = [0., 1., 1.],(2, 1) = [.2500000000, 0., 1.], (2, 2) = [.2500000000, .2500000000, -5.606325097], (2, 3) = [.2500000000, .5000000000, 1.], (2, 4) = [.2500000000, .7500000000, -5.606325097], (2, 5) = [.2500000000, 1., 1.],(3, 1) = [.5000000000, 0., 1.], (3, 2) = [.5000000000, .2500000000, 1.], (3, 3) = [.5000000000, .5000000000, 12.], (3, 4) = [.5000000000, .7500000000, 1.], (3, 5) = [.5000000000, 1., 1.],(4, 1) = [.7500000000, 0., 1.], (4, 2) = [.7500000000, .2500000000, -5.606325097], (4, 3) = [.7500000000, .5000000000, 1.], (4, 4) = [.7500000000, .7500000000, -5.606325097], (4, 5) = [.7500000000, 1., 1.],(5, 1) = [1., 0., 1.], (5, 2) = [1., .2500000000, 1.], (5, 3) = [1., .5000000000, 1.], (5, 4) = [1., .7500000000, 1.], (5, 5) = [1., 1., 1.]});
+P2 := Matrix(5, 5, {(1, 1) = [0., 0., -1.], (1, 2) = [0., .2500000000, -1.], (1, 3) = [0., .5000000000, -1.], (1, 4) = [0., .7500000000, -1.], (1, 5) = [0., 1., -1.],(2, 1) = [.2500000000, 0., -1.], (2, 2) = [.2500000000, .2500000000, 5.606325097], (2, 3) = [.2500000000, .5000000000, -1.], (2, 4) = [.2500000000, .7500000000, 5.606325097], (2, 5) = [.2500000000, 1., -1.],(3, 1) = [.5000000000, 0., -1.], (3, 2) = [.5000000000, .2500000000, -1.], (3, 3) = [.5000000000, .5000000000, -12.], (3, 4) = [.5000000000, .7500000000, -1.], (3, 5) = [.5000000000, 1., -1.],(4, 1) = [.7500000000, 0., -1.], (4, 2) = [.7500000000, .2500000000, 5.606325097], (4, 3) = [.7500000000, .5000000000, -1.], (4, 4) = [.7500000000, .7500000000, 5.606325097], (4, 5) = [.7500000000, 1., -1.],(5, 1) = [1., 0., -1.], (5, 2) = [1., .2500000000, -1.], (5, 3) = [1., .5000000000, -1.], (5, 4) = [1., .7500000000, -1.], (5, 5) = [1., 1., -1.]});

main.cpp

@@ -5,6 +5,8 @@
 #include "QList"
 #include "QRegularExpression"
 #include "gauss_map.h"
+#include "bvh.h"
+#include "utils.h"
 
 SrfMesh getMesh(const RationalSurface<float> &s, int sampleLevel) {
     SrfMesh res;
@@ -82,7 +84,7 @@ int main() {
 //    };
 
     ifstream fin;
-    fin.open(R"(E:\qt\OpenGLDemo\surfacePlayer\assets\intersectTest\casea3\surfaces.txt)");
+    fin.open(R"(intersectTest\case2\surfaces.txt)");
 
     string str;
     string tmp;
@@ -90,7 +92,7 @@ int main() {
         getline(fin, tmp);
         str += tmp + "\n";
     }
-    cout << str << endl;
+//    cout << str << endl;
 
     auto infos = QString(str.c_str()).split(";");
     auto wData = infos[0];
@@ -138,38 +140,48 @@ int main() {
     fin.close();
 
     // ====================== 测试 =======================
-    vector <vector<glm::vec3>> s_evaluation;
-    vector <vector<glm::vec3>> f_evaluation;
+    vector<vector<glm::vec3>> s_evaluation;
+    vector<vector<glm::vec3>> f_evaluation;
     // 曲面s和f的切向量。zd*-sznmj
-    vector <vector<glm::vec3>> s_tangent_v;
-    vector <vector<glm::vec3>> f_tangent_u;
-    const vector <vector<glm::vec3>> f_tangent_v;
+    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); // [0, 63]
-    auto mesh2 = getMesh(f, 7);
-
+    const vector<vector<glm::vec3>> s_normal;
+    const vector<vector<glm::vec3>> f_normal;
+
+    auto mesh1 = getMesh(s, 8); // [0, 63]
+    auto mesh2 = getMesh(f, 8);
+
+    BVH bvh1(mesh1.evaluation);
+    BVH bvh2(mesh2.evaluation);
+    bvh1.build();
+    bvh2.build();
+    auto intersectBoxPairs = getOverlapLeafNodes(bvh1, bvh2);
+    printf("box pairs size: %lld\n", intersectBoxPairs.size());
     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;
+//    for (auto pair: intersectBoxPairs) {
+//        singularityJudger.judge({bvh1.tree[pair.first].idx_u, bvh1.tree[pair.first].idx_u + 1},
+//                                {bvh1.tree[pair.first].idx_v, bvh1.tree[pair.first].idx_v + 1},
+//                                {bvh2.tree[pair.second].idx_u, bvh2.tree[pair.second].idx_u + 1},
+//                                {bvh2.tree[pair.second].idx_v, bvh2.tree[pair.second].idx_v + 1});
 //    }
-//    singularityJudger.judge({12, 48}, {12, 48}, {12, 48}, {12, 48});
-    singularityJudger.judge({0, 63}, {0, 63}, {0, 63}, {0, 63});
+
+    double time_cost = utils::get_time();
+//    singularityJudger.judge(intersectBoxPairs, {0, 255}, {0, 255}, {0, 255}, {0, 255});
+    singularityJudger.judge(intersectBoxPairs, {0, 127}, {0, 127}, {0, 127}, {0, 127});
+//    singularityJudger.judge(intersectBoxPairs, {0, 63}, {0, 63}, {0, 63}, {0, 63});
+    time_cost = utils::get_time() - time_cost;
+    printf("time cost: %lf\n", time_cost);
 
     // ================== 测试对整个曲面，gauss能排除多少（或保留多少）==================
-    GaussMap gaussMap1(mesh1.normal);
-    GaussMap gaussMap2(mesh2.normal);
-    gaussMap1.build();
-    gaussMap2.build();
-    auto pairs = getOverlapLeafNodes(gaussMap1, gaussMap2);
-    printf("Gauss Map: keep %lld samples in totally %lld boxes\n", pairs.size(),
-           mesh1.normal.size() * mesh1.normal[0].size() * mesh2.normal.size() * mesh2.normal[0].size());
+//    GaussMap gaussMap1(mesh1.normal);
+//    GaussMap gaussMap2(mesh2.normal);
+//    gaussMap1.build();
+//    gaussMap2.build();
+//    auto pairs = getOverlapLeafNodes(gaussMap1, gaussMap2);
+//    printf("Gauss Map: keep %lld samples in totally %lld boxes\n", pairs.size(),
+//           mesh1.normal.size() * mesh1.normal[0].size() * mesh2.normal.size() * mesh2.normal[0].size());
 
     return 0;
 }

src/SingularityJudger.cpp

@@ -1,19 +1,17 @@
 #include "SingularityJudger.h"
 
 #include <utility>
-#include "gauss_map.h"
 #include "loop_detector.h"
 #include "C2C4.h"
+#include "set"
+#include "utils.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);
-    gaussMap1.build();
-    gaussMap2.build();
+
     if (!isGaussMapsOverlapped(gaussMap1, gaussMap2, focusRange_u1, focusRange_u2, focusRange_v1, focusRange_v2)) {
         // 一定没有交
         return false;
@@ -26,48 +24,120 @@ bool SingularityJudger::judge(pair<int, int> focusRange_u1, pair<int, int> focus
                                     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;
-                printf("(%d, %d) ", focusRange_u1.first + i, focusRange_v1.first + j);
-                // 非零，有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;
-                    printf("singular point: s:(%d,%d), f:(%d,%d)\n", focusRange_u1.first + i, focusRange_v1.first + j, focusRange_u2.first + accordPointOnF.first, focusRange_v2.first + accordPointOnF.second);
-                } else judgeRes[i][j] = 1; // 圆环
-            } else {
-                judgeRes[i][j] = 0; // 0表示啥也没有
-            }
-        }
-    }
+//    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;
+//                printf("(%d, %d) ", focusRange_u1.first + i, focusRange_v1.first + j);
+//                // 非零，有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;
+//                    printf("singular point: s:(%d,%d), f:(%d,%d)\n", focusRange_u1.first + i, focusRange_v1.first + j,
+//                           focusRange_u2.first + accordPointOnF.first, focusRange_v2.first + accordPointOnF.second);
+//                } 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_) {}
+        : srf1(srf1_), srf2(srf2_), mesh1(mesh1_), mesh2(mesh2_),
+          gaussMap1(GaussMap(mesh1_.normal)), gaussMap2(GaussMap(mesh2_.normal)) {
+    gaussMapTimeCost = utils::get_time();
+    gaussMap1.build();
+    gaussMap2.build();
+    gaussMapTimeCost = utils::get_time() - gaussMapTimeCost;
+}
 
 void SingularityJudger::judge2(pair<int, int> focusRange_u1, pair<int, int> focusRange_v1, pair<int, int> focusRange_u2,
                                pair<int, int> focusRange_v2) {
-    GaussMap gaussMap1(mesh1.normal);
-    GaussMap gaussMap2(mesh2.normal);
-    gaussMap1.build();
-    gaussMap2.build();
+
     printf("gauss map: %d\n",
            isGaussMapsOverlapped(gaussMap1, gaussMap2, focusRange_u1, focusRange_u2, focusRange_v1, focusRange_v2));
 
-
 //    gaussMap1.printQuadTree();
 
     C2C4 c2C4(mesh1, mesh2, srf1, srf2);
     printf("c2c4: %d\n",
            c2C4.c2OrC4(focusRange_u1.first, focusRange_v1.first, focusRange_u2.first, focusRange_v2.first).first);
 }
+
+bool hasPair(vector<pair<pair<int, int>, pair<int, int>>> pairs, pair<int, int> p) {
+    for (auto e: pairs) {
+        if (p == e.first)return true;
+    }
+    return false;
+}
+
+
+void SingularityJudger::
+judge(const vector<pair<pair<int, int>, pair<int, int>>> &intersectBoxPairs,
+      pair<int, int> focusRange_u1, pair<int, int> focusRange_v1,
+      pair<int, int> focusRange_u2, pair<int, int> focusRange_v2) {
+    size_t originBoxCnt = intersectBoxPairs.size();
+    int keepCnt = 0;
+    map<pair<int, int>, set<pair<int, int>>> pairMap;
+    double time_record = utils::get_time();
+    for (const auto &boxPair: intersectBoxPairs) {
+        if (isGaussMapsOverlapped(gaussMap1, gaussMap2, {boxPair.first.first, boxPair.first.first},
+                                  {boxPair.second.first, boxPair.second.first},
+                                  {boxPair.first.second, boxPair.first.second},
+                                  {boxPair.second.second, boxPair.second.second})) {
+            pairMap[boxPair.first].insert(boxPair.second);
+            keepCnt++;
+        }
+    }
+    gaussMapTimeCost += utils::get_time() - time_record;
+    printf("time cost of Gauss Map: %lf\n", gaussMapTimeCost);
+    printf("the gauss map kept %d boxes in totally %lld.\n", keepCnt, originBoxCnt);
+//    if (!isGaussMapsOverlapped(gaussMap1, gaussMap2, focusRange_u1, focusRange_u2, focusRange_v1, focusRange_v2)) {
+//        // 一定没有交
+//        return;
+//    }
+    // 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, srf1);
+    loopDetector.detect(intersectBoxPairs, 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) {
+                if (!hasPair(intersectBoxPairs, {focusRange_u1.first + i, focusRange_v1.first + j})) {
+                    // 不在相交box pairs中
+                    continue;
+                }
+                hasLoop = true;
+
+                printf("(%d, %d) ", focusRange_u1.first + i, focusRange_v1.first + j);
+                // 非零，有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;
+//                    printf("singular point: s:(%d,%d), f:(%d,%d)\n", focusRange_u1.first + i, focusRange_v1.first + j,
+//                           focusRange_u2.first + accordPointOnF.first, focusRange_v2.first + accordPointOnF.second);
+//                } else judgeRes[i][j] = 1; // 圆环
+            } else {
+                judgeRes[i][j] = 0; // 0表示啥也没有
+            }
+        }
+    }
+}

src/bvh.cpp

@@ -0,0 +1,88 @@
+#include "bvh.h"
+
+void BVH::recursiveBuild(int level, int idx, int idx_u, int idx_v) {
+    AABB bound;
+
+    tree[idx].level = level;
+    tree[idx].idx_u = idx_u;
+    tree[idx].idx_v = idx_v;
+    
+    if (level == maxLevel) {
+        // 叶子结点
+        // miniaabb
+        for (int i = 0; i < 2; i++)
+            for (int j = 0; j < 2; j++)
+                bound = Union(bound, evaluations[idx_u + i][idx_v + j]);
+
+        tree[idx].bound = bound;
+        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].bound);
+        tree[idx].bound = bound;
+        tree[idx].firstChild = firstChild;
+    }
+}
+
+BVH::BVH(const std::vector<std::vector<glm::vec3>> &evaluations_): evaluations(evaluations_){
+    maxLevel = int(log2(int(evaluations_.size()) - 1) + 1);
+    tree.resize(int((pow(4, maxLevel) - 1)) / 3);
+}
+
+void BVH::build() {
+    recursiveBuild(1, 0, 0, 0);
+}
+
+void recursiveGetOverlapLeafNodes(const BVH &bvh1, const BVH &bvh2, int idx1, int idx2,
+                                  std::vector<std::pair<int, int>> &pairs) {
+    BVHNode A = bvh1.tree[idx1];
+    BVHNode B = bvh2.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.bound, B.bound)) 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(bvh1, bvh2, A.firstChild + i, idx2, pairs);
+    } else if (A.firstChild == -1 && B.firstChild != -1) {
+        // A是叶子结点，B是中间结点
+        for (int i = 0; i < 4; i++) recursiveGetOverlapLeafNodes(bvh1, bvh2, idx1, B.firstChild + i, pairs);
+    } else {
+        // 都是中间结点
+        if (AABBSize(A.bound) > AABBSize(B.bound)) {
+            // A的包围盒更大
+            for (int i = 0; i < 4; i++) recursiveGetOverlapLeafNodes(bvh1, bvh2, A.firstChild + i, idx2, pairs);
+        } else {
+            // B的包围盒更大
+            for (int i = 0; i < 4; i++) recursiveGetOverlapLeafNodes(bvh1, bvh2, idx1, B.firstChild + i, pairs);
+        }
+    }
+}
+
+std::vector<std::pair<std::pair<int, int>, std::pair<int,int>>> getOverlapLeafNodes(const BVH &bvh1, const BVH &bvh2) {
+    std::vector<std::pair<int, int>> resPairs;
+    recursiveGetOverlapLeafNodes(bvh1, bvh2, 0, 0, resPairs);
+    std::vector<std::pair<std::pair<int, int>, std::pair<int,int>>> twoDimensionalPairs(resPairs.size());
+    for(int i = 0; i < resPairs.size(); i++) {
+        twoDimensionalPairs[i] = {{bvh1.tree[resPairs[i].first].idx_u, bvh1.tree[resPairs[i].first].idx_v},
+                                  {bvh2.tree[resPairs[i].second].idx_u, bvh2.tree[resPairs[i].second].idx_v}};
+    }
+    return twoDimensionalPairs;
+}

src/gauss_map.cpp

@@ -13,6 +13,7 @@ GaussMap:: GaussMap(const std::vector<std::vector<glm::vec3>> &normals_): normal
 
 void GaussMap::recursiveBuild(int level, int idx, int idx_u, int idx_v) {
     AABB bound;
+    tree[idx].level = level;
     if (level == maxLevel) {
         // 叶子节点
         for (int i = 0; i < 2; i++)
@@ -23,8 +24,7 @@ void GaussMap::recursiveBuild(int level, int idx, int idx_u, int idx_v) {
                 bound = Union(bound, normal);
             }
 
-        tree[idx].nBound = bound;
-        tree[idx].level = level;
+        tree[idx].bound = bound;
         tree[idx].firstChild = -1;
         return;
     } else {
@@ -35,9 +35,8 @@ void GaussMap::recursiveBuild(int level, int idx, int idx_u, int idx_v) {
         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;
+            bound = Union(bound, tree[firstChild + i].bound);
+        tree[idx].bound = bound;
         tree[idx].firstChild = firstChild;
     }
 }
@@ -73,8 +72,8 @@ 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;
+            auto pMax = tree[idx].bound.pMax;
+            auto pMin = tree[idx].bound.pMin;
             printf("<<%g, %g, %g>,<%g, %g, %g>>  ", pMin.x, pMin.y, pMin.z, pMax.x, pMax.y, pMax.z);
         }
         printf("\n =============$$$$$$$============= \n");
@@ -99,7 +98,7 @@ void recursiveGetOverlapLeafNodes(const GaussMap &gm1, const GaussMap &gm2, int
     };
 
     // 两个包围盒不相交，返回
-    if (!IsOverlap(A.nBound, B.nBound)) return;
+    if (!IsOverlap(A.bound, B.bound)) return;
     // 相交
     if (A.firstChild == -1 && B.firstChild == -1) {
         // 两者都是叶子节点
@@ -112,7 +111,7 @@ void recursiveGetOverlapLeafNodes(const GaussMap &gm1, const GaussMap &gm2, int
         for (int i = 0; i < 4; i++) recursiveGetOverlapLeafNodes(gm1, gm2, idx1, B.firstChild + i, pairs);
     } else {
         // 都是中间结点
-        if (AABBSize(A.nBound) > AABBSize(B.nBound)) {
+        if (AABBSize(A.bound) > AABBSize(B.bound)) {
             // A的包围盒更大
             for (int i = 0; i < 4; i++) recursiveGetOverlapLeafNodes(gm1, gm2, A.firstChild + i, idx2, pairs);
         } else {
@@ -166,7 +165,7 @@ bool isGaussMapsOverlapped(const GaussMap &gm1, const GaussMap &gm2, std::pair<i
         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);
+                                 gm.tree[getStartIdxOfLayerN(commonMaxLayer) + getChildNodeIdx(i, j)].bound);
             }
         }
         return bounding;

src/loop_detector.cpp

@@ -47,6 +47,62 @@ void LoopDetector::initOrientedDistance() {
     }
 }
 
+void LoopDetector::initOrientedDisAndVecFields(const vector<pair<pair<int, int>, pair<int, int>>> &intersectBoxPairs) {
+
+    selectedPointsIdx = vector<vector<pair<int, int>>>(s_subPatchEdgeSampleCnt_u,
+                                                       vector<pair<int, int>>(s_subPatchEdgeSampleCnt_v, {-1, -1}));
+    vectorFields = vector<vector<glm::vec2>>(s_subPatchEdgeSampleCnt_u, vector<glm::vec2>(s_subPatchEdgeSampleCnt_v));
+    int neighborIdxes[4][2] = {{1,  0},
+                               {0,  1},
+                               {-1, 0},
+                               {0,  -1}};
+    for (const auto &boxPair: intersectBoxPairs) {
+        float minDis = FLT_MAX;
+        int minDisFIdx_u = -1, minDisFIdx_v = -1;
+        for (int i = 0; i < 2; i++) {
+            for (int j = 0; j < 2; j++) {
+                for (auto neighborIdx: neighborIdxes) {
+                    int idxI = boxPair.first.first + i + neighborIdx[0];
+                    int idxJ = boxPair.first.second + j + neighborIdx[1];
+                    if (idxI - s_subPatchIdxRange_u.first < 0 ||
+                        idxI - s_subPatchIdxRange_u.first >= s_subPatchEdgeSampleCnt_u ||
+                        idxJ - s_subPatchIdxRange_v.first < 0 ||
+                        idxJ - s_subPatchIdxRange_v.first >= s_subPatchEdgeSampleCnt_v ||
+                        selectedPointsIdx[idxI - s_subPatchIdxRange_u.first][idxJ - s_subPatchIdxRange_v.first] !=
+                        pair<int, int>(-1, -1)) {
+                        // 不重复计算
+                        continue;
+                    }
+                    for (int k = 0; k < f_subPatchEdgeSampleCnt_u; k++) {
+                        for (int l = 0; l < f_subPatchEdgeSampleCnt_v; l++) {
+                            auto dis = glm::distance(
+                                    s_evaluation[idxI][idxJ],
+                                    f_evaluation[f_subPatchIdxRange_u.first + k][f_subPatchIdxRange_v.first + l]);
+                            // 确定f上的对应点时，依据的最小欧氏距离，而不是有向距离
+                            if (dis < minDis) {
+                                minDis = dis;
+                                minDisFIdx_u = k;
+                                minDisFIdx_v = l;
+                            }
+                        }
+                    }
+
+                    selectedPointsIdx[idxI - s_subPatchIdxRange_u.first][idxJ - s_subPatchIdxRange_v.first]
+                            = {minDisFIdx_u, minDisFIdx_v};
+                    // 计算vector fields
+                    auto fPtRelatedIdx = pair<int, int>(minDisFIdx_u, minDisFIdx_v);
+                    auto n2 = f_normal[f_subPatchIdxRange_u.first + fPtRelatedIdx.first][f_subPatchIdxRange_v.first +
+                                                                                         fPtRelatedIdx.second];
+                    auto ps_pu = s_tangent_u[idxI][idxJ];
+                    auto ps_pv = s_tangent_v[idxI][idxJ];
+                    vectorFields[idxI - s_subPatchIdxRange_u.first][idxJ - s_subPatchIdxRange_v.first]
+                            = glm::vec2(glm::dot(n2, ps_pu), glm::dot(n2, ps_pv));
+                }
+            }
+        }
+    }
+}
+
 // Yawei Ma的文章中的做法
 void LoopDetector::initVectorField0() {
     vectorFields = vector<vector<glm::vec2>>(s_subPatchEdgeSampleCnt_u, vector<glm::vec2>(s_subPatchEdgeSampleCnt_v));
@@ -119,12 +175,12 @@ void LoopDetector::getRotationNumber() {
     auto edgeSampleCnt = s_evaluation.size();
     auto uParamCellSize = (*(s.knots_u.end() - 1) - *s.knots_u.begin()) / (float) edgeSampleCnt;
     auto vParamCellSize = (*(s.knots_v.end() - 1) - *s.knots_v.begin()) / (float) edgeSampleCnt;
-    vector<vector<float>> pu_pt = {{0.,  1.},
-                                   {-1., 0.}};
-    vector<vector<float>> pv_pt = {{1., 0.},
-                                   {0., -1.}};
+//    vector<vector<float>> pu_pt = {{0.,  1.},
+//                                   {-1., 0.}};
+//    vector<vector<float>> pv_pt = {{1., 0.},
+//                                   {0., -1.}};
     rotationNumbers = vector<vector<int>>(s_subPatchEdgeSampleCnt_u,
-                                                     vector<int>(s_subPatchEdgeSampleCnt_v));
+                                          vector<int>(s_subPatchEdgeSampleCnt_v));
 
     // 以小格子为单位遍历
     for (int i = 0; i < s_subPatchEdgeSampleCnt_u - 1; i++) {
@@ -163,11 +219,59 @@ void LoopDetector::getRotationNumber() {
             }
             rotationNumbers[i][j] = int(
                     round(((F[0][1] + F[1][1] - F[0][0] - F[1][0]) * uParamCellSize +
-                           (G[0][0] + G[0][1] - G[1][0] - G[1][1]) * vParamCellSize) / 2.));
+                           (G[0][0] + G[0][1] - G[1][0] - G[1][1]) * vParamCellSize)) / 2.);
         }
     }
-    int a = 3;
-    int b = 1;
+//    int a = 3;
+//    int b = 1;
+}
+
+void LoopDetector::getRotationNumber(const vector<pair<pair<int, int>, pair<int, int>>> &intersectBoxPairs) {
+    // vectorFields的参数域是与s（第一个曲面）一致的
+    auto edgeSampleCnt = s_evaluation.size();
+    auto uParamCellSize = (*(s.knots_u.end() - 1) - *s.knots_u.begin()) / (float) edgeSampleCnt;
+    auto vParamCellSize = (*(s.knots_v.end() - 1) - *s.knots_v.begin()) / (float) edgeSampleCnt;
+    rotationNumbers = vector<vector<int>>(s_subPatchEdgeSampleCnt_u,
+                                          vector<int>(s_subPatchEdgeSampleCnt_v));
+    // 以小格子为单位遍历
+    for (const auto &boxPair: intersectBoxPairs) {
+        vector<vector<float>> F(2, vector<float>(2));
+        vector<vector<float>> G(2, vector<float>(2));
+        for (int biasI = 0; biasI < 2; biasI++) {
+            for (int biasJ = 0; biasJ < 2; biasJ++) {
+                auto idxI = boxPair.first.first + biasI - s_subPatchIdxRange_u.first;
+                auto idxJ = boxPair.first.second + biasJ - s_subPatchIdxRange_v.first;
+                auto v = vectorFields[idxI][idxJ];
+                auto vSquareSum = v.x * v.x + v.y * v.y;
+                auto pTheta_pChi = -v.y / vSquareSum;
+                auto pTheta_pPhi = v.x / vSquareSum;
+//                    auto pChi_pu =
+//                            (vectorFields[idxI + 1][idxJ].x - vectorFields[idxI - 1][idxJ].x) / (2 * uParamCellSize);
+//                    auto pChi_pv =
+//                            (vectorFields[idxI][idxJ + 1].x - vectorFields[idxI][idxJ - 1].x) / (2 * vParamCellSize);
+                glm::vec2 pV_pu;
+                glm::vec2 pV_pv;
+                if (idxI == 0)
+                    pV_pu = (vectorFields[idxI + 1][idxJ] - vectorFields[idxI][idxJ]) / uParamCellSize;
+                else if (idxI == s_subPatchEdgeSampleCnt_u - 1)
+                    pV_pu = (vectorFields[idxI][idxJ] - vectorFields[idxI - 1][idxJ]) / uParamCellSize;
+                else
+                    pV_pu = (vectorFields[idxI + 1][idxJ] - vectorFields[idxI - 1][idxJ]) / (2 * uParamCellSize);
+                if (idxJ == 0)
+                    pV_pv = (vectorFields[idxI][idxJ + 1] - vectorFields[idxI][idxJ]) / vParamCellSize;
+                else if (idxJ == s_subPatchEdgeSampleCnt_v - 1)
+                    pV_pv = (vectorFields[idxI][idxJ] - vectorFields[idxI][idxJ - 1]) / vParamCellSize;
+                else
+                    pV_pv = (vectorFields[idxI][idxJ + 1] - vectorFields[idxI][idxJ - 1]) / (2 * vParamCellSize);
+                F[biasI][biasJ] = pTheta_pChi * pV_pu.x + pTheta_pPhi * pV_pu.y;
+                G[biasI][biasJ] = pTheta_pChi * pV_pv.x + pTheta_pPhi * pV_pv.y;
+            }
+        }
+        rotationNumbers[boxPair.first.first - s_subPatchIdxRange_u.first][boxPair.first.second -
+                                                                          s_subPatchIdxRange_v.first] = int(
+                round(((F[0][1] + F[1][1] - F[0][0] - F[1][0]) * uParamCellSize +
+                       (G[0][0] + G[0][1] - G[1][0] - G[1][1]) * vParamCellSize)) / 2.);
+    }
 }
 
 vector<pair<int, int>> LoopDetector::detect(pair<int, int> _s_subPatchIdxRange_u, pair<int, int> _s_subPatchIdxRange_v,
@@ -190,5 +294,25 @@ vector<pair<int, int>> LoopDetector::detect(pair<int, int> _s_subPatchIdxRange_u
     return {};
 }
 
+vector<pair<int, int>> LoopDetector::detect(const vector<pair<pair<int, int>, pair<int, int>>> &intersectBoxPairs,
+                                            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;
+
+    initOrientedDisAndVecFields(intersectBoxPairs);
+    int a = 1;
+    getRotationNumber(intersectBoxPairs);
+    return {};
+}
 
 

src/utils.cpp

@@ -0,0 +1,31 @@
+#include "utils.h"
+
+#if IN_UNIX
+double utils::get_time() {
+    struct timeval tv{};
+    double t;
+
+    gettimeofday(&tv, (struct timezone *) nullptr);
+    t = tv.tv_sec + (double) tv.tv_usec * 1e-6;
+
+    return t;
+}
+#else
+double utils::get_time() {
+    LARGE_INTEGER timer;
+    static LARGE_INTEGER fre;
+    static int init = 0;
+    double t;
+
+    if (init != 1) {
+        QueryPerformanceFrequency(&fre);
+        init = 1;
+    }
+
+    QueryPerformanceCounter(&timer);
+
+    t = timer.QuadPart * 1. / fre.QuadPart;
+
+    return t;
+}
+#endif