GJJ
/
SingularityJudger
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity
Browse Source

修复了一些bug，大幅提高了求解速度

master
Dtouch 2 years ago
parent
384d3ed18e
commit
be57891363
7 changed files with 270 additions and 115 deletions
Whitespace
Show all changes Ignore whitespace when comparing lines Ignore changes in amount of whitespace Ignore changes in whitespace at EOL
Split View
Diff Options
Show Stats Download Patch File Download Diff File
  1. 9
      include/SingularityJudger.h
  2. 14
      include/gauss_map.h
  3. 7
      include/loop_detector.h
  4. 66
      main.cpp
  5. 109
      src/SingularityJudger.cpp
  6. 54
      src/gauss_map.cpp
  7. 126
      src/loop_detector.cpp

9
include/SingularityJudger.h
View File

@ -6,6 +6,7 @@
#include "srf_mesh.h"
#include "gauss_map.h"
#include "tinynurbs/tinynurbs.h"
#include "loop_detector.h"
using namespace std;
using namespace tinynurbs;
@ -42,8 +43,16 @@ public:
vector<vector<char>> judgeRes;
private:
constexpr static int dirs[4][2] = {{-1, 0},
{0, 1},
{1, 0},
{0, -1}};
GaussMap gaussMap1;
GaussMap gaussMap2;
void
dfs(set<pair<int, int>> &book, vector<vector<pair<int, int>>> &cellGroups, const LoopDetector &loopDetector, int i,
int x, int y);
};

14
include/gauss_map.h
View File

@ -10,30 +10,34 @@
#include "glm/glm.hpp"
#include "aabb.h"
#include "map"
#include "set"
// 一个节点就表示一个球上矩形面片
struct GaussMapNode {
// 四个顶点的法向量值确定一个平面的法向量值
AABB bound;
int level, firstChild;
int idx_u, idx_v;
};
class GaussMap {
public:
int maxLevel;
const std::vector<std::vector<glm::vec3>> & normals;
const std::vector<std::vector<glm::vec3>> &normals;
std::vector<GaussMapNode> tree;
void recursiveBuild(int level, int idx, int idx_u, int idx_v);
explicit GaussMap(const std::vector<std::vector<glm::vec3>>& normals_);
explicit GaussMap(const std::vector<std::vector<glm::vec3>> &normals_);
void build();
void printQuadTree();
};
std::vector<std::pair<int, int>> getOverlapLeafNodes(const GaussMap &gm1, const GaussMap &gm2);
//std::vector<std::pair<int, int>> getOverlapLeafNodes(const GaussMap &gm1, const GaussMap &gm2);
std::map<std::pair<int, int>, std::set<std::pair<int, int>>>
getOverlapLeafNodes(const GaussMap &gm1, const GaussMap &gm2);
void recursiveGetOverlapLeafNodes(const GaussMap &gm1, const GaussMap &gm2, int idx1, int idx2,
std::vector<std::pair<int, int>> &pairs);
@ -50,6 +54,9 @@ bool isGaussMapsOverlapped(const GaussMap &gm1, const GaussMap &gm2, std::pair<i
std::pair<int, int> idxRange_v1, std::pair<int, int> idxRange_u2,
std::pair<int, int> idxRange_v2);
bool
isLeafNodesOverlapped(const GaussMap &gm1, const GaussMap &gm2, std::pair<int, int> idx1, std::pair<int, int> idx2);
/**
* 判断两个曲面的Gauss Map在指定的、各自的参数范围内有没有交
* @param range_u1 第一个gauss map的参数u范围
@ -74,6 +81,7 @@ bool isGaussMapsOverlapped(const GaussMap &gm1, const GaussMap &gm2, std::pair<f
* @param layer 层号。根节点为第一层
*/
int getStartIdxOfLayerN(int layer);
/**
* 根据u、v参数域中的小矩形的位置，判断它在BVH叶节点层（也就是最后一层）中的位置
*/

7
include/loop_detector.h
View File

