From 60670bfd3b2b0a7ae830169e26bf06ed2950a0ad Mon Sep 17 00:00:00 2001 From: Dtouch Date: Mon, 24 Apr 2023 14:20:36 +0800 Subject: [PATCH] fine tune, mark old version functions --- intersectTest/case1/surfaces.txt | 3 + intersectTest/caseb5/curve.txt | 0 intersectTest/caseb5/surfaces.txt | 3 + main.cpp | 251 +++++++++++++++++++++++++++++- src/SingularityJudger.cpp | 45 +++--- src/loop_detector.cpp | 11 +- 6 files changed, 287 insertions(+), 26 deletions(-) create mode 100644 intersectTest/case1/surfaces.txt create mode 100644 intersectTest/caseb5/curve.txt create mode 100644 intersectTest/caseb5/surfaces.txt diff --git a/intersectTest/case1/surfaces.txt b/intersectTest/case1/surfaces.txt new file mode 100644 index 0000000..fa7472b --- /dev/null +++ b/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]}); \ No newline at end of file diff --git a/intersectTest/caseb5/curve.txt b/intersectTest/caseb5/curve.txt new file mode 100644 index 0000000..e69de29 diff --git a/intersectTest/caseb5/surfaces.txt b/intersectTest/caseb5/surfaces.txt new file mode 100644 index 0000000..5708cc2 --- /dev/null +++ b/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.]}); \ No newline at end of file diff --git a/main.cpp b/main.cpp index 3384229..38b3982 100644 --- a/main.cpp +++ b/main.cpp @@ -7,11 +7,165 @@ #include "gauss_map.h" #include "bvh.h" #include "utils.h" +#include "unordered_set" +#include "unordered_map" + +array2 getPtsFromStr(QString srfData); +void printCtrPtsAsQuadruples(const RationalSurface &s); +void printCtrPtsAsQuadruples(const Surface &s); + +void sampleTimeTest(const RationalSurface &s_, int sampleLevel) { + // 由于ParaSolid那边曲面参数为double类型，这里也要保证是double类型（控制变量） + RationalSurface s; + vector knots_u(s_.knots_u.size()); + vector knots_v(s_.knots_v.size()); + array2 weights(s_.weights.rows(), s_.weights.cols()); + array2> 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(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"<(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(endMomScdDer - startMomScdDer).count()); + +} + +void sampleTimeTestNonRational(const RationalSurface &s_, int sampleLevel) { + // 由于ParaSolid那边曲面参数为double类型，这里也要保证是double类型（控制变量） + Surface s; + vector knots_u(s_.knots_u.size()); + vector knots_v(s_.knots_v.size()); + array2> 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(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"<(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(endMomScdDer - startMomScdDer).count()); + +} SrfMesh getMesh(const RationalSurface &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>(res.sampleCntOnSingleDir, vector(res.sampleCntOnSingleDir)); res.tangent_u = vector>(res.sampleCntOnSingleDir, vector(res.sampleCntOnSingleDir)); res.tangent_v = vector>(res.sampleCntOnSingleDir, vector(res.sampleCntOnSingleDir)); @@ -22,26 +176,62 @@ SrfMesh getMesh(const RationalSurface &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"< getPtsFromStr(QString srfData); +int dirs[4][2] = {{-1, 0}, + {0, 1}, + {1, 0}, + {0, -1}}; + +void dfs(const map, set>>& pairMap, +// unordered_map, char> &book, + set>& book, +// vector, pair>>> &boxGroups, + vector>> &boxGroups, + int x, int y, int iOfGroup) { +// book[{x, y}] = 1; +// book.insert(pair, char>(pair(x, y), 1)); + book.insert({x,y}); + boxGroups[iOfGroup].emplace_back(pair(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 s; RationalSurface 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> s_evaluation; vector> f_evaluation; @@ -112,6 +304,9 @@ int main() { const vector> s_normal; const vector> 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, set>> pairMap; +// for (const auto &boxPair: intersectBoxPairs) { +// pairMap[boxPair.first].insert(boxPair.second); +// } +//// vector, pair>>> groups; +// vector>> groups; +//// unordered_map, char> book; +// set>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 getPtsFromStr(QString srfData) { } array2 res = {(size_t) dimU, (size_t) dimV, points}; return res; +} + +void printCtrPtsAsQuadruples(const RationalSurface &s) { + cout< &s) { + cout< focusRange_u1, pair focusRange_v1, pair focusRange_u2, pair focusRange_v2) { // TODO gauss map to exclude patch pair that must not make loop or singular intersection @@ -59,6 +59,7 @@ SingularityJudger(RationalSurface &srf1_, RationalSurface &srf2_, gaussMapTimeCost = utils::get_time() - gaussMapTimeCost; } +// old version void SingularityJudger::judge2(pair focusRange_u1, pair focusRange_v1, pair focusRange_u2, pair focusRange_v2) { @@ -137,28 +138,34 @@ judge(const vector, pair>> &intersectBoxPairs, /** * 合并相邻盒子 */ - 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); - } - 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); } -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(pair(x, y)); diff --git a/src/loop_detector.cpp b/src/loop_detector.cpp index d285159..afc2863 100644 --- a/src/loop_detector.cpp +++ b/src/loop_detector.cpp @@ -14,6 +14,7 @@ LoopDetector(const vector> &s_evaluation_, const vector>>(s_subPatchEdgeSampleCnt_u, @@ -47,8 +48,8 @@ void LoopDetector::initOrientedDistance() { } } +// 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)); @@ -156,7 +157,6 @@ void LoopDetector::initOrientedDisAndVecFields(const map, set, set>(s_subPatchEdgeSampleCnt_u, vector(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>(s_subPatchEdgeSampleCnt_u, vector(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>> &intersectBoxPairs) { // vectorFields的参数域是与s（第一个曲面）一致的 auto edgeSampleCnt = s_evaluation.size(); @@ -382,13 +386,14 @@ void LoopDetector::getRotationNumber(const map, set> LoopDetector::detect(pair _s_subPatchIdxRange_u, pair _s_subPatchIdxRange_v, pair _f_subPatchIdxRange_u, pair _f_subPatchIdxRange_v) {