closest point

include/alg.h

@@ -0,0 +1,5 @@
+#pragma once
+
+class Math {
+  static void newtonIteration();
+};

include/loop_detector.h

@@ -14,10 +14,11 @@ public:
                  const vector<vector<glm::vec3>> &s_tangent_u_, const vector<vector<glm::vec3>> &f_tangent_u_,
                  const vector<vector<glm::vec3>> &s_tangent_v_, const vector<vector<glm::vec3>> &f_tangent_v_,
                  const vector<vector<glm::vec3>> &s_normal_, const vector<vector<glm::vec3>> &f_normal_,
-                 const tinynurbs::RationalSurface<float> &s_);
+                 const tinynurbs::RationalSurface<float> &s_, const tinynurbs::RationalSurface<float> &f_);
 
     // 两个曲面
     const tinynurbs::RationalSurface<float> &s;
+    const tinynurbs::RationalSurface<float> &f;
     // 细分采样的总层数（最高与BVH的总层数相等）。从第二层开始，每一层都表示把原来的一个小sub patch分成了更小的四份
     int maxSplitLayer;
 
@@ -36,8 +37,9 @@ public:
 
     // 有向距离计算结果。有向距离通过采样的方式计算，其结果与s曲面sub patch中的采样网格规模相同，即与s_evaluation大小一致
 //    vector<vector<glm::vec3>> distance;
-    // 确定有向距离后，对应的f曲面上的最短距离点的位置
+    // 确定有向距离后，对应的f曲面上的最短距离点的位置(由于是格点，都是估计值)
     vector<vector<pair<int, int>>> selectedPointsIdx;
+    vector<vector<pair<float, float>>> selectedPointsParam;
     // vector fields, 即有向距离关于u、v的导数
     vector<vector<glm::vec2>> vectorFields;
     vector<vector<int>> rotationNumbers;
@@ -71,6 +73,7 @@ private:
     void initOrientedDistance();
 
     void initOrientedDisAndVecFields(const vector<pair<pair<int, int>, pair<int, int>>> &intersectBoxPairs);
+    void initOrientedDisAndVecFields_deprecated(const map<pair<int, int>, set<pair<int, int>>> &intersectBoxPairs);
     void initOrientedDisAndVecFields(const map<pair<int, int>, set<pair<int, int>>> &intersectBoxPairs);
 
     // Yawei Ma的文章中的做法

include/utils.h

@@ -1,23 +1,17 @@
-#ifndef UNTITLED1_UTILS_H
-#define UNTITLED1_UTILS_H
+#pragma once
+#include <chrono>
+#include <string>
+#include <iostream>
+using namespace std;
 
-
-#define IN_UNIX 0  // 确定当前运行的操作系统（需要通过系统调用获得时间）
 namespace utils {
-#if IN_UNIX
-
-#include <sys/time.h>
-#include <ctime>
-
-double get_time();
-
-#else
-#include <windows.h>
+    class Timer {
+    public:
+        explicit Timer(const string& name);
+        chrono::time_point<chrono::steady_clock> record;
+        string eventName;
+        void restart();
+        void end() const;
+    };
 
-    double get_time();
-
-#endif
 }
-
-
-#endif //UNTITLED1_UTILS_H

intersectTest/singularCase2/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]]);
+P1 := Matrix(4, 4, {(1, 1) = [-1.880907654762268, -0.4340604543685913, -0.8796463012695312], (1, 2) = [-1.5740504264831543, 0.11143946647644043, -0.09973067045211792], (1, 3) = [-1.26719331741333, 0.6569393873214722, 0.6801850199699402], (1, 4) = [-0.9603360891342163, 1.202439308166504, 1.4601006507873535], (2, 1) = [-0.9338264465332031, -0.6901867389678955, -1.0731310844421387], (2, 2) = [-0.6269692182540894, -0.14468681812286377, -0.29321545362472534], (2, 3) = [-0.32011204957962036, 0.40081310272216797, 0.48670023679733276], (2, 4) = [-0.01325485110282898, 0.9463130235671997, 1.266615867614746], (3, 1) = [0.01325485110282898, -0.9463130235671997, -1.266615867614746], (3, 2) = [0.32011204957962036, -0.40081310272216797, -0.48670023679733276], (3, 3) = [0.6269692182540894, 0.14468681812286377, 0.29321545362472534], (3, 4) = [0.9338264465332031, 0.6901867389678955, 1.0731310844421387], (4, 1) = [0.9603360891342163, -1.202439308166504, -1.4601006507873535], (4, 2) = [1.26719331741333, -0.6569393873214722, -0.6801850199699402], (4, 3) = [1.5740504264831543, -0.11143946647644043, 0.09973067045211792], (4, 4) = [1.880907654762268, 0.4340604543685913, 0.8796463012695312]});
+P2 := Matrix(4, 4, {(1, 1) = [-1.918592095375061, -0.7532667517662048, -0.6415555477142334], (1, 2) = [-1.5457870960235596, 0.35084420442581177, -0.2782987356185913], (1, 3) = [-1.2389299869537354, 0.8963441252708435, 0.5016169548034668], (1, 4) = [-1.1675983667373657, -0.5531795024871826, 2.769587993621826], (2, 1) = [-0.9715108275413513, -1.0093929767608643, -0.8350403308868408], (2, 2) = [-0.5987059473991394, 0.09471790492534637, -0.47178351879119873], (2, 3) = [-0.291848748922348, 0.6402178406715393, 0.3081321716308594], (2, 4) = [-0.2205171436071396, -0.8093057870864868, 2.5761032104492188], (3, 1) = [-0.02442954108119011, -1.265519380569458, -1.0285251140594482], (3, 2) = [0.3483753502368927, -0.16140837967395782, -0.6652683019638062], (3, 3) = [0.6552324891090393, 0.3840915560722351, 0.11464738845825195], (3, 4) = [0.7265641689300537, -1.065432071685791, 2.3826184272766113], (4, 1) = [0.9226517081260681, -1.5216455459594727, -1.2220098972320557], (4, 2) = [1.2954566478729248, -0.41753464937210083, -0.8587530851364136], (4, 3) = [1.602313756942749, 0.12796525657176971, -0.07883739471435547], (4, 4) = [1.6736453771591187, -1.3215583562850952, 2.189133644104004]});

intersectTest/singularCase3/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]]);
+P1 := Matrix(4, 4, {(1, 1) = [-1.5, -1.5, 0.3999999761581421], (1, 2) = [-1.5, -0.5, -0.29999998211860657], (1, 3) = [-1.5, 0.5, -0.29999998211860657], (1, 4) = [-1.5, 1.5, 1.7505102157592773], (2, 1) = [-0.5, -1.5, 0.3999999761581421], (2, 2) = [-0.5, -0.5, -0.29999998211860657], (2, 3) = [-0.5, 0.5, -0.29999998211860657], (2, 4) = [-0.5, 1.5, 2.1999800205230713], (3, 1) = [0.5, -1.5, 0.3999999761581421], (3, 2) = [0.5, -0.5, -0.29999998211860657], (3, 3) = [0.5, 0.5, -0.29999998211860657], (3, 4) = [0.5, 1.5, 2.1999800205230713], (4, 1) = [1.5, -1.5, 0.3999999761581421], (4, 2) = [1.5, -0.5, -0.29999998211860657], (4, 3) = [1.5, 0.5, -0.29999998211860657], (4, 4) = [1.5, 1.5, 2.1999800205230713]});
+P2 := Matrix(4, 4, {(1, 1) = [-1.5, -1.5, 0.0], (1, 2) = [-1.5, -0.5, 0.0], (1, 3) = [-1.5, 0.5, 0.0], (1, 4) = [-1.5, 1.5, 0.0], (2, 1) = [-0.5, -1.5, 0.0], (2, 2) = [-0.5, -0.5, 0.0], (2, 3) = [-0.5, 0.5, 0.0], (2, 4) = [-0.5, 1.5, 0.0], (3, 1) = [0.5, -1.5, 0.0], (3, 2) = [0.5, -0.5, 0.0], (3, 3) = [0.5, 0.5, 0.0], (3, 4) = [0.5, 1.5, 0.0], (4, 1) = [1.5, -1.5, 0.0], (4, 2) = [1.5, -0.5, 0.0], (4, 3) = [1.5, 0.5, 0.0], (4, 4) = [1.5, 1.5, 0.0]});

main.cpp

@@ -188,11 +188,12 @@ int main() {
     RationalSurface<float> s;
     RationalSurface<float> f;
 
-    int level = 9;
+    int level = 7;
     printf("level: %d, sample cnt: %d * %d\n", level, int(pow(2, level - 1)), int(pow(2, level - 1)));
 
     ifstream fin;
-    fin.open(R"(intersectTest\myCase\surfaces.txt)");
+//    fin.open(R"(intersectTest\singularCase3\surfaces.txt)");
+    fin.open(R"(intersectTest\casea3\surfaces.txt)");
 
     string str;
     string tmp;
@@ -297,15 +298,13 @@ int main() {
      * end of test
      */
 
-    double time_cost = utils::get_time();
+    utils::Timer timerSJ("Winding number detection");
     SingularityJudger singularityJudger(s, f, mesh1, mesh2);
+    timerSJ.end();
 
     pair<int, int> cellIdxFullRange = {0, pow(2, level - 1) - 1};
     auto cellsWithCriticalPts = singularityJudger.judge(intersectBoxPairs, cellIdxFullRange, cellIdxFullRange, cellIdxFullRange, cellIdxFullRange);
 
-    time_cost = utils::get_time() - time_cost;
-    printf("time cost of singularityJudger: %lf\n", time_cost);
-
     // ================== 测试对整个曲面，gauss能排除多少（或保留多少）==================
 //    GaussMap gaussMap1(mesh1.normal);
 //    GaussMap gaussMap2(mesh2.normal);

src/SingularityJudger.cpp

@@ -53,16 +53,15 @@ SingularityJudger::
 SingularityJudger(RationalSurface<float> &srf1_, RationalSurface<float> &srf2_, SrfMesh &mesh1_, SrfMesh &mesh2_)
     : srf1(srf1_), srf2(srf2_), mesh1(mesh1_), mesh2(mesh2_),
       gaussMap1(GaussMap(mesh1_.normal)), gaussMap2(GaussMap(mesh2_.normal)) {
-    gaussMapTimeCost = utils::get_time();
+    // gaussMapTimeCost = utils::get_time();
     //    gaussMap1.build();
     //    gaussMap2.build();
-    gaussMapTimeCost = utils::get_time() - gaussMapTimeCost;
+    // 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) {
-
     printf("gauss map: %d\n",
            isGaussMapsOverlapped(gaussMap1, gaussMap2, focusRange_u1, focusRange_u2, focusRange_v1, focusRange_v2));
 
@@ -88,10 +87,10 @@ judge(const vector<pair<pair<int, int>, pair<int, int>>> &intersectBoxPairs,
     size_t originBoxCnt = intersectBoxPairs.size();
     int keepCnt = 0;
     map<pair<int, int>, set<pair<int, int> > > pairMap;
-    double time_record = utils::get_time();
     //    auto gmPairMap = getOverlapLeafNodes(gaussMap1, gaussMap2);
     for (const auto &boxPair: intersectBoxPairs) {
-        if (isLeafNodesOverlapped(gaussMap1, gaussMap2, boxPair.first, boxPair.second)) {
+        //        if (isLeafNodesOverlapped(gaussMap1, gaussMap2, boxPair.first, boxPair.second)) {
+        if (true) {
             pairMap[boxPair.first].insert(boxPair.second);
             keepCnt++;
         }
@@ -119,19 +118,16 @@ judge(const vector<pair<pair<int, int>, pair<int, int>>> &intersectBoxPairs,
         //            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);
 
     // 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);
+                              mesh2.tangent_v, mesh1.normal, mesh2.normal, srf1, srf2);
     loopDetector.detect(pairMap, 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));
 
     for (auto targetCell: loopDetector.targetCells) {
-//        printf("{%d, %d}, ", targetCell.first, targetCell.second);
         auto uvIdx = targetCell.first;
         printf("At parameter index {%d, %d}, the Poincare index is %d\n", uvIdx.first, uvIdx.second, targetCell.second);
     }

src/alg.cpp

@@ -0,0 +1,5 @@
+#include "alg.h"
+
+void Math::newtonIteration() {
+
+}

src/loop_detector.cpp

@@ -8,15 +8,15 @@ LoopDetector(const vector<vector<glm::vec3>> &s_evaluation_, const vector<vector
              const vector<vector<glm::vec3> > &s_tangent_u_, const vector<vector<glm::vec3> > &f_tangent_u_,
              const vector<vector<glm::vec3> > &s_tangent_v_, const vector<vector<glm::vec3> > &f_tangent_v_,
              const vector<vector<glm::vec3> > &s_normal_, const vector<vector<glm::vec3> > &f_normal_,
-             const tinynurbs::RationalSurface<float> &s_)
+             const tinynurbs::RationalSurface<float> &s_, const tinynurbs::RationalSurface<float> &f_)
     : s_evaluation(s_evaluation_), f_evaluation(f_evaluation_),
       s_tangent_u(s_tangent_u_), s_tangent_v(s_tangent_v_),
       f_tangent_u(f_tangent_u_), f_tangent_v(f_tangent_v_),
-          s_normal(s_normal_), f_normal(f_normal_), s(s_) {}
+      s_normal(s_normal_), f_normal(f_normal_), s(s_), f(f_) {
+}
 
 // old version
 void LoopDetector::initOrientedDistance() {
-
     selectedPointsIdx = vector<vector<pair<int, int> > >(s_subPatchEdgeSampleCnt_u,
                                                          vector<pair<int, int> >(s_subPatchEdgeSampleCnt_v));
     //    distance = vector<vector<glm::vec3>>(s_subPatchEdgeSampleCnt_u, vector<glm::vec3>(s_subPatchEdgeSampleCnt_v));
@@ -49,14 +49,17 @@ void LoopDetector::initOrientedDistance() {
 }
 
 // old version
-void LoopDetector::initOrientedDisAndVecFields(const vector<pair<pair<int, int>, pair<int, int>>> &intersectBoxPairs) {
+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},
+    int neighborIdxes[4][2] = {
+        {1, 0},
         {0, 1},
         {-1, 0},
-                               {0,  -1}};
+        {0, -1}
+    };
     for (const auto &boxPair: intersectBoxPairs) {
         float minDis = FLT_MAX;
         int minDisFIdx_u = -1, minDisFIdx_v = -1;
@@ -104,17 +107,23 @@ void LoopDetector::initOrientedDisAndVecFields(const vector<pair<pair<int, int>,
     }
 }
 
-void LoopDetector::initOrientedDisAndVecFields(const map<pair<int, int>, set<pair<int, int>>> &intersectBoxPairs) {
-
+// without newton iteration and with wrong approximate
+void LoopDetector::initOrientedDisAndVecFields_deprecated(
+    const map<pair<int, int>, set<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},
+    auto vectorFieldsCalculated = vector<vector<int> >(s_subPatchEdgeSampleCnt_u, vector<int>(0));
+    int neighborIdxes[4][2] = {
+        {1, 0},
         {0, 1},
         {-1, 0},
-                               {0,  -1}};
-    int dirs[2][2][2] = {{{-1, 1},  {1, 1}},
-                         {{-1, -1}, {1, -1}}}; // 必须是左上，右上，左下，右下的顺序
+        {0, -1}
+    };
+    int dirs[2][2][2] = {
+        {{-1, 1}, {1, 1}},
+        {{-1, -1}, {1, -1}}
+    }; // 必须是左上，右上，左下，右下的顺序
     for (const auto &boxPair: intersectBoxPairs) {
         for (int i = 0; i < 2; i++) {
             for (int j = 0; j < 2; j++) {
@@ -137,7 +146,6 @@ void LoopDetector::initOrientedDisAndVecFields(const map<pair<int, int>, set<pai
                             for (int l = 0; l < 2; l++) {
                                 for (int step_u = 0; step_u <= 2; step_u++) {
                                     for (int step_v = 0; step_v + step_u <= 2; step_v++) {
-
                                         int fIdx_u = relatedBox.first + k + step_u * dirs[k][l][0];
                                         int fIdx_v = relatedBox.second + l + step_v * dirs[k][l][1];
                                         if (fIdx_u < 0 || fIdx_u >= f_evaluation.size() ||
@@ -168,6 +176,72 @@ void LoopDetector::initOrientedDisAndVecFields(const map<pair<int, int>, set<pai
     }
 }
 
+void LoopDetector::initOrientedDisAndVecFields(const map<pair<int, int>, set<pair<int, int> > > &intersectBoxPairs) {
+    selectedPointsIdx = vector(s_subPatchEdgeSampleCnt_u, vector<pair<int, int> >(s_subPatchEdgeSampleCnt_v, {-1, -1}));
+    vectorFields = vector(s_subPatchEdgeSampleCnt_u, vector<glm::vec2>(s_subPatchEdgeSampleCnt_v));
+    auto vectorFieldsCalculated = vector(s_subPatchEdgeSampleCnt_u, vector<int>(0));
+    int dirs[2][2][2] = {
+        {{-1, -1}, {-1, 1}},
+        {{1, -1}, {1, 1}}
+    }; // 必须是左上，右上，左下，右下的顺序
+    for (const auto &boxPair: intersectBoxPairs) {
+        for (int i = 0; i < 2; i++) {
+            for (int j = 0; j < 2; j++) {
+                float minDis = FLT_MAX;
+                int minDisFIdx_u = -1, minDisFIdx_v = -1;
+                int idxI = boxPair.first.first + i;
+                int idxJ = boxPair.first.second + j;
+                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 (auto relatedBox: boxPair.second) {
+                    for (int k = 0; k < 2; k++) {
+                        for (int l = 0; l < 2; l++) {
+                            for (int step_u = 0; step_u <= 4; step_u++) {
+                                for (int step_v = 0; step_v + step_u <= 4; step_v++) {
+                                    int fIdx_u = relatedBox.first + k + step_u * dirs[k][l][0];
+                                    int fIdx_v = relatedBox.second + l + step_v * dirs[k][l][1];
+                                    if (fIdx_u < 0 || fIdx_u >= f_evaluation.size() ||
+                                        fIdx_v < 0 || fIdx_v >= f_evaluation[0].size())
+                                        continue;
+                                    auto dis = glm::distance(s_evaluation[idxI][idxJ],
+                                                             f_evaluation[fIdx_u][fIdx_v]);
+                                    if (dis < minDis) {
+                                        minDis = dis;
+                                        minDisFIdx_u = fIdx_u;
+                                        minDisFIdx_v = fIdx_v;
+                                    }
+                                }
+                            }
+                        }
+                    }
+                }
+                selectedPointsIdx[idxI - s_subPatchIdxRange_u.first][idxJ - s_subPatchIdxRange_v.first]
+                        = {minDisFIdx_u - f_subPatchIdxRange_u.first, minDisFIdx_v - f_subPatchIdxRange_v.first};
+                pair<int, int> selectedIdx = {minDisFIdx_u - f_subPatchIdxRange_u.first, minDisFIdx_v - f_subPatchIdxRange_v.first};
+                pair<float, float> fURange = {f.knots_u[0], f.knots_u[f.knots_u.size() - 1]};
+                pair<float, float> fVRange = {f.knots_v[0], f.knots_v[f.knots_v.size() - 1]};
+                pair<float, float> selectedParam = {
+                    fURange.first + (fURange.second - fURange.first) / static_cast<float>(f_evaluation.size() - 1) * static_cast<float>(selectedIdx.first),
+                    fVRange.first + (fVRange.second - fVRange.first) / static_cast<float>(f_evaluation[0].size() - 1) * static_cast<float>(selectedIdx.second)
+                };
+                auto n2 = f_normal[minDisFIdx_u][minDisFIdx_v];
+                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的文章中的做法
 // old version
 void LoopDetector::initVectorField0() {
@@ -261,7 +335,7 @@ void LoopDetector::getRotationNumber() {
                     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 pTheta_pPsi = v.x / vSquareSum;
                     //                    auto pChi_pu =
                     //                            (vectorFields[idxI + 1][idxJ].x - vectorFields[idxI - 1][idxJ].x) / (2 * uParamCellSize);
                     //                    auto pChi_pv =
@@ -280,8 +354,8 @@ void LoopDetector::getRotationNumber() {
                         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;
+                    F[biasI][biasJ] = pTheta_pChi * pV_pu.x + pTheta_pPsi * pV_pu.y;
+                    G[biasI][biasJ] = pTheta_pChi * pV_pv.x + pTheta_pPsi * pV_pv.y;
                 }
             }
             rotationNumbers[i][j] = int(
@@ -311,7 +385,7 @@ void LoopDetector::getRotationNumber(const vector<pair<pair<int, int>, pair<int,
                 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 pTheta_pPsi = v.x / vSquareSum;
                 //                    auto pChi_pu =
                 //                            (vectorFields[idxI + 1][idxJ].x - vectorFields[idxI - 1][idxJ].x) / (2 * uParamCellSize);
                 //                    auto pChi_pv =
@@ -330,8 +404,8 @@ void LoopDetector::getRotationNumber(const vector<pair<pair<int, int>, pair<int,
                     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;
+                F[biasI][biasJ] = pTheta_pChi * pV_pu.x + pTheta_pPsi * pV_pu.y;
+                G[biasI][biasJ] = pTheta_pChi * pV_pv.x + pTheta_pPsi * pV_pv.y;
             }
         }
         rotationNumbers[boxPair.first.first - s_subPatchIdxRange_u.first][boxPair.first.second -
@@ -351,9 +425,6 @@ void LoopDetector::getRotationNumber(const map<pair<int, int>, set<pair<int, int
     targetCells = {};
     // 以小格子为单位遍历
     for (const auto &boxPair: intersectBoxPairs) {
-        if (boxPair.first.first == 127 && boxPair.first.second == 0) {
-            int haha = 1;
-        }
         float F[2][2], G[2][2];
         for (int biasI = 0; biasI < 2; biasI++) {
             for (int biasJ = 0; biasJ < 2; biasJ++) {
@@ -361,8 +432,12 @@ void LoopDetector::getRotationNumber(const map<pair<int, int>, set<pair<int, int
                 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;
+                float pTheta_pChi = -1e6;
+                float pTheta_pPsi = 1e6;
+                if (vSquareSum != 0) {
+                    pTheta_pChi = -v.y / vSquareSum;
+                    pTheta_pPsi = v.x / vSquareSum;
+                }
                 //                    auto pChi_pu =
                 //                            (vectorFields[idxI + 1][idxJ].x - vectorFields[idxI - 1][idxJ].x) / (2 * uParamCellSize);
                 //                    auto pChi_pv =
@@ -381,19 +456,18 @@ void LoopDetector::getRotationNumber(const map<pair<int, int>, set<pair<int, int
                     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;
+                F[biasI][biasJ] = pTheta_pChi * pV_pu.x + pTheta_pPsi * pV_pu.y;
+                G[biasI][biasJ] = pTheta_pChi * pV_pv.x + pTheta_pPsi * pV_pv.y;
             }
         }
 
         int rotationNumber = 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.));
         rotationNumbers[boxPair.first.first - s_subPatchIdxRange_u.first][boxPair.first.second -
                                                                           s_subPatchIdxRange_v.first] = rotationNumber;
         if (rotationNumber != 0) {
             targetCells.insert({{boxPair.first.first, boxPair.first.second}, rotationNumber});
         }
-
     }
 }
 
@@ -437,5 +511,3 @@ vector<pair<int, int>> LoopDetector::detect(const map<pair<int, int>, set<pair<i
     getRotationNumber(intersectBoxPairs);
     return {};
 }
-
-

src/utils.cpp

@@ -1,31 +1,16 @@
 #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;
+utils::Timer::Timer(const string &name) {
+    eventName = name;
+    record = chrono::steady_clock::now();
 }
-#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;
+void utils::Timer::restart() {
+    record = chrono::steady_clock::now();
 }
 
-    QueryPerformanceCounter(&timer);
-
-    t = timer.QuadPart * 1. / fre.QuadPart;
-
-    return t;
+void utils::Timer::end() const {
+    auto nowTime = chrono::steady_clock::now();
+    auto duration = chrono::duration_cast<chrono::milliseconds>(nowTime-record);
+    cout << "Time cost of even { " << eventName << " } is " << duration.count() << "ms." <<endl;
 }
-#endif