@ -3,6 +3,8 @@
#include "tinynurbs/tinynurbs.h"
#include "glm/glm.hpp"
#include "set"
#include "map"
using namespace std;
@ -39,12 +41,13 @@ public:
// vector fields, 即有向距离关于u、v的导数
vector<vector<glm::vec2>> vectorFields;
vector<vector<int>> rotationNumbers;
vector<pair<int, int>> targetCells;
// int subPatchEdgeSampleCnt; // 在一个sub patch的采样网格中，边上的采样点个数
// void init(tinynurbs::RationalSurface<float> _s, tinynurbs::RationalSurface<float> _f, int _maxSplitLayer);
vector<pair<int, int>> detect(pair<int, int> _s_subPatchIdxRange_u, pair<int, int> _s_subPatchIdxRange_v,
pair<int, int> _f_subPatchIdxRange_u, pair<int, int> _f_subPatchIdxRange_v);
vector<pair<int, int>> detect(const vector<pair<pair<int, int>, pair<int, int>>> &intersectBoxPairs,
vector<pair<int, int>> detect(const map<pair<int, int>, set<pair<int, int>>> &intersectBoxPairs,
pair<int, int> _s_subPatchIdxRange_u, pair<int, int> _s_subPatchIdxRange_v,
pair<int, int> _f_subPatchIdxRange_u, pair<int, int> _f_subPatchIdxRange_v);
@ -68,6 +71,7 @@ private:
void initOrientedDistance();
void initOrientedDisAndVecFields(const vector<pair<pair<int, int>, pair<int, int>>> &intersectBoxPairs);
void initOrientedDisAndVecFields(const map<pair<int, int>, set<pair<int, int>>> &intersectBoxPairs);
// Yawei Ma的文章中的做法
void initVectorField0();
@ -77,6 +81,7 @@ private:
void getRotationNumber();
void getRotationNumber(const vector<pair<pair<int, int>, pair<int, int>>> &intersectBoxPairs);
void getRotationNumber(const map<pair<int, int>, set<pair<int, int>>> &intersectBoxPairs);
};

66
main.cpp
View File

@ -41,47 +41,9 @@ array2<glm::vec3> getPtsFromStr(QString srfData);
int main() {
RationalSurface<float> s;
RationalSurface<float> 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,
// }
// };
int level = 6;
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\case2\surfaces.txt)");
@ -150,29 +112,23 @@ int main() {
const vector<vector<glm::vec3>> s_normal;
const vector<vector<glm::vec3>> f_normal;
auto mesh1 = getMesh(s, 8); // [0, 63]
auto mesh2 = getMesh(f, 8);
auto mesh1 = getMesh(s, level);
auto mesh2 = getMesh(f, level);
BVH bvh1(mesh1.evaluation);
BVH bvh2(mesh2.evaluation);
bvh1.build();
bvh2.build();
auto intersectBoxPairs = getOverlapLeafNodes(bvh1, bvh2);
auto intersectBoxPairs = getOverlapLeafNodes(bvh1, bvh2); // [{{u1, v1}, {u2, v2}}]
printf("box pairs size: %lld\n", intersectBoxPairs.size());
double time_cost = utils::get_time();
SingularityJudger singularityJudger(s, f, mesh1, mesh2);
// for (auto pair: intersectBoxPairs) {
// singularityJudger.judge({bvh1.tree[pair.first].idx_u, bvh1.tree[pair.first].idx_u + 1},
// {bvh1.tree[pair.first].idx_v, bvh1.tree[pair.first].idx_v + 1},
// {bvh2.tree[pair.second].idx_u, bvh2.tree[pair.second].idx_u + 1},
// {bvh2.tree[pair.second].idx_v, bvh2.tree[pair.second].idx_v + 1});
// }
double time_cost = utils::get_time();
// singularityJudger.judge(intersectBoxPairs, {0, 255}, {0, 255}, {0, 255}, {0, 255});
singularityJudger.judge(intersectBoxPairs, {0, 127}, {0, 127}, {0, 127}, {0, 127});
// singularityJudger.judge(intersectBoxPairs, {0, 63}, {0, 63}, {0, 63}, {0, 63});
pair<int, int> cellIdxFullRange = {0, pow(2, level - 1) - 1};
singularityJudger.judge(intersectBoxPairs, cellIdxFullRange, cellIdxFullRange, cellIdxFullRange, cellIdxFullRange);
time_cost = utils::get_time() - time_cost;
printf("time cost: %lf\n", time_cost);
printf("time cost of singularityJudger: %lf\n", time_cost);
// ================== 测试对整个曲面，gauss能排除多少（或保留多少）==================
// GaussMap gaussMap1(mesh1.normal);

109
src/SingularityJudger.cpp
View File

@ -1,7 +1,6 @@
#include "SingularityJudger.h"
#include <utility>
#include "loop_detector.h"
#include "C2C4.h"
#include "set"
#include "utils.h"
@ -55,8 +54,8 @@ SingularityJudger(RationalSurface<float> &srf1_, RationalSurface<float> &srf2_,
: srf1(srf1_), srf2(srf2_), mesh1(mesh1_), mesh2(mesh2_),
gaussMap1(GaussMap(mesh1_.normal)), gaussMap2(GaussMap(mesh2_.normal)) {
gaussMapTimeCost = utils::get_time();
gaussMap1.build();
gaussMap2.build();
// gaussMap1.build();
// gaussMap2.build();
gaussMapTimeCost = utils::get_time() - gaussMapTimeCost;
}
@ -73,8 +72,8 @@ void SingularityJudger::judge2(pair<int, int> focusRange_u1, pair<int, int> focu
c2C4.c2OrC4(focusRange_u1.first, focusRange_v1.first, focusRange_u2.first, focusRange_v2.first).first);
}
bool hasPair(vector<pair<pair<int, int>, pair<int, int>>> pairs, pair<int, int> p) {
for (auto e: pairs) {
bool hasPair(const map<pair<int, int>, set<pair<int, int>>> &pairs, pair<int, int> p) {
for (const auto &e: pairs) {
if (p == e.first)return true;
}
return false;
@ -89,55 +88,83 @@ judge(const vector<pair<pair<int, int>, pair<int, int>>> &intersectBoxPairs,
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 (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})) {
if (isLeafNodesOverlapped(gaussMap1, gaussMap2, boxPair.first, boxPair.second)) {
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);
// if (!isGaussMapsOverlapped(gaussMap1, gaussMap2, focusRange_u1, focusRange_u2, focusRange_v1, focusRange_v2)) {
// // 一定没有交
// return;
// }
// 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);
loopDetector.detect(intersectBoxPairs, focusRange_u1, focusRange_v1, focusRange_u2, focusRange_v2);
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));
// 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) {
if (!hasPair(intersectBoxPairs, {focusRange_u1.first + i, focusRange_v1.first + j})) {
// 不在相交box pairs中
continue;
}
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表示啥也没有
}
for (auto targetCell: loopDetector.targetCells) {
printf("{%d, %d}, ", targetCell.first, targetCell.second);
}
printf("\ntarget cells cnt: %lld\n", loopDetector.targetCells.size());
/**
* 合并相邻盒子
*/
vector<vector<pair<int, int>>> cellGroups;
set<pair<int, int>> book;
int groupNum = 0;
for (auto targetCell: loopDetector.targetCells) {
if (book.find(targetCell) != book.end()) continue;
cellGroups.emplace_back();
dfs(book, cellGroups, loopDetector, groupNum++, targetCell.first, targetCell.second);
}
for(const auto& cellGroup: cellGroups)
{
printf("{");
for(int i = 0; i < cellGroup.size(); i++) {
printf("{%d, %d}", cellGroup[i].first, cellGroup[i].second);
if(i != cellGroup.size()-1)printf(", ");
}
printf("},\n");
}
// TODO c2 determination for tangency case
// C2C4 c2C4(mesh1, mesh2, srf1, srf2);
}
void SingularityJudger::dfs(set<pair<int, int>> &book, vector<vector<pair<int, int>>> &cellGroups,
const LoopDetector &loopDetector, int i, int x, int y) {
book.insert({x, y});
cellGroups[i].emplace_back(pair<int, int>(x, y));
for (auto dir: dirs) {
auto nx = x + dir[0], ny = y + dir[1];
if (book.find({nx, ny}) != book.end() || loopDetector.rotationNumbers[nx][ny] == 0)continue;
dfs(book, cellGroups, loopDetector, i, nx, ny);
}
}

54
src/gauss_map.cpp
View File

@ -6,7 +6,7 @@
#include <utility>
GaussMap:: GaussMap(const std::vector<std::vector<glm::vec3>> &normals_): normals(normals_) {
GaussMap::GaussMap(const std::vector<std::vector<glm::vec3>> &normals_) : normals(normals_) {
maxLevel = int(log2(int(normals_.size()) - 1) + 1);
tree.resize(int((pow(4, maxLevel) - 1)) / 3);
}
@ -14,13 +14,15 @@ GaussMap:: GaussMap(const std::vector<std::vector<glm::vec3>> &normals_): normal
void GaussMap::recursiveBuild(int level, int idx, int idx_u, int idx_v) {
AABB bound;
tree[idx].level = level;
tree[idx].idx_u = idx_u;
tree[idx].idx_v = idx_v;
if (level == maxLevel) {
// 叶子节点
for (int i = 0; i < 2; i++)
for (int j = 0; j < 2; 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;
if (glm::dot(normal, glm::vec3(1, 1, 1)) < 0)normal = -normal;
bound = Union(bound, normal);
}
@ -81,10 +83,21 @@ void GaussMap::printQuadTree() {
}
std::vector<std::pair<int, int>> getOverlapLeafNodes(const GaussMap &gm1, const GaussMap &gm2) {
//std::vector<std::pair<int, int>> getOverlapLeafNodes(const GaussMap &gm1, const GaussMap &gm2) {
// std::vector<std::pair<int, int>> resPairs;
// recursiveGetOverlapLeafNodes(gm1, gm2, 0, 0, resPairs);
// return resPairs;
//}
std::map<std::pair<int, int>, std::set<std::pair<int, int>>>
getOverlapLeafNodes(const GaussMap &gm1, const GaussMap &gm2) {
std::vector<std::pair<int, int>> resPairs;
recursiveGetOverlapLeafNodes(gm1, gm2, 0, 0, resPairs);
return resPairs;
std::map<std::pair<int, int>, std::set<std::pair<int, int>>> pairMap;
for (auto & resPair : resPairs) {
pairMap[{gm1.tree[resPair.first].idx_u, gm2.tree[resPair.first].idx_v}].
insert({gm1.tree[resPair.second].idx_u, gm2.tree[resPair.second].idx_v});
}
return pairMap;
}
void recursiveGetOverlapLeafNodes(const GaussMap &gm1, const GaussMap &gm2, int idx1, int idx2,
@ -164,8 +177,11 @@ bool isGaussMapsOverlapped(const GaussMap &gm1, const GaussMap &gm2, std::pair<i
AABB bounding;
for (int i = idxRange_u.first; i <= idxRange_u.second; ++i) {
for (int j = idxRange_v.first; j <= idxRange_v.second; ++j) {
auto tmp1 = getStartIdxOfLayerN(commonMaxLayer);
auto tmp2 = getChildNodeIdx(i, j);
auto globalIdx = getStartIdxOfLayerN(commonMaxLayer) + getChildNodeIdx(i, j);
bounding = Union(bounding,
gm.tree[getStartIdxOfLayerN(commonMaxLayer) + getChildNodeIdx(i, j)].bound);
gm.tree[globalIdx].bound);
}
}
return bounding;
@ -187,13 +203,33 @@ int getChildNodeIdx(int ix, int iy) {
int idx = 0;
while (ix > 0) {
int log2XCeil = int(log2f(ix));
idx += int(pow(4, log2XCeil));
idx += 2 * int(pow(4, log2XCeil));
ix -= int(pow(2, log2XCeil));
}
while (iy > 0) {
int log2XCeil = int(log2f(iy));
idx += 2 * int(pow(4, log2XCeil));
idx += int(pow(4, log2XCeil));
iy -= int(pow(2, log2XCeil));
}
return idx;
}
}
bool
isLeafNodesOverlapped(const GaussMap &gm1, const GaussMap &gm2, std::pair<int, int> idx1, std::pair<int, int> idx2) {
auto getLeafBox = [](const GaussMap &gm, const std::pair<int, int> &idx) {
AABB bounding;
if(idx.first==6 && idx.second == 7) {
int a = 1;
}
for (int i = 0; i < 2; i++) {
for (int j = 0; j < 2; j++) {
auto normal = gm.normals[idx.first + i][idx.second + j];
if (glm::dot(normal, glm::vec3(1, 1, 1)) < 0)normal = -normal;
bounding = Union(bounding, normal);
}
}
return bounding;
};
return IsOverlap(getLeafBox(gm1, idx1), getLeafBox(gm2, idx2));
}

126
src/loop_detector.cpp
View File

@ -103,6 +103,71 @@ void LoopDetector::initOrientedDisAndVecFields(const vector<pair<pair<int, int>,
}
}
void LoopDetector::initOrientedDisAndVecFields(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},
{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++) {
for (auto neighborIdx: neighborIdxes) {
float minDis = FLT_MAX;
int minDisFIdx_u = -1, minDisFIdx_v = -1;
int idxI = boxPair.first.first + i + neighborIdx[0];
int idxJ = boxPair.first.second + j + neighborIdx[1];
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 <= 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() ||
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};
// 计算vector fields
auto n2 = f_normal[minDisFIdx_u][minDisFIdx_v];
auto ps_pu = s_tangent_u[idxI][idxJ];
auto ps_pv = s_tangent_v[idxI][idxJ];
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的文章中的做法
void LoopDetector::initVectorField0() {
vectorFields = vector<vector<glm::vec2>>(s_subPatchEdgeSampleCnt_u, vector<glm::vec2>(s_subPatchEdgeSampleCnt_v));
@ -185,8 +250,7 @@ void LoopDetector::getRotationNumber() {
// 以小格子为单位遍历
for (int i = 0; i < s_subPatchEdgeSampleCnt_u - 1; i++) {
for (int j = 0; j < s_subPatchEdgeSampleCnt_v - 1; j++) {
vector<vector<float>> F(2, vector<float>(2));
vector<vector<float>> G(2, vector<float>(2));
float F[2][2], G[2][2];
for (int biasI = 0; biasI < 2; biasI++) {
for (int biasJ = 0; biasJ < 2; biasJ++) {
auto idxI = i + biasI;
@ -235,8 +299,7 @@ void LoopDetector::getRotationNumber(const vector<pair<pair<int, int>, pair<int,
vector<int>(s_subPatchEdgeSampleCnt_v));
// 以小格子为单位遍历
for (const auto &boxPair: intersectBoxPairs) {
vector<vector<float>> F(2, vector<float>(2));
vector<vector<float>> G(2, vector<float>(2));
float F[2][2], G[2][2];
for (int biasI = 0; biasI < 2; biasI++) {
for (int biasJ = 0; biasJ < 2; biasJ++) {
auto idxI = boxPair.first.first + biasI - s_subPatchIdxRange_u.first;
@ -274,6 +337,58 @@ void LoopDetector::getRotationNumber(const vector<pair<pair<int, int>, pair<int,
}
}
void LoopDetector::getRotationNumber(const map<pair<int, int>, set<pair<int, int>>> &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<vector<int>>(s_subPatchEdgeSampleCnt_u,
vector<int>(s_subPatchEdgeSampleCnt_v));
targetCells = {};
// 以小格子为单位遍历
for (const auto &boxPair: intersectBoxPairs) {
float F[2][2], G[2][2];
for (int biasI = 0; biasI < 2; biasI++) {
for (int biasJ = 0; biasJ < 2; biasJ++) {
auto idxI = boxPair.first.first + biasI - s_subPatchIdxRange_u.first;
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;
// 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;
}
}
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.);
rotationNumbers[boxPair.first.first - s_subPatchIdxRange_u.first][boxPair.first.second -
s_subPatchIdxRange_v.first] = rotationNumber;
if(rotationNumber != 0) {
targetCells.emplace_back(boxPair.first.first, boxPair.first.second);
}
}
}
vector<pair<int, int>> LoopDetector::detect(pair<int, int> _s_subPatchIdxRange_u, pair<int, int> _s_subPatchIdxRange_v,
pair<int, int> _f_subPatchIdxRange_u,
pair<int, int> _f_subPatchIdxRange_v) {
@ -294,7 +409,7 @@ vector<pair<int, int>> LoopDetector::detect(pair<int, int> _s_subPatchIdxRange_u
return {};
}
vector<pair<int, int>> LoopDetector::detect(const vector<pair<pair<int, int>, pair<int, int>>> &intersectBoxPairs,
vector<pair<int, int>> LoopDetector::detect(const map<pair<int, int>, set<pair<int, int>>> &intersectBoxPairs,
pair<int, int> _s_subPatchIdxRange_u, pair<int, int> _s_subPatchIdxRange_v,
pair<int, int> _f_subPatchIdxRange_u,
pair<int, int> _f_subPatchIdxRange_v) {
@ -310,7 +425,6 @@ vector<pair<int, int>> LoopDetector::detect(const vector<pair<pair<int, int>, pa
f_subPatchEdgeSampleCnt_v = f_subPatchIdxRange_v.second - f_subPatchIdxRange_v.first + 2;
initOrientedDisAndVecFields(intersectBoxPairs);
int a = 1;
getRotationNumber(intersectBoxPairs);
return {};
}

Write Preview
Loading…
Cancel
Save