diff --git a/include/alg.h b/include/alg.h new file mode 100644 index 0000000..a7261e6 --- /dev/null +++ b/include/alg.h @@ -0,0 +1,5 @@ +#pragma once + +class Math { + static void newtonIteration(); +}; diff --git a/include/loop_detector.h b/include/loop_detector.h index be396dd..0776995 100644 --- a/include/loop_detector.h +++ b/include/loop_detector.h @@ -14,10 +14,11 @@ public: const vector> &s_tangent_u_, const vector> &f_tangent_u_, const vector> &s_tangent_v_, const vector> &f_tangent_v_, const vector> &s_normal_, const vector> &f_normal_, - const tinynurbs::RationalSurface &s_); + const tinynurbs::RationalSurface &s_, const tinynurbs::RationalSurface &f_); // 两个曲面 const tinynurbs::RationalSurface &s; + const tinynurbs::RationalSurface &f; // 细分采样的总层数（最高与BVH的总层数相等）。从第二层开始，每一层都表示把原来的一个小sub patch分成了更小的四份 int maxSplitLayer; @@ -36,8 +37,9 @@ public: // 有向距离计算结果。有向距离通过采样的方式计算，其结果与s曲面sub patch中的采样网格规模相同，即与s_evaluation大小一致 // vector> distance; - // 确定有向距离后，对应的f曲面上的最短距离点的位置 + // 确定有向距离后，对应的f曲面上的最短距离点的位置(由于是格点，都是估计值) vector>> selectedPointsIdx; + vector>> selectedPointsParam; // vector fields, 即有向距离关于u、v的导数 vector> vectorFields; vector> rotationNumbers; @@ -71,6 +73,7 @@ private: void initOrientedDistance(); void initOrientedDisAndVecFields(const vector, pair>> &intersectBoxPairs); + void initOrientedDisAndVecFields_deprecated(const map, set>> &intersectBoxPairs); void initOrientedDisAndVecFields(const map, set>> &intersectBoxPairs); // Yawei Ma的文章中的做法 diff --git a/include/utils.h b/include/utils.h index b62d2ed..d93d585 100644 --- a/include/utils.h +++ b/include/utils.h @@ -1,23 +1,17 @@ -#ifndef UNTITLED1_UTILS_H -#define UNTITLED1_UTILS_H +#pragma once +#include +#include +#include +using namespace std; - -#define IN_UNIX 0 // 确定当前运行的操作系统（需要通过系统调用获得时间） namespace utils { -#if IN_UNIX - -#include -#include - -double get_time(); - -#else -#include + class Timer { + public: + explicit Timer(const string& name); + chrono::time_point record; + string eventName; + void restart(); + void end() const; + }; - double get_time(); - -#endif } - - -#endif //UNTITLED1_UTILS_H diff --git a/intersectTest/myCase/surfaces.txt b/intersectTest/singularCase/surfaces.txt similarity index 100% rename from intersectTest/myCase/surfaces.txt rename to intersectTest/singularCase/surfaces.txt diff --git a/intersectTest/singularCase2/surfaces.txt b/intersectTest/singularCase2/surfaces.txt new file mode 100644 index 0000000..75eb323 --- /dev/null +++ b/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]}); diff --git a/intersectTest/singularCase3/surfaces.txt b/intersectTest/singularCase3/surfaces.txt new file mode 100644 index 0000000..696d78d --- /dev/null +++ b/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]}); diff --git a/main.cpp b/main.cpp index 5e11efb..2b4dc58 100644 --- a/main.cpp +++ b/main.cpp @@ -188,11 +188,12 @@ int main() { RationalSurface s; RationalSurface 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 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); diff --git a/src/SingularityJudger.cpp b/src/SingularityJudger.cpp index 3ded92f..d77c032 100644 --- a/src/SingularityJudger.cpp +++ b/src/SingularityJudger.cpp @@ -19,61 +19,60 @@ bool SingularityJudger::judge(pair focusRange_u1, pair focus 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(focusRange_u1, focusRange_v1, focusRange_u2, focusRange_v2); - judgeRes = vector>(loopDetector.s_subPatchEdgeSampleCnt_u - 1, - vector(loopDetector.s_subPatchEdgeSampleCnt_v - 1, 0)); + judgeRes = vector >(loopDetector.s_subPatchEdgeSampleCnt_u - 1, + vector(loopDetector.s_subPatchEdgeSampleCnt_v - 1, 0)); // TODO c2 determination for tangency case bool hasLoop = false; -// C2C4 c2C4(mesh1, mesh2, srf1, srf2); -// for (int i = 0; i < judgeRes.size(); i++) { -// for (int j = 0; j < judgeRes[0].size(); j++) { -// if (loopDetector.rotationNumbers[i][j] != 0) { -// -// hasLoop = true; -// 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 &srf1_, RationalSurface &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(); -// gaussMap1.build(); -// gaussMap2.build(); - gaussMapTimeCost = utils::get_time() - gaussMapTimeCost; + : 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; } // old version void SingularityJudger::judge2(pair focusRange_u1, pair focusRange_v1, pair focusRange_u2, pair focusRange_v2) { - printf("gauss map: %d\n", isGaussMapsOverlapped(gaussMap1, gaussMap2, focusRange_u1, focusRange_u2, focusRange_v1, focusRange_v2)); -// gaussMap1.printQuadTree(); + // 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(const map, set>> &pairs, pair p) { +bool hasPair(const map, set > > &pairs, pair p) { for (const auto &e: pairs) { if (p == e.first)return true; } @@ -82,56 +81,53 @@ bool hasPair(const map, set>> &pairs, pair, int> SingularityJudger:: -judge(const vector, pair>> &intersectBoxPairs, +judge(const vector, pair > > &intersectBoxPairs, pair focusRange_u1, pair focusRange_v1, pair focusRange_u2, pair focusRange_v2) { size_t originBoxCnt = intersectBoxPairs.size(); int keepCnt = 0; - map, set>> pairMap; - double time_record = utils::get_time(); -// auto gmPairMap = getOverlapLeafNodes(gaussMap1, gaussMap2); + map, set > > pairMap; + // 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++; } -// 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++; -// } -// if (isGaussMapsOverlapped(gaussMap1, gaussMap2, {boxPair.first.first, boxPair.first.first}, -// {boxPair.first.second, boxPair.first.second}, -// {boxPair.second.first, boxPair.second.first}, -// {boxPair.second.second, boxPair.second.second})) { -// pairMap[boxPair.first].insert(boxPair.second); -// keepCnt++; -// } -// if (true) { -// pairMap[boxPair.first].insert(boxPair.second); -// keepCnt++; -// } -// if (gmPairMap.find(boxPair.first) != gmPairMap.end() && -// gmPairMap[boxPair.first].find(boxPair.second) != gmPairMap[boxPair.first].end()) { -// pairMap[boxPair.first].insert(boxPair.second); -// keepCnt++; -// } + // 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++; + // } + // if (isGaussMapsOverlapped(gaussMap1, gaussMap2, {boxPair.first.first, boxPair.first.first}, + // {boxPair.first.second, boxPair.first.second}, + // {boxPair.second.first, boxPair.second.first}, + // {boxPair.second.second, boxPair.second.second})) { + // pairMap[boxPair.first].insert(boxPair.second); + // keepCnt++; + // } + // if (true) { + // pairMap[boxPair.first].insert(boxPair.second); + // keepCnt++; + // } + // if (gmPairMap.find(boxPair.first) != gmPairMap.end() && + // gmPairMap[boxPair.first].find(boxPair.second) != gmPairMap[boxPair.first].end()) { + // 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); // 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>(loopDetector.s_subPatchEdgeSampleCnt_u - 1, - vector(loopDetector.s_subPatchEdgeSampleCnt_v - 1, 0)); + judgeRes = vector >(loopDetector.s_subPatchEdgeSampleCnt_u - 1, + vector(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); } @@ -140,35 +136,35 @@ judge(const vector, pair>> &intersectBoxPairs, /** * 合并相邻盒子 */ -// double timeClassify = utils::get_time(); -// vector>> cellGroups; -// set> 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"); -// } + // double timeClassify = utils::get_time(); + // vector>> cellGroups; + // set> 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); + // C2C4 c2C4(mesh1, mesh2, srf1, srf2); return loopDetector.targetCells; } -void SingularityJudger:: dfs(set> &book, vector>> &cellGroups, +void SingularityJudger::dfs(set > &book, vector > > &cellGroups, const LoopDetector &loopDetector, int i, int x, int y) { book.insert({x, y}); cellGroups[i].emplace_back(x, y); @@ -177,4 +173,4 @@ void SingularityJudger:: dfs(set> &book, vector> &s_evaluation_, const vector> &f_evaluation_, - const vector> &s_tangent_u_, const vector> &f_tangent_u_, - const vector> &s_tangent_v_, const vector> &f_tangent_v_, - const vector> &s_normal_, const vector> &f_normal_, - const tinynurbs::RationalSurface &s_) - : 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_) {} +LoopDetector(const vector > &s_evaluation_, const vector > &f_evaluation_, + const vector > &s_tangent_u_, const vector > &f_tangent_u_, + const vector > &s_tangent_v_, const vector > &f_tangent_v_, + const vector > &s_normal_, const vector > &f_normal_, + const tinynurbs::RationalSurface &s_, const tinynurbs::RationalSurface &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_), f(f_) { +} // old version void LoopDetector::initOrientedDistance() { - - selectedPointsIdx = vector>>(s_subPatchEdgeSampleCnt_u, - vector>(s_subPatchEdgeSampleCnt_v)); -// distance = vector>(s_subPatchEdgeSampleCnt_u, vector(s_subPatchEdgeSampleCnt_v)); + selectedPointsIdx = vector > >(s_subPatchEdgeSampleCnt_u, + vector >(s_subPatchEdgeSampleCnt_v)); + // distance = vector>(s_subPatchEdgeSampleCnt_u, vector(s_subPatchEdgeSampleCnt_v)); for (int i = 0; i < s_subPatchEdgeSampleCnt_u; i++) { for (int j = 0; j < s_subPatchEdgeSampleCnt_v; j++) { @@ -28,8 +28,8 @@ void LoopDetector::initOrientedDistance() { for (int k = 0; k < f_subPatchEdgeSampleCnt_u; k++) { for (int l = 0; l < f_subPatchEdgeSampleCnt_v; l++) { auto dis = glm::distance( - s_evaluation[s_subPatchIdxRange_u.first + i][s_subPatchIdxRange_v.first + j], - f_evaluation[f_subPatchIdxRange_u.first + k][f_subPatchIdxRange_v.first + l]); + s_evaluation[s_subPatchIdxRange_u.first + i][s_subPatchIdxRange_v.first + j], + f_evaluation[f_subPatchIdxRange_u.first + k][f_subPatchIdxRange_v.first + l]); // 确定f上的对应点时，依据的最小欧氏距离，而不是有向距离 if (dis < minDis) { minDis = dis; @@ -38,25 +38,28 @@ void LoopDetector::initOrientedDistance() { } } } -// auto N = glm::normalize(s_evaluation[i][j] - f_evaluation[minDisFIdx_u][minDisFIdx_v]); -// // 直接normalize得到的单位向量总是与SF同向，但实际上N应当总是与f的一侧法向量成锐角 -// auto testAngle = N * glm::cross(f_tangent_u[minDisFIdx_u][minDisFIdx_v], f_tangent_v[minDisFIdx_u][minDisFIdx_v]); -// N = testAngle.x + testAngle.y + testAngle.z > 0 ? N : -N; -// distance[i][j] = N * (s_evaluation[i][j] - f_evaluation[minDisFIdx_u][minDisFIdx_v]); + // auto N = glm::normalize(s_evaluation[i][j] - f_evaluation[minDisFIdx_u][minDisFIdx_v]); + // // 直接normalize得到的单位向量总是与SF同向，但实际上N应当总是与f的一侧法向量成锐角 + // auto testAngle = N * glm::cross(f_tangent_u[minDisFIdx_u][minDisFIdx_v], f_tangent_v[minDisFIdx_u][minDisFIdx_v]); + // N = testAngle.x + testAngle.y + testAngle.z > 0 ? N : -N; + // distance[i][j] = N * (s_evaluation[i][j] - f_evaluation[minDisFIdx_u][minDisFIdx_v]); selectedPointsIdx[i][j] = {minDisFIdx_u, minDisFIdx_v}; } } } // old version -void LoopDetector::initOrientedDisAndVecFields(const vector, pair>> &intersectBoxPairs) { - selectedPointsIdx = vector>>(s_subPatchEdgeSampleCnt_u, - vector>(s_subPatchEdgeSampleCnt_v, {-1, -1})); - vectorFields = vector>(s_subPatchEdgeSampleCnt_u, vector(s_subPatchEdgeSampleCnt_v)); - int neighborIdxes[4][2] = {{1, 0}, - {0, 1}, - {-1, 0}, - {0, -1}}; +void LoopDetector::initOrientedDisAndVecFields( + const vector, pair > > &intersectBoxPairs) { + selectedPointsIdx = vector > >(s_subPatchEdgeSampleCnt_u, + vector >(s_subPatchEdgeSampleCnt_v, {-1, -1})); + vectorFields = vector >(s_subPatchEdgeSampleCnt_u, vector(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; @@ -77,8 +80,8 @@ void LoopDetector::initOrientedDisAndVecFields(const vector, 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]); + s_evaluation[idxI][idxJ], + f_evaluation[f_subPatchIdxRange_u.first + k][f_subPatchIdxRange_v.first + l]); // 确定f上的对应点时，依据的最小欧氏距离，而不是有向距离 if (dis < minDis) { minDis = dis; @@ -93,7 +96,7 @@ void LoopDetector::initOrientedDisAndVecFields(const vector, // 计算vector fields auto fPtRelatedIdx = pair(minDisFIdx_u, minDisFIdx_v); auto n2 = f_normal[f_subPatchIdxRange_u.first + fPtRelatedIdx.first][f_subPatchIdxRange_v.first + - fPtRelatedIdx.second]; + 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] @@ -104,17 +107,23 @@ void LoopDetector::initOrientedDisAndVecFields(const vector, } } -void LoopDetector::initOrientedDisAndVecFields(const map, set>> &intersectBoxPairs) { - - selectedPointsIdx = vector>>(s_subPatchEdgeSampleCnt_u, - vector>(s_subPatchEdgeSampleCnt_v, {-1, -1})); - vectorFields = vector>(s_subPatchEdgeSampleCnt_u, vector(s_subPatchEdgeSampleCnt_v)); - 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}}}; // 必须是左上，右上，左下，右下的顺序 +// without newton iteration and with wrong approximate +void LoopDetector::initOrientedDisAndVecFields_deprecated( + const map, set > > &intersectBoxPairs) { + selectedPointsIdx = vector > >(s_subPatchEdgeSampleCnt_u, + vector >(s_subPatchEdgeSampleCnt_v, {-1, -1})); + vectorFields = vector >(s_subPatchEdgeSampleCnt_u, vector(s_subPatchEdgeSampleCnt_v)); + auto vectorFieldsCalculated = vector >(s_subPatchEdgeSampleCnt_u, vector(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}} + }; // 必须是左上，右上，左下，右下的顺序 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, set= f_evaluation.size() || @@ -168,10 +176,76 @@ void LoopDetector::initOrientedDisAndVecFields(const map, set, set > > &intersectBoxPairs) { + selectedPointsIdx = vector(s_subPatchEdgeSampleCnt_u, vector >(s_subPatchEdgeSampleCnt_v, {-1, -1})); + vectorFields = vector(s_subPatchEdgeSampleCnt_u, vector(s_subPatchEdgeSampleCnt_v)); + auto vectorFieldsCalculated = vector(s_subPatchEdgeSampleCnt_u, vector(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(-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 selectedIdx = {minDisFIdx_u - f_subPatchIdxRange_u.first, minDisFIdx_v - f_subPatchIdxRange_v.first}; + pair fURange = {f.knots_u[0], f.knots_u[f.knots_u.size() - 1]}; + pair fVRange = {f.knots_v[0], f.knots_v[f.knots_v.size() - 1]}; + pair selectedParam = { + fURange.first + (fURange.second - fURange.first) / static_cast(f_evaluation.size() - 1) * static_cast(selectedIdx.first), + fVRange.first + (fVRange.second - fVRange.first) / static_cast(f_evaluation[0].size() - 1) * static_cast(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() { - vectorFields = vector>(s_subPatchEdgeSampleCnt_u, vector(s_subPatchEdgeSampleCnt_v)); + vectorFields = vector >(s_subPatchEdgeSampleCnt_u, vector(s_subPatchEdgeSampleCnt_v)); for (int i = 0; i < s_subPatchEdgeSampleCnt_u; i++) { for (int j = 0; j < s_subPatchEdgeSampleCnt_v; j++) { auto fPointIdx = selectedPointsIdx[i][j]; @@ -181,9 +255,9 @@ void LoopDetector::initVectorField0() { if (isnan(N.x) || isnan(N.y) || isnan(N.z)) { // 绝了，曲面s和f的交线正好经过采样点，此时就无法定义向量N // 这里取两个曲面在交点的法向量的平均做近似 -// auto tmp1 = glm::cross(s_tangent_u[i][j], s_tangent_v[i][j]); -// auto tmp2 = glm::cross(f_tangent_u[fPointIdx.first][fPointIdx.second], -// f_tangent_v[fPointIdx.first][fPointIdx.second]); + // auto tmp1 = glm::cross(s_tangent_u[i][j], s_tangent_v[i][j]); + // auto tmp2 = glm::cross(f_tangent_u[fPointIdx.first][fPointIdx.second], + // f_tangent_v[fPointIdx.first][fPointIdx.second]); N = glm::normalize(s_normal[s_subPatchIdxRange_u.first + i][s_subPatchIdxRange_v.first + j] + f_normal[f_subPatchIdxRange_u.first + fPointIdx.first] [f_subPatchIdxRange_v.first + fPointIdx.second]); @@ -197,14 +271,14 @@ void LoopDetector::initVectorField0() { auto u2 = glm::normalize(glm::cross(N, u1)); // s，f曲面在两个方向上的偏导 auto ps_pu = s_tangent_u[s_subPatchIdxRange_u.first + i][s_subPatchIdxRange_v.first + j] - , ps_pv = s_tangent_v[s_subPatchIdxRange_u.first + i][s_subPatchIdxRange_v.first + j]; + , ps_pv = s_tangent_v[s_subPatchIdxRange_u.first + i][s_subPatchIdxRange_v.first + j]; auto pf_pp = f_tangent_u[f_subPatchIdxRange_u.first + fPointIdx.first] - [f_subPatchIdxRange_v.first + fPointIdx.second], pf_pq = f_tangent_v[f_subPatchIdxRange_u.first + - fPointIdx.first] - [f_subPatchIdxRange_v.first + fPointIdx.second]; + [f_subPatchIdxRange_v.first + fPointIdx.second], pf_pq = f_tangent_v[f_subPatchIdxRange_u.first + + fPointIdx.first] + [f_subPatchIdxRange_v.first + fPointIdx.second]; // 构造Aij矩阵，见 APPENDIX (A7) auto inverseAij = glm::inverse( - glm::mat2x2(glm::dot(u1, pf_pp), glm::dot(u1, pf_pq), glm::dot(u2, pf_pp), glm::dot(u2, pf_pq)) + glm::mat2x2(glm::dot(u1, pf_pp), glm::dot(u1, pf_pq), glm::dot(u2, pf_pp), glm::dot(u2, pf_pq)) ); auto mBmn = glm::mat2x2(glm::dot(u1, ps_pu), glm::dot(u1, ps_pv), glm::dot(u2, ps_pu), glm::dot(u2, ps_pv)); auto tmp = glm::vec2(glm::dot(N, pf_pp), glm::dot(N, pf_pq)); @@ -223,13 +297,13 @@ void LoopDetector::initVectorField0() { // G A Kriezis的文章中的做法 // old version void LoopDetector::initVectorField() { - vectorFields = vector>(s_subPatchEdgeSampleCnt_u, vector(s_subPatchEdgeSampleCnt_v)); + vectorFields = vector >(s_subPatchEdgeSampleCnt_u, vector(s_subPatchEdgeSampleCnt_v)); for (int i = 0; i < s_subPatchEdgeSampleCnt_u; i++) { for (int j = 0; j < s_subPatchEdgeSampleCnt_v; j++) { auto fPtRelatedIdx = selectedPointsIdx[i][j]; auto n2 = f_normal[f_subPatchIdxRange_u.first + fPtRelatedIdx.first][f_subPatchIdxRange_v.first + - fPtRelatedIdx.second]; -// glm::vec3 n2(1.,0.,0.); + fPtRelatedIdx.second]; + // glm::vec3 n2(1.,0.,0.); auto ps_pu = s_tangent_u[s_subPatchIdxRange_u.first + i][s_subPatchIdxRange_v.first + j]; auto ps_pv = s_tangent_v[s_subPatchIdxRange_u.first + i][s_subPatchIdxRange_v.first + j]; vectorFields[i][j] = glm::vec2(glm::dot(n2, ps_pu), glm::dot(n2, ps_pv)); @@ -243,12 +317,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> pu_pt = {{0., 1.}, -// {-1., 0.}}; -// vector> pv_pt = {{1., 0.}, -// {0., -1.}}; - rotationNumbers = vector>(s_subPatchEdgeSampleCnt_u, - vector(s_subPatchEdgeSampleCnt_v)); + // vector> pu_pt = {{0., 1.}, + // {-1., 0.}}; + // vector> pv_pt = {{1., 0.}, + // {0., -1.}}; + rotationNumbers = vector >(s_subPatchEdgeSampleCnt_u, + vector(s_subPatchEdgeSampleCnt_v)); // 以小格子为单位遍历 for (int i = 0; i < s_subPatchEdgeSampleCnt_u - 1; i++) { @@ -261,11 +335,11 @@ 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 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); + auto pTheta_pPsi = 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) @@ -280,27 +354,27 @@ 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( - 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.); + 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.); } } -// int a = 3; -// int b = 1; + // int a = 3; + // int b = 1; } // old version -void LoopDetector::getRotationNumber(const vector, pair>> &intersectBoxPairs) { +void LoopDetector::getRotationNumber(const vector, pair > > &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>(s_subPatchEdgeSampleCnt_u, - vector(s_subPatchEdgeSampleCnt_v)); + rotationNumbers = vector >(s_subPatchEdgeSampleCnt_u, + vector(s_subPatchEdgeSampleCnt_v)); // 以小格子为单位遍历 for (const auto &boxPair: intersectBoxPairs) { float F[2][2], G[2][2]; @@ -311,11 +385,11 @@ void LoopDetector::getRotationNumber(const vector, pair, pair, set>> &intersectBoxPairs) { +void LoopDetector::getRotationNumber(const map, set > > &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>(s_subPatchEdgeSampleCnt_u, - vector(s_subPatchEdgeSampleCnt_v)); + rotationNumbers = vector >(s_subPatchEdgeSampleCnt_u, + vector(s_subPatchEdgeSampleCnt_v)); 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,12 +432,16 @@ void LoopDetector::getRotationNumber(const map, set, set> LoopDetector::detect(pair _s_subPatchIdxRange_u, pair _s_subPatchIdxRange_v, - pair _f_subPatchIdxRange_u, - pair _f_subPatchIdxRange_v) { +vector > LoopDetector::detect(pair _s_subPatchIdxRange_u, pair _s_subPatchIdxRange_v, + pair _f_subPatchIdxRange_u, + pair _f_subPatchIdxRange_v) { s_subPatchIdxRange_u = _s_subPatchIdxRange_u; s_subPatchIdxRange_v = _s_subPatchIdxRange_v; f_subPatchIdxRange_u = _f_subPatchIdxRange_u; @@ -418,10 +492,10 @@ vector> LoopDetector::detect(pair _s_subPatchIdxRange_u return {}; } -vector> LoopDetector::detect(const map, set>> &intersectBoxPairs, - pair _s_subPatchIdxRange_u, pair _s_subPatchIdxRange_v, - pair _f_subPatchIdxRange_u, - pair _f_subPatchIdxRange_v) { +vector > LoopDetector::detect(const map, set > > &intersectBoxPairs, + pair _s_subPatchIdxRange_u, pair _s_subPatchIdxRange_v, + pair _f_subPatchIdxRange_u, + pair _f_subPatchIdxRange_v) { s_subPatchIdxRange_u = _s_subPatchIdxRange_u; s_subPatchIdxRange_v = _s_subPatchIdxRange_v; f_subPatchIdxRange_u = _f_subPatchIdxRange_u; @@ -437,5 +511,3 @@ vector> LoopDetector::detect(const map, set(nowTime-record); + cout << "Time cost of even { " << eventName << " } is " << duration.count() << "ms." <