diff --git a/CMakeLists.txt b/CMakeLists.txt index 8c3a7fe..001d44f 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -7,6 +7,10 @@ include_directories(include) # 引入glm tinynurbs依赖 include_directories(E:/CLib/tinynurbs/include E:/CLib/glm) +find_package(QT NAMES Qt6 Qt5 REQUIRED COMPONENTS Widgets) +find_package(Qt${QT_VERSION_MAJOR} REQUIRED COMPONENTS Widgets) + + add_executable(SingularityJudger main.cpp include/loop_detector.h src/loop_detector.cpp include/gauss_map.h src/gauss_map.cpp @@ -14,3 +18,5 @@ add_executable(SingularityJudger main.cpp include/C2C4.h src/C2C4.cpp include/Range.h src/Range.cpp src/SingularityJudger.cpp include/SingularityJudger.h src/srf_mesh.cpp include/srf_mesh.h) + +target_link_libraries(SingularityJudger PRIVATE Qt${QT_VERSION_MAJOR}::Widgets) \ No newline at end of file diff --git a/include/SingularityJudger.h b/include/SingularityJudger.h index aba42f6..3150068 100644 --- a/include/SingularityJudger.h +++ b/include/SingularityJudger.h @@ -22,9 +22,12 @@ public: RationalSurface &srf2; // 关注mesh的采样网格中的哪个范围 - // 假设mesh.sampleLevel = 5; + // 假设mesh.sampleLevel = 5, 则sampleRange 属于[0, 15]; bool judge(pair focusRange_u1, pair focusRange_v1, pair focusRange_u2, pair focusRange_v2); + // 没有loop detection的版本 + void judge2(pair focusRange_u1, pair focusRange_v1, + pair focusRange_u2, pair focusRange_v2); vector> judgeRes; diff --git a/include/gauss_map.h b/include/gauss_map.h index ab634d4..a5bd12b 100644 --- a/include/gauss_map.h +++ b/include/gauss_map.h @@ -21,7 +21,6 @@ struct GaussMapNode { class GaussMap { public: int maxLevel; - int leafSampleCnt; const std::vector> & normals; tinynurbs::RationalSurface srf; std::vector tree; diff --git a/include/loop_detector.h b/include/loop_detector.h index 749c835..99e2c0c 100644 --- a/include/loop_detector.h +++ b/include/loop_detector.h @@ -11,11 +11,11 @@ public: 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 vector> &s_normal_, const vector> &f_normal_, + const tinynurbs::RationalSurface &s_); // 两个曲面 - tinynurbs::RationalSurface s; - tinynurbs::RationalSurface f; + const tinynurbs::RationalSurface &s; // 细分采样的总层数(最高与BVH的总层数相等)。从第二层开始,每一层都表示把原来的一个小sub patch分成了更小的四份 int maxSplitLayer; @@ -60,11 +60,14 @@ public: private: // 整个曲面一条边上的采样点个数 - int edgeSampleCnt; void initOrientedDistance(); + // Yawei Ma的文章中的做法 + void initVectorField0(); + + // G A Kriezis的文章中的做法 void initVectorField(); void getRotationNumber(); diff --git a/main.cpp b/main.cpp index 2a72e40..304b91c 100644 --- a/main.cpp +++ b/main.cpp @@ -1,7 +1,12 @@ #include #include "SingularityJudger.h" +#include "fstream" +#include "QString" +#include "QList" +#include "QRegularExpression" +#include "gauss_map.h" -SrfMesh getMesh(const RationalSurface& s, int sampleLevel){ +SrfMesh getMesh(const RationalSurface &s, int sampleLevel) { SrfMesh res; res.sampleLevel = sampleLevel; res.sampleCntOnSingleDir = int(pow(2, sampleLevel - 1) + 1); @@ -15,86 +20,193 @@ 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++) { + for (int i = 0; i < res.sampleCntOnSingleDir; i++) { auto u = s_first_u + s_step_u * float(i); - for(int j = 0; j < res.sampleCntOnSingleDir; j++) { + for (int j = 0; j < res.sampleCntOnSingleDir; 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); res.tangent_u[i][j] = der(1, 0); res.tangent_v[i][j] = der(0, 1); - res.normal[i][j] = glm::cross(res.tangent_u[i][j], res.tangent_v[i][j]); + res.normal[i][j] = glm::normalize(glm::cross(res.tangent_u[i][j], res.tangent_v[i][j])); } } + return res; } +array2 getPtsFromStr(QString srfData); int main() { RationalSurface s; RationalSurface f; - s.degree_u = 3; - s.degree_v = 3; - s.knots_u = {0, 0, 0, 0, 1, 1, 1, 1}; - s.knots_v = {0, 0, 0, 0, 1, 1, 1, 1}; - s.control_points = {4, 4, { - glm::vec3(0, 0.3, 0.9), glm::vec3(0, 0.6, 1), glm::vec3(0, 0.9, 1.1), glm::vec3(0, 1.2, 1), - glm::vec3(0.33, 0.3, 0.12), glm::vec3(0.33, 0.6, 0.12), glm::vec3(0.33, 0.9, 0.12), - glm::vec3(0.33, 1.2, 0.12), - glm::vec3(0.66, 0.3, 0.12), glm::vec3(0.66, 0.6, 0.12), glm::vec3(0.66, 0.9, 0.12), - glm::vec3(0.66, 1.2, 0.12), - glm::vec3(1, 0.3, 0.8), glm::vec3(1, 0.6, 1), glm::vec3(1, 0.9, 1.1), glm::vec3(1, 1.2, 1) - }}; - s.weights = {4, 4, - { - 1, 1, 1, 1, - 1, 1, 1, 1, - 1, 1, 1, 1, - 1, 1, 1, 1, - } - }; - - f.degree_u = 3; - f.degree_v = 3; - f.knots_u = {0, 0, 0, 0, 1, 1, 1, 1}; - f.knots_v = {0, 0, 0, 0, 1, 1, 1, 1}; - f.control_points = {4, 4, { - glm::vec3(0, 0.2, 0.9), glm::vec3(0, 0.5, 1.8), glm::vec3(0, 0.8, 1.1), glm::vec3(0, 1.2, 1), - glm::vec3(0.33, 0.2, 0.12), glm::vec3(0.33, 0.5, 0.42), glm::vec3(0.33, 0.9, -0.62), - glm::vec3(0.33, 1.1, -1.756), - glm::vec3(0.66, 0.2, 0.12), glm::vec3(0.66, 0.5, 0.42), glm::vec3(0.66, 0.9, -0.62), - glm::vec3(0.66, 1.0, -1.756), - glm::vec3(1, 0.2, 0.8), glm::vec3(1, 0.5, 1), glm::vec3(1, 0.9, 1.1), glm::vec3(1, 1.2, 1) - }}; - f.weights = {4, 4, - { - 1, 1, 1, 1, - 1, 1, 1, 1, - 1, 1, 1, 1, - 1, 1, 1, 1, - } - }; - - vector> s_evaluation; - vector> f_evaluation; - // 曲面s和f的切向量。 - vector> s_tangent_u; - vector> s_tangent_v; - vector> f_tangent_u; - const vector> f_tangent_v; +// s.degree_u = 3; +// s.degree_v = 3; +// s.knots_u = {0, 0, 0, 0, 1, 1, 1, 1}; +// s.knots_v = {0, 0, 0, 0, 1, 1, 1, 1}; +// s.control_points = {4, 4, { +// glm::vec3(0, 0.3, 0.9), glm::vec3(0, 0.6, 1), glm::vec3(0, 0.9, 1.1), glm::vec3(0, 1.2, 1), +// glm::vec3(0.33, 0.3, 0.12), glm::vec3(0.33, 0.6, 0.12), glm::vec3(0.33, 0.9, 0.12), +// glm::vec3(0.33, 1.2, 0.12), +// glm::vec3(0.66, 0.3, 0.12), glm::vec3(0.66, 0.6, 0.12), glm::vec3(0.66, 0.9, 0.12), +// glm::vec3(0.66, 1.2, 0.12), +// glm::vec3(1, 0.3, 0.8), glm::vec3(1, 0.6, 1), glm::vec3(1, 0.9, 1.1), glm::vec3(1, 1.2, 1) +// }}; +// s.weights = {4, 4, +// { +// 1, 1, 1, 1, +// 1, 1, 1, 1, +// 1, 1, 1, 1, +// 1, 1, 1, 1, +// } +// }; +// +// f.degree_u = 3; +// f.degree_v = 3; +// f.knots_u = {0, 0, 0, 0, 1, 1, 1, 1}; +// f.knots_v = {0, 0, 0, 0, 1, 1, 1, 1}; +// f.control_points = {4, 4, { +// glm::vec3(0, 0.2, 0.9), glm::vec3(0, 0.5, 1.8), glm::vec3(0, 0.8, 1.1), glm::vec3(0, 1.2, 1), +// glm::vec3(0.33, 0.2, 0.12), glm::vec3(0.33, 0.5, 0.42), glm::vec3(0.33, 0.9, -0.62), +// glm::vec3(0.33, 1.1, -1.756), +// glm::vec3(0.66, 0.2, 0.12), glm::vec3(0.66, 0.5, 0.42), glm::vec3(0.66, 0.9, -0.62), +// glm::vec3(0.66, 1.0, -1.756), +// glm::vec3(1, 0.2, 0.8), glm::vec3(1, 0.5, 1), glm::vec3(1, 0.9, 1.1), glm::vec3(1, 1.2, 1) +// }}; +// f.weights = {4, 4, +// { +// 1, 1, 1, 1, +// 1, 1, 1, 1, +// 1, 1, 1, 1, +// 1, 1, 1, 1, +// } +// }; + + ifstream fin; + fin.open(R"(E:\qt\OpenGLDemo\surfacePlayer\assets\intersectTest\casea3\surfaces.txt)"); + + string str; + string tmp; + while (fin.good()) { + getline(fin, tmp); + str += tmp + "\n"; + } + cout << str << endl; + + auto infos = QString(str.c_str()).split(";"); + auto wData = infos[0]; + auto srf1Data = infos[1]; + auto srf2Data = infos[2]; + + auto points1 = getPtsFromStr(srf1Data); + auto points2 = getPtsFromStr(srf2Data); + + wData = wData.remove(0, wData.indexOf("Matrix") + 6).trimmed(); + wData.remove(0, 3); + wData.remove(wData.size() - 3, 3); + auto wInfos = wData.split(QRegularExpression("\\]( )*,( )*\\[")); + + array2 weights; + vector wArray(points1.cols() * points1.rows()); + for (int i = 0; i < points1.cols(); i++) { + auto wsInV = wInfos[i].split(QRegularExpression("( )*,( )*")); + for (int j = 0; j < points1.rows(); j++) { + wArray[i * points1.rows() + j] = wsInV[j].toFloat(); + } + } + weights = {points1.cols(), points2.rows(), wArray}; + + // knot + std::vector knot_u(2 * points1.cols(), 1.); + for (int i = 0; i < knot_u.size() / 2; i++)knot_u[i] = 0; + std::vector knot_v(2 * points1.rows(), 1.); + for (int i = 0; i < knot_v.size() / 2; i++)knot_v[i] = 0; + + s.degree_u = knot_u.size() / 2 - 1; + s.degree_v = knot_v.size() / 2 - 1; + s.knots_u = knot_u; + s.knots_v = knot_v; + s.control_points = points1; + s.weights = weights; + + f.degree_u = knot_u.size() / 2 - 1; + f.degree_v = knot_v.size() / 2 - 1; + f.knots_u = knot_u; + f.knots_v = knot_v; + f.control_points = points2; + f.weights = weights; + + fin.close(); + + // ====================== 测试 ======================= + vector > s_evaluation; + vector > f_evaluation; + // 曲面s和f的切向量。zd*-sznmj + vector > s_tangent_v; + vector > f_tangent_u; + const vector > f_tangent_v; // 曲面s和f的法向量 - const vector> s_normal; - const vector> f_normal; + const vector > s_normal; + const vector > f_normal; - auto mesh1 = getMesh(s, 7); + auto mesh1 = getMesh(s, 7); // [0, 63] auto mesh2 = getMesh(f, 7); SingularityJudger singularityJudger(s, f, mesh1, mesh2); - auto hasSingularity = singularityJudger.judge({2,7}, {4,12}, {3, 9}, {10, 16}); - for(auto line: singularityJudger.judgeRes) { - for(auto el: line) { - cout< getPtsFromStr(QString srfData) { + srfData = srfData.remove(0, srfData.indexOf("Matrix") + 6).trimmed(); + // 去掉首尾的括号 + srfData.remove(0, 1); + srfData.remove(srfData.size() - 1, 1); + auto srfInfos = srfData.split("{"); + auto srfDimInfos = srfInfos[0].split(QRegularExpression(",( )*")); + auto dimU = srfDimInfos[0].toInt(); + auto dimV = srfDimInfos[1].toInt(); + auto srfPtInfo = srfInfos[1].remove(srfInfos[1].indexOf("}"), 1); + auto srfPtsInfos = srfPtInfo.split(QRegularExpression(",( )*\\(")); + std::vector points(dimU * dimV, glm::vec3()); + for (int i = 0; i < srfPtsInfos.size(); i++) { + auto pt = srfPtsInfos[i].split("=")[1].trimmed(); + // 去掉首尾的方括号 + pt.remove(0, 1); + pt.remove(pt.size() - 1, 1); +// int iu = i / dimU; +// int iv = i % dimU; + auto coords = pt.split(QRegularExpression(",( )*")); + for (int k = 0; k < 3; k++) { + float num; + if (coords[k].contains("/")) { + auto nums = coords[k].split("/"); + num = nums[0].trimmed().toFloat() / nums[1].trimmed().toFloat(); + } else { + num = coords[k].toFloat(); + } + if (k == 0) points[i].x = num; + else if (k == 1) points[i].y = num; + else points[i].z = num; + } + } + array2 res = {(size_t) dimU, (size_t) dimV, points}; + return res; +} \ No newline at end of file diff --git a/src/C2C4.cpp b/src/C2C4.cpp index e27849b..0aba64c 100644 --- a/src/C2C4.cpp +++ b/src/C2C4.cpp @@ -31,30 +31,30 @@ pair C2C4::c2OrC4(int patchIdx_u1, int patchIdx_v1, int patchIdx_u2, if (i == 0 && j == 0) { der1_u = mesh1.tangent_u[patchIdx_u1][patchIdx_v1]; der1_v = mesh1.tangent_v[patchIdx_u1][patchIdx_v1]; - der2_u = mesh2.tangent_u[patchIdx_u2][patchIdx_v2]; - der2_v = mesh2.tangent_v[patchIdx_u2][patchIdx_v2]; + der2_u = -mesh2.tangent_u[patchIdx_u2][patchIdx_v2]; + der2_v = -mesh2.tangent_v[patchIdx_u2][patchIdx_v2]; } else if (i == 0 && j == reSampleCnt_v - 1) { der1_u = mesh1.tangent_u[patchIdx_u1][patchIdx_v1 + 1]; der1_v = mesh1.tangent_v[patchIdx_u1][patchIdx_v1 + 1]; - der2_u = mesh2.tangent_u[patchIdx_u2][patchIdx_v2 + 1]; - der2_v = mesh2.tangent_v[patchIdx_u2][patchIdx_v2 + 1]; + der2_u = -mesh2.tangent_u[patchIdx_u2][patchIdx_v2 + 1]; + der2_v = -mesh2.tangent_v[patchIdx_u2][patchIdx_v2 + 1]; } else if (i == reSampleCnt_u - 1 && j == 0) { der1_u = mesh1.tangent_u[patchIdx_u1 + 1][patchIdx_v1]; der1_v = mesh1.tangent_v[patchIdx_u1 + 1][patchIdx_v1]; - der2_u = mesh2.tangent_u[patchIdx_u2 + 1][patchIdx_v2]; - der2_v = mesh2.tangent_v[patchIdx_u2 + 1][patchIdx_v2]; + der2_u = -mesh2.tangent_u[patchIdx_u2 + 1][patchIdx_v2]; + der2_v = -mesh2.tangent_v[patchIdx_u2 + 1][patchIdx_v2]; } else if (i == reSampleCnt_u - 1 || j == reSampleCnt_v - 1) { der1_u = mesh1.tangent_u[patchIdx_u1 + 1][patchIdx_v1 + 1]; der1_v = mesh1.tangent_v[patchIdx_u1 + 1][patchIdx_v1 + 1]; - der2_u = mesh2.tangent_u[patchIdx_u2 + 1][patchIdx_v2 + 1]; - der2_v = mesh2.tangent_v[patchIdx_u2 + 1][patchIdx_v2 + 1]; + der2_u = -mesh2.tangent_u[patchIdx_u2 + 1][patchIdx_v2 + 1]; + der2_v = -mesh2.tangent_v[patchIdx_u2 + 1][patchIdx_v2 + 1]; } else { float v1 = v1Base + j * patchEdge_v1 / (reSampleCnt_v - 1); float v2 = v2Base + j * patchEdge_v2 / (reSampleCnt_v - 1); auto der1 = tinynurbs::surfaceDerivatives(srf1, 1, u1, v1); auto der2 = tinynurbs::surfaceDerivatives(srf2, 1, u2, v2); - der1_u = der1[2], der1_v = der1[1]; - der2_u = -der2[2], der2_v = -der2[1]; + der1_u = glm::normalize(der1(1, 0)), der1_v = glm::normalize(der1(0, 1)); + der2_u = -glm::normalize(der2(1, 0)), der2_v = -glm::normalize(der2(0, 1)); } jacobian[0][0].a = min(jacobian[0][0].a, der1_u.x), jacobian[0][0].b = max(jacobian[0][0].b, der1_u.x); jacobian[1][0].a = min(jacobian[1][0].a, der1_u.y), jacobian[1][0].b = max(jacobian[1][0].b, der1_u.y); diff --git a/src/SingularityJudger.cpp b/src/SingularityJudger.cpp index 0fcaf70..6c1e161 100644 --- a/src/SingularityJudger.cpp +++ b/src/SingularityJudger.cpp @@ -12,13 +12,15 @@ bool SingularityJudger::judge(pair focusRange_u1, pair focus // gauss map 只需要法向量信息,不需要evaluations GaussMap gaussMap1(mesh1.normal); GaussMap gaussMap2(mesh2.normal); + gaussMap1.build(); + gaussMap2.build(); if (!isGaussMapsOverlapped(gaussMap1, gaussMap2, focusRange_u1, focusRange_u2, focusRange_v1, focusRange_v2)) { // 一定没有交 return false; } // 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); + 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)); @@ -30,7 +32,7 @@ bool SingularityJudger::judge(pair focusRange_u1, pair focus 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曲面上的对应面片的坐标 @@ -39,6 +41,7 @@ bool SingularityJudger::judge(pair focusRange_u1, pair focus 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表示啥也没有 @@ -50,4 +53,21 @@ bool SingularityJudger::judge(pair focusRange_u1, pair focus SingularityJudger:: SingularityJudger(RationalSurface &srf1_, RationalSurface &srf2_, SrfMesh &mesh1_, SrfMesh &mesh2_) -:srf1(srf1_), srf2(srf2_), mesh1(mesh1_), mesh2(mesh2_){} + : srf1(srf1_), srf2(srf2_), mesh1(mesh1_), mesh2(mesh2_) {} + +void SingularityJudger::judge2(pair focusRange_u1, pair focusRange_v1, pair focusRange_u2, + pair focusRange_v2) { + GaussMap gaussMap1(mesh1.normal); + GaussMap gaussMap2(mesh2.normal); + gaussMap1.build(); + gaussMap2.build(); + printf("gauss map: %d\n", + isGaussMapsOverlapped(gaussMap1, gaussMap2, focusRange_u1, focusRange_u2, focusRange_v1, focusRange_v2)); + + +// gaussMap1.printQuadTree(); + + C2C4 c2C4(mesh1, mesh2, srf1, srf2); + printf("c2c4: %d\n", + c2C4.c2OrC4(focusRange_u1.first, focusRange_v1.first, focusRange_u2.first, focusRange_v2.first).first); +} diff --git a/src/gauss_map.cpp b/src/gauss_map.cpp index 0cfe199..17617e9 100644 --- a/src/gauss_map.cpp +++ b/src/gauss_map.cpp @@ -6,10 +6,8 @@ #include - -GaussMap::GaussMap(const std::vector> &normals_): normals(normals_) { - leafSampleCnt = int(normals_.size()); - maxLevel = int(log2(leafSampleCnt - 1) + 1); +GaussMap:: GaussMap(const std::vector> &normals_): normals(normals_) { + maxLevel = int(log2(int(normals_.size()) - 1) + 1); tree.resize(int((pow(4, maxLevel) - 1)) / 3); } @@ -18,8 +16,13 @@ void GaussMap::recursiveBuild(int level, int idx, int idx_u, int idx_v) { if (level == maxLevel) { // 叶子节点 for (int i = 0; i < 2; i++) - for (int j = 0; j < 2; j++) - bound = Union(bound, normals[idx_u + i][idx_v + j]); + for (int j = 0; j < 2; j++){ + auto normal = normals[idx_u + i][idx_v + j]; +// printf("normal: (%f, %f, %f)\n", normal.x, normal.y, normal.z); + if(glm::dot(normal, glm::vec3(1,1,1)) < 0)normal = -normal; + bound = Union(bound, normal); + } + tree[idx].nBound = bound; tree[idx].level = level; tree[idx].firstChild = -1; diff --git a/src/loop_detector.cpp b/src/loop_detector.cpp index d9a95aa..4ad5875 100644 --- a/src/loop_detector.cpp +++ b/src/loop_detector.cpp @@ -7,12 +7,12 @@ LoopDetector:: 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 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_normal(s_normal_), f_normal(f_normal_), s(s_) {} void LoopDetector::initOrientedDistance() { @@ -47,7 +47,8 @@ void LoopDetector::initOrientedDistance() { } } -void LoopDetector::initVectorField() { +// Yawei Ma的文章中的做法 +void LoopDetector::initVectorField0() { 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++) { @@ -74,9 +75,10 @@ void LoopDetector::initVectorField() { 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], pf_pq = f_tangent_v[f_subPatchIdxRange_u.first + + fPointIdx.first] [f_subPatchIdxRange_v.first + fPointIdx.second]; // 构造Aij矩阵,见 APPENDIX (A7) auto inverseAij = glm::inverse( @@ -96,31 +98,76 @@ void LoopDetector::initVectorField() { int b = 2; } +// G A Kriezis的文章中的做法 +void LoopDetector::initVectorField() { + 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.); + 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)); + } + } +} + void LoopDetector::getRotationNumber() { + // 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; + 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++) { for (int j = 0; j < s_subPatchEdgeSampleCnt_v - 1; j++) { - auto rotationNumber = 0.; - + vector> F(2, vector(2)); + vector> G(2, vector(2)); for (int biasI = 0; biasI < 2; biasI++) { for (int biasJ = 0; biasJ < 2; biasJ++) { - auto v = vectorFields[i + biasI][j + biasJ]; - rotationNumber += atan(double(v.x / v.y)); + auto idxI = i + biasI; + auto idxJ = j + biasJ; + auto v = vectorFields[idxI][idxJ]; + auto vSquareSum = v.x * v.x + v.y * v.y; + auto pTheta_pChi = -v.y / vSquareSum; + auto pTheta_pPhi = v.x / vSquareSum; +// auto pChi_pu = +// (vectorFields[idxI + 1][idxJ].x - vectorFields[idxI - 1][idxJ].x) / (2 * uParamCellSize); +// auto pChi_pv = +// (vectorFields[idxI][idxJ + 1].x - vectorFields[idxI][idxJ - 1].x) / (2 * vParamCellSize); + glm::vec2 pV_pu; + glm::vec2 pV_pv; + if (idxI == 0) + pV_pu = (vectorFields[idxI + 1][idxJ] - vectorFields[idxI][idxJ]) / uParamCellSize; + else if (idxI == s_subPatchEdgeSampleCnt_u - 1) + pV_pu = (vectorFields[idxI][idxJ] - vectorFields[idxI - 1][idxJ]) / uParamCellSize; + else + pV_pu = (vectorFields[idxI + 1][idxJ] - vectorFields[idxI - 1][idxJ]) / (2 * uParamCellSize); + if (idxJ == 0) + pV_pv = (vectorFields[idxI][idxJ + 1] - vectorFields[idxI][idxJ]) / vParamCellSize; + else if (idxJ == s_subPatchEdgeSampleCnt_v - 1) + pV_pv = (vectorFields[idxI][idxJ] - vectorFields[idxI][idxJ - 1]) / vParamCellSize; + else + pV_pv = (vectorFields[idxI][idxJ + 1] - vectorFields[idxI][idxJ - 1]) / (2 * vParamCellSize); + F[biasI][biasJ] = pTheta_pChi * pV_pu.x + pTheta_pPhi * pV_pu.y; + G[biasI][biasJ] = pTheta_pChi * pV_pv.x + pTheta_pPhi * pV_pv.y; } } - auto v00 = vectorFields[i][j]; - auto v01 = vectorFields[i][j + 1]; - auto v11 = vectorFields[i + 1][j + 1]; - auto v10 = vectorFields[i + 1][j]; - - rotationNumber = 2 * atan(double(v11.x / v11.y)) - 2 * atan(double(v00.x / v00.y)); - rotationNumber /= (2 * PI); - - rotationNumbers[i][j] = int(round(rotationNumber)); - - printf("rotationNumber: %lf\n", rotationNumber); + 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.)); } } + int a = 3; + int b = 1; } vector> LoopDetector::detect(pair _s_subPatchIdxRange_u, pair _s_subPatchIdxRange_v, @@ -143,3 +190,5 @@ vector> LoopDetector::detect(pair _s_subPatchIdxRange_u return {}; } + + diff --git a/test.txt b/test.txt new file mode 100644 index 0000000..af2f461 --- /dev/null +++ b/test.txt @@ -0,0 +1,3 @@ +123rrt +uyttt8i8 fg ,, +fds 546 lk \ No newline at end of file