fine tune, mark old version functions

intersectTest/case1/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,1,1,1,1,1],[1,1,1,1,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/caseb5/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.]});

main.cpp

@@ -7,11 +7,165 @@
 #include "gauss_map.h"
 #include "bvh.h"
 #include "utils.h"
+#include "unordered_set"
+#include "unordered_map"
+
+array2<glm::vec3> getPtsFromStr(QString srfData);
+void printCtrPtsAsQuadruples(const RationalSurface<float> &s);
+void printCtrPtsAsQuadruples(const Surface<double> &s);
+
+void sampleTimeTest(const RationalSurface<float> &s_, int sampleLevel) {
+    // 由于ParaSolid那边曲面参数为double类型，这里也要保证是double类型（控制变量）
+    RationalSurface<double> s;
+    vector<double> knots_u(s_.knots_u.size());
+    vector<double> knots_v(s_.knots_v.size());
+    array2<double> weights(s_.weights.rows(), s_.weights.cols());
+    array2<glm::vec<3, double>> control_points(s_.control_points.rows(), s_.control_points.cols());
+
+    for(int i = 0; i < s_.knots_u.size(); i++) knots_u[i] = (double)s_.knots_u[i];
+    for(int i = 0; i < s_.knots_v.size(); i++) knots_v[i] = (double)s_.knots_v[i];
+    for(int i = 0; i < s_.weights.rows(); i++) {
+        for(int j = 0; j < s_.weights.cols(); j++) {
+            weights(i, j) = (double)s_.weights(i, j);
+        }
+    }
+    for(int i = 0; i < s_.control_points.rows(); i++) {
+        for(int j = 0; j < s_.control_points.cols(); j++) {
+            control_points(i, j).x = (double)s_.control_points(i, j).x;
+            control_points(i, j).y = (double)s_.control_points(i, j).y;
+            control_points(i, j).z = (double)s_.control_points(i, j).z;
+        }
+    }
+
+    s.knots_u = knots_u;
+    s.knots_v = knots_v;
+    s.weights = weights;
+    s.control_points = control_points;
+    s.degree_u = s_.degree_u;
+    s.degree_v = s_.degree_v;
+    auto sampleCnt = int(pow(2, sampleLevel - 1) + 1);
+    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) / double (sampleCnt - 1);
+    auto s_step_v = (*(s.knots_v.end() - 1) - s_first_v) / double (sampleCnt - 1);
+    // 为了分别测试赋值和求梯度的时间，这里把它们分开写了
+
+    auto startMomEval = std::chrono::steady_clock::now();
+    for (int i = 0; i < sampleCnt; i++) {
+        auto u = s_first_u + s_step_u * float(i);
+        for (int j = 0; j < sampleCnt; j++) {
+            auto v = s_first_v + s_step_v * float(j);
+            auto eval = tinynurbs::surfacePoint(s, u, v);
+//            printf("(%d, %d) --> (%g, %g, %g)\n", i, j, eval.x, eval.y, eval.z);
+        }
+    }
+    auto endMomEval = std::chrono::steady_clock::now();
+    printf("time cost of evaluation: %lf ms\n", std::chrono::duration<double, std::milli>(endMomEval - startMomEval).count());
+
+    auto startMomDer = std::chrono::steady_clock::now();
+    for (int i = 0; i < sampleCnt; i++) {
+        auto u = s_first_u + s_step_u * float(i);
+        for (int j = 0; j < sampleCnt; j++) {
+            auto v = s_first_v + s_step_v * float(j);
+            auto der = tinynurbs::surfaceDerivatives(s, 1, u, v);
+//            if(der(0, 0) == res.evaluation[i][j]) cout<<"amazing"<<endl;
+//            else cout<<"what??? ("<<res.evaluation[i][j].x<<" "<<res.evaluation[i][j].y<<" "<<res.evaluation[i][j].z<<") | ("<<der(0, 0).x<<" "<<der(0, 0).y<<" "<<der(0, 0).z<<")"<<endl;
+//            res.tangent_u[i][j] = der(1, 0);
+//            res.tangent_v[i][j] = der(0, 1);
+        }
+    }
+    auto endMomDer = std::chrono::steady_clock::now();
+    printf("time cost of derivatives: %lf ms\n", std::chrono::duration<double, std::milli>(endMomDer - startMomDer).count());
+
+   auto startMomScdDer = std::chrono::steady_clock::now();
+    for (int i = 0; i < sampleCnt; i++) {
+        auto u = s_first_u + s_step_u * float(i);
+        for (int j = 0; j < sampleCnt; j++) {
+            auto v = s_first_v + s_step_v * float(j);
+            auto der = tinynurbs::surfaceDerivatives(s, 2, u, v);
+        }
+    }
+    auto endMomScdDer = std::chrono::steady_clock::now();
+    printf("time cost of second derivatives: %lf ms\n", std::chrono::duration<double, std::milli>(endMomScdDer - startMomScdDer).count());
+
+}
+
+void sampleTimeTestNonRational(const RationalSurface<float> &s_, int sampleLevel) {
+    // 由于ParaSolid那边曲面参数为double类型，这里也要保证是double类型（控制变量）
+    Surface<double> s;
+    vector<double> knots_u(s_.knots_u.size());
+    vector<double> knots_v(s_.knots_v.size());
+    array2<glm::vec<3, double>> control_points(s_.control_points.rows(), s_.control_points.cols());
+
+    for(int i = 0; i < s_.knots_u.size(); i++) knots_u[i] = (double)s_.knots_u[i];
+    for(int i = 0; i < s_.knots_v.size(); i++) knots_v[i] = (double)s_.knots_v[i];
+    for(int i = 0; i < s_.control_points.rows(); i++) {
+        for(int j = 0; j < s_.control_points.cols(); j++) {
+            control_points(i, j).x = (double)s_.control_points(i, j).x;
+            control_points(i, j).y = (double)s_.control_points(i, j).y;
+            control_points(i, j).z = (double)s_.control_points(i, j).z;
+        }
+    }
+
+    s.knots_u = knots_u;
+    s.knots_v = knots_v;
+    s.control_points = control_points;
+    s.degree_u = s_.degree_u;
+    s.degree_v = s_.degree_v;
+    auto sampleCnt = int(pow(2, sampleLevel - 1) + 1);
+    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) / double (sampleCnt - 1);
+    auto s_step_v = (*(s.knots_v.end() - 1) - s_first_v) / double (sampleCnt - 1);
+    // 为了分别测试赋值和求梯度的时间，这里把它们分开写了
+
+    printCtrPtsAsQuadruples(s);
+
+    auto startMomEval = std::chrono::steady_clock::now();
+    for (int i = 0; i < sampleCnt; i++) {
+        auto u = s_first_u + s_step_u * float(i);
+        for (int j = 0; j < sampleCnt; j++) {
+            auto v = s_first_v + s_step_v * float(j);
+            auto eval = tinynurbs::surfacePoint(s, u, v);
+//            res.evaluation[i][j] = tinynurbs::surfacePoint(s, u, v);
+        }
+    }
+    auto endMomEval = std::chrono::steady_clock::now();
+    printf("time cost of evaluation: %lf ms\n", std::chrono::duration<double, std::milli>(endMomEval - startMomEval).count());
+
+    auto startMomDer = std::chrono::steady_clock::now();
+    for (int i = 0; i < sampleCnt; i++) {
+        auto u = s_first_u + s_step_u * float(i);
+        for (int j = 0; j < sampleCnt; j++) {
+            auto v = s_first_v + s_step_v * float(j);
+            auto der = tinynurbs::surfaceDerivatives(s, 1, u, v);
+//            if(der(0, 0) == res.evaluation[i][j]) cout<<"amazing"<<endl;
+//            else cout<<"what??? ("<<res.evaluation[i][j].x<<" "<<res.evaluation[i][j].y<<" "<<res.evaluation[i][j].z<<") | ("<<der(0, 0).x<<" "<<der(0, 0).y<<" "<<der(0, 0).z<<")"<<endl;
+//            res.tangent_u[i][j] = der(1, 0);
+//            res.tangent_v[i][j] = der(0, 1);
+        }
+    }
+    auto endMomDer = std::chrono::steady_clock::now();
+    printf("time cost of derivatives: %lf ms\n", std::chrono::duration<double, std::milli>(endMomDer - startMomDer).count());
+
+
+    auto startMomScdDer = std::chrono::steady_clock::now();
+    for (int i = 0; i < sampleCnt; i++) {
+        auto u = s_first_u + s_step_u * float(i);
+        for (int j = 0; j < sampleCnt; j++) {
+            auto v = s_first_v + s_step_v * float(j);
+            auto der = tinynurbs::surfaceDerivatives(s, 2, u, v);
+        }
+    }
+    auto endMomScdDer = std::chrono::steady_clock::now();
+    printf("time cost of second derivatives: %lf ms\n", std::chrono::duration<double, std::milli>(endMomScdDer - startMomScdDer).count());
+
+}
 
 SrfMesh getMesh(const RationalSurface<float> &s, int sampleLevel) {
     SrfMesh res;
     res.sampleLevel = sampleLevel;
-    res.sampleCntOnSingleDir = int(pow(2, sampleLevel - 1) + 1);
+    auto sampleCnt = 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));
@@ -22,26 +176,62 @@ SrfMesh getMesh(const RationalSurface<float> &s, int sampleLevel) {
     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++) {
+    // evaluation correction test
+    auto u = 0.35f;
+    auto v = 0.35f;
+    auto pt = tinynurbs::surfacePoint(s, u, v);
+    printf("pt on u = %g, v = %g is (%f, %f, %f)\n", u, v, pt.x, pt.y, pt.z);
+    auto der = tinynurbs::surfaceDerivatives(s, 1, u, v);
+    auto du = der(1, 0);
+    auto dv = der(0, 1);
+    printf("derivatives of u when u = %g, v = %g is (%f, %f, %f)\n", u, v, du.x, du.y, du.z);
+    printf("derivatives of v when u = %g, v = %g is (%f, %f, %f)\n", u, v, dv.x, dv.y, dv.z);
+    // end of correction test
+
+    for (int i = 0; i < sampleCnt; i++) {
         auto u = s_first_u + s_step_u * float(i);
-        for (int j = 0; j < res.sampleCntOnSingleDir; j++) {
+        for (int j = 0; j < sampleCnt; 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);
+//            if(der(0, 0) == res.evaluation[i][j]) cout<<"amazing"<<endl;
+//            else cout<<"what??? ("<<res.evaluation[i][j].x<<" "<<res.evaluation[i][j].y<<" "<<res.evaluation[i][j].z<<") | ("<<der(0, 0).x<<" "<<der(0, 0).y<<" "<<der(0, 0).z<<")"<<endl;
             res.tangent_u[i][j] = der(1, 0);
             res.tangent_v[i][j] = der(0, 1);
             res.normal[i][j] = glm::normalize(glm::cross(res.tangent_u[i][j], res.tangent_v[i][j]));
         }
     }
+
     return res;
 }
 
-array2<glm::vec3> getPtsFromStr(QString srfData);
+int dirs[4][2] = {{-1, 0},
+                  {0,  1},
+                  {1,  0},
+                  {0,  -1}};
+
+void dfs(const map<pair<int, int>, set<pair<int, int>>>& pairMap,
+//         unordered_map<pair<int, int>, char> &book,
+         set<pair<int, int>>& book,
+//         vector<vector<pair<pair<int, int>, pair<int, int>>>> &boxGroups,
+         vector<vector<pair<int, int>>> &boxGroups,
+         int x, int y, int iOfGroup) {
+//    book[{x, y}] = 1;
+//    book.insert(pair<pair<int, int>, char>(pair<int, int>(x, y), 1));
+    book.insert({x,y});
+    boxGroups[iOfGroup].emplace_back(pair<int, int>(x, y));
+    for (auto dir: dirs) {
+        auto nx = x + dir[0], ny = y + dir[1];
+        if (book.find({nx, ny}) != book.end() || pairMap.find({nx, ny}) == pairMap.end())continue;
+        dfs(pairMap, book, boxGroups, nx, ny, iOfGroup);
+    }
+}
 
 int main() {
     RationalSurface<float> s;
     RationalSurface<float> f;
 
+
     int level = 6;
     printf("level: %d, sample cnt: %d * %d\n", level, int(pow(2, level - 1)), int(pow(2, level - 1)));
 
@@ -101,6 +291,8 @@ int main() {
 
     fin.close();
 
+    printCtrPtsAsQuadruples(f);
+
     // ====================== 测试 =======================
     vector<vector<glm::vec3>> s_evaluation;
     vector<vector<glm::vec3>> f_evaluation;
@@ -112,6 +304,9 @@ int main() {
     const vector<vector<glm::vec3>> s_normal;
     const vector<vector<glm::vec3>> f_normal;
 
+//    sampleTimeTestNonRational(s, level);
+    sampleTimeTest(s, level);
+
     auto mesh1 = getMesh(s, level);
     auto mesh2 = getMesh(f, level);
 
@@ -121,6 +316,31 @@ int main() {
     bvh2.build();
     auto intersectBoxPairs = getOverlapLeafNodes(bvh1, bvh2); // [{{u1, v1}, {u2, v2}}]
     printf("box pairs size: %lld\n", intersectBoxPairs.size());
+
+    /**
+     * 测试对bvh结果用dfs
+     */
+//    double timeClassify = utils::get_time();
+//    map<pair<int, int>, set<pair<int, int>>> pairMap;
+//    for (const auto &boxPair: intersectBoxPairs) {
+//        pairMap[boxPair.first].insert(boxPair.second);
+//    }
+////    vector<vector<pair<pair<int, int>, pair<int, int>>>> groups;
+//    vector<vector<pair<int, int>>> groups;
+////    unordered_map<pair<int, int>, char> book;
+//    set<pair<int, int>>book;
+//    int groupNum = 0;
+//    for(auto boxPair: intersectBoxPairs) {
+//        if(book.find(boxPair.first) != book.end())continue;
+//        groups.emplace_back();
+//        dfs(pairMap, book, groups, boxPair.first.first, boxPair.first.second, groupNum++);
+//    }
+//    timeClassify = utils::get_time() - timeClassify;
+//    printf("time cost to group boxes from BVH: %lf\n", timeClassify);
+    /**
+     * end of test
+     */
+
     double time_cost = utils::get_time();
     SingularityJudger singularityJudger(s, f, mesh1, mesh2);
 
@@ -177,4 +397,27 @@ array2<glm::vec3> getPtsFromStr(QString srfData) {
     }
     array2<glm::vec3> res = {(size_t) dimU, (size_t) dimV, points};
     return res;
+}
+
+void printCtrPtsAsQuadruples(const RationalSurface<float> &s) {
+    cout<<endl;
+    for(int i = 0; i < s.control_points.rows(); i++) {
+        for(int j = 0; j < s.control_points.cols(); j++) {
+            auto pt = s.control_points(i, j);
+            auto w = s.weights(i, j);
+            cout<<pt.x * w<<", "<<pt.y * w<<", "<<pt.z * w<<", "<<s.weights(i, j)<<", "<<endl;
+        }
+    }
+    cout<<s.control_points.rows() * s.control_points.cols() * 4 << endl;
+}
+
+void printCtrPtsAsQuadruples(const Surface<double> &s) {
+    cout<<endl;
+    for(int i = 0; i < s.control_points.rows(); i++) {
+        for(int j = 0; j < s.control_points.cols(); j++) {
+            auto pt = s.control_points(i, j);
+            cout<<pt.x<<", "<<pt.y<<", "<<pt.z<<", "<<endl;
+        }
+    }
+    cout<<s.control_points.rows() * s.control_points.cols() * 4 << endl;
 }

src/SingularityJudger.cpp

@@ -5,7 +5,7 @@
 #include "set"
 #include "utils.h"
 
-
+// old version
 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
@@ -59,6 +59,7 @@ SingularityJudger(RationalSurface<float> &srf1_, RationalSurface<float> &srf2_,
     gaussMapTimeCost = utils::get_time() - gaussMapTimeCost;
 }
 
+// old version
 void SingularityJudger::judge2(pair<int, int> focusRange_u1, pair<int, int> focusRange_v1, pair<int, int> focusRange_u2,
                                pair<int, int> focusRange_v2) {
 
@@ -137,28 +138,34 @@ judge(const vector<pair<pair<int, int>, pair<int, int>>> &intersectBoxPairs,
     /**
      * 合并相邻盒子
      */
-    vector<vector<pair<int, int>>> cellGroups;
-    set<pair<int, int>> book;
-    int groupNum = 0;
-    for (auto targetCell: loopDetector.targetCells) {
-        if (book.find(targetCell) != book.end()) continue;
-        cellGroups.emplace_back();
-        dfs(book, cellGroups, loopDetector, groupNum++, targetCell.first, targetCell.second);
-    }
-    for(const auto& cellGroup: cellGroups)
-    {
-        printf("{");
-        for(int i = 0; i < cellGroup.size(); i++) {
-            printf("{%d, %d}", cellGroup[i].first, cellGroup[i].second);
-            if(i != cellGroup.size()-1)printf(", ");
-        }
-        printf("},\n");
-    }
+//    double timeClassify = utils::get_time();
+//    vector<vector<pair<int, int>>> cellGroups;
+//    set<pair<int, int>> book;
+//    int groupNum = 0;
+//    for (auto targetCell: loopDetector.targetCells) {
+//        if (book.find(targetCell) != book.end()) continue;
+//        cellGroups.emplace_back();
+//        dfs(book, cellGroups, loopDetector, groupNum++, targetCell.first, targetCell.second);
+//    }
+//    timeClassify = utils::get_time() - timeClassify;
+//    printf("time cost to group boxes: %lf\n", timeClassify);
+//    /**
+//     * 输出盒子分组后结果
+//     */
+//    for(const auto& cellGroup: cellGroups)
+//    {
+//        printf("{");
+//        for(int i = 0; i < cellGroup.size(); i++) {
+//            printf("{%d, %d}", cellGroup[i].first, cellGroup[i].second);
+//            if(i != cellGroup.size()-1)printf(", ");
+//        }
+//        printf("},\n");
+//    }
     // TODO c2 determination for tangency case
 //    C2C4 c2C4(mesh1, mesh2, srf1, srf2);
 }
 
-void SingularityJudger::dfs(set<pair<int, int>> &book, vector<vector<pair<int, int>>> &cellGroups,
+void SingularityJudger:: dfs(set<pair<int, int>> &book, vector<vector<pair<int, int>>> &cellGroups,
                             const LoopDetector &loopDetector, int i, int x, int y) {
     book.insert({x, y});
     cellGroups[i].emplace_back(pair<int, int>(x, y));

src/loop_detector.cpp

@@ -14,6 +14,7 @@ LoopDetector(const vector<vector<glm::vec3>> &s_evaluation_, const vector<vector
           f_tangent_u(f_tangent_u_), f_tangent_v(f_tangent_v_),
           s_normal(s_normal_), f_normal(f_normal_), s(s_) {}
 
+// old version
 void LoopDetector::initOrientedDistance() {
 
     selectedPointsIdx = vector<vector<pair<int, int>>>(s_subPatchEdgeSampleCnt_u,
@@ -47,8 +48,8 @@ void LoopDetector::initOrientedDistance() {
     }
 }
 
+// old version
 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));
@@ -156,7 +157,6 @@ void LoopDetector::initOrientedDisAndVecFields(const map<pair<int, int>, set<pai
                     }
                     selectedPointsIdx[idxI - s_subPatchIdxRange_u.first][idxJ - s_subPatchIdxRange_v.first]
                             = {minDisFIdx_u - f_subPatchIdxRange_u.first, minDisFIdx_v - f_subPatchIdxRange_v.first};
-                    // 计算vector fields
                     auto n2 = f_normal[minDisFIdx_u][minDisFIdx_v];
                     auto ps_pu = s_tangent_u[idxI][idxJ];
                     auto ps_pv = s_tangent_v[idxI][idxJ];
@@ -169,6 +169,7 @@ void LoopDetector::initOrientedDisAndVecFields(const map<pair<int, int>, set<pai
 }
 
 // Yawei Ma的文章中的做法
+// old version
 void LoopDetector::initVectorField0() {
     vectorFields = vector<vector<glm::vec2>>(s_subPatchEdgeSampleCnt_u, vector<glm::vec2>(s_subPatchEdgeSampleCnt_v));
     for (int i = 0; i < s_subPatchEdgeSampleCnt_u; i++) {
@@ -220,6 +221,7 @@ void LoopDetector::initVectorField0() {
 }
 
 // G A Kriezis的文章中的做法
+// old version
 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++) {
@@ -235,6 +237,7 @@ void LoopDetector::initVectorField() {
     }
 }
 
+// old version
 void LoopDetector::getRotationNumber() {
     // vectorFields的参数域是与s（第一个曲面）一致的
     auto edgeSampleCnt = s_evaluation.size();
@@ -290,6 +293,7 @@ void LoopDetector::getRotationNumber() {
 //    int b = 1;
 }
 
+// old version
 void LoopDetector::getRotationNumber(const vector<pair<pair<int, int>, pair<int, int>>> &intersectBoxPairs) {
     // vectorFields的参数域是与s（第一个曲面）一致的
     auto edgeSampleCnt = s_evaluation.size();
@@ -382,13 +386,14 @@ void LoopDetector::getRotationNumber(const map<pair<int, int>, set<pair<int, int
                                         (G[0][0] + G[0][1] - G[1][0] - G[1][1]) * vParamCellSize)) / 2.);
         rotationNumbers[boxPair.first.first - s_subPatchIdxRange_u.first][boxPair.first.second -
                                                                           s_subPatchIdxRange_v.first] = rotationNumber;
-        if(rotationNumber != 0) {
+        if (rotationNumber != 0) {
             targetCells.emplace_back(boxPair.first.first, boxPair.first.second);
         }
 
     }
 }
 
+// old version
 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) {