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

origin

master
Dtouch 3 years ago
commit
a1efe768ea
6 changed files with 359 additions and 0 deletions
Whitespace
Show all changes Ignore whitespace when comparing lines Ignore changes in amount of whitespace Ignore changes in whitespace at EOL
Unified View
Diff Options
Show Stats Download Patch File Download Diff File
  1. 36
      .gitignore
  2. 9
      CMakeLists.txt
  3. 3
      README.md
  4. 64
      include/loop_detector.h
  5. 77
      main.cpp
  6. 170
      src/loop_detector.cpp

36
.gitignore
View File

@ -0,0 +1,36 @@
# ---> C++
# Prerequisites
*.d
# Compiled Object files
*.slo
*.lo
*.o
*.obj
# Precompiled Headers
*.gch
*.pch
# Compiled Dynamic libraries
*.so
*.dylib
*.dll
# Fortran module files
*.mod
*.smod
# Compiled Static libraries
*.lai
*.la
*.a
*.lib
# Executables
*.exe
*.out
*.app
cmake-build-debug/
.idea/

9
CMakeLists.txt
View File

@ -0,0 +1,9 @@
cmake_minimum_required(VERSION 3.21)
project(LoopDetector)
set(CMAKE_CXX_STANDARD 14)
add_executable(LoopDetector main.cpp src/loop_detector.cpp include/loop_detector.h)
include_directories(E:/CLib/tinynurbs/include E:/CLib/glm)
include_directories(include)

3
README.md
View File

@ -0,0 +1,3 @@
# GaussMap
A Gauss mapping tool for Nurbs surfaces on CPUs.

64
include/loop_detector.h
View File

@ -0,0 +1,64 @@
#ifndef LOOPDETECTOR_LOOP_DETECTOR_H
#define LOOPDETECTOR_LOOP_DETECTOR_H
#include "tinynurbs/tinynurbs.h"
#include "glm/glm.hpp"
using namespace std;
class LoopDetector {
public:
// 两个曲面
tinynurbs::RationalSurface<float> s;
tinynurbs::RationalSurface<float> f;
// 细分采样的总层数（最高与BVH的总层数相等）。从第二层开始，每一层都表示把原来的一个小sub patch分成了更小的四份
int maxSplitLayer;
// 采样点的求值和求切向量先设为public，因为这些因该是算好的，LoopDetector不用重复计算
// 曲面s和f的采样点。
vector<vector<glm::vec3>> s_evaluation;
vector<vector<glm::vec3>> f_evaluation;
// 曲面s和f的切向量。
vector<vector<glm::vec3>> s_tangent_u;
vector<vector<glm::vec3>> s_tangent_v;
vector<vector<glm::vec3>> f_tangent_u;
vector<vector<glm::vec3>> f_tangent_v;
// 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);
private:
// 需要做Loop检测的sub patch在最后一层上的下标的范围(每个范围都真包含于[0, 2^(maxSplitLayer-1)-1])
pair<int, int> s_subPatchIdxRange_u;
pair<int, int> s_subPatchIdxRange_v;
pair<int, int> f_subPatchIdxRange_u;
pair<int, int> f_subPatchIdxRange_v;
// subPatch每条边上的采样点个数
// 这四个值根据subPatchRange的范围得到，不可被外部修改
int s_subPatchEdgeSampleCnt_u;
int s_subPatchEdgeSampleCnt_v;
int f_subPatchEdgeSampleCnt_u;
int f_subPatchEdgeSampleCnt_v;
// 整个曲面一条边上的采样点个数
int edgeSampleCnt;
// 有向距离计算结果。有向距离通过采样的方式计算，其结果与s曲面sub patch中的采样网格规模相同，即与s_evaluation大小一致
// vector<vector<glm::vec3>> distance;
// 确定有向距离后，对应的f曲面上的最短距离点的位置
vector<vector<pair<int, int>>> selectedPointsIdx;
// vector fields, 即有向距离关于u、v的导数
vector<vector<glm::vec2>> vectorFields;
void initEvaluation();
void initOrientedDistance();
void initVectorField();
void getRotationNumber();
};
#endif //LOOPDETECTOR_LOOP_DETECTOR_H

77
main.cpp
View File

@ -0,0 +1,77 @@
#include <iostream>
#include "loop_detector.h"
int main() {
LoopDetector loopDetector;
tinynurbs::RationalSurface<float> s;
tinynurbs::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,
}
};
loopDetector.s = s;
loopDetector.f = f;
loopDetector.maxSplitLayer = 6;
// 需要做Loop检测的sub patch在最后一层上的下标的范围(每个范围都真包含于[0, 2^(maxSplitLayer-1)-1])
// 这里范围真包含于[0, 31]
loopDetector.detect({3, 11}, {4, 11}, {2, 7}, {6, 15});
glm::vec3 a(2, 3, 4);
glm::vec3 b(3, 7, 1);
auto nab = glm::normalize(a - b);
auto res = nab * (a - b);
auto ab = a * b;
cout << res.x << ", " << res.y << ", " << res.z << endl;
cout << nab.x << ", " << nab.y << ", " << nab.z << endl;
cout << ab.x << ", " << ab.y << ", " << ab.z << endl;
std::cout << "Hello, World!" << std::endl;
auto m = glm::mat2x2(1, 2, 3, 4);
cout << m[0][0] << " " << m[0][1] << endl;
auto dotRes = glm::dot(a, b);
cout << dotRes << endl;
cout<<isnan(glm::normalize(glm::vec3(0, 0, 0)).x)<<endl;
cout<<isnan(1.)<<endl;
// auto m1 = glm::mat2(1., 2.);
cout<<round(3.7)<<" "<<round(3.5657567658)<<" "<<round(3.2)<<" "<<round(3.0)<<endl;
return 0;
}

170
src/loop_detector.cpp
View File

@ -0,0 +1,170 @@
#include "loop_detector.h"
#include <utility>
const float PI = 3.1415927;
void LoopDetector::initEvaluation() {
// 初始化subPatch上所有采样点的值，以及它们的切向量，并存储
// 当该工具集成到系统中后，这一部分应该是可以直接全部去掉的
// 因为在GPU上已经做过对每个点的最细粒度的evaluation了，只需要从整个曲面把需要的subPatch找出来即可
// subPatch每条边上的采样点个数。边上的格子个数=range.second-range.first+1，采样点个数=格子个数+1
s_subPatchEdgeSampleCnt_u = s_subPatchIdxRange_u.second - s_subPatchIdxRange_u.first + 2;
s_subPatchEdgeSampleCnt_v = s_subPatchIdxRange_v.second - s_subPatchIdxRange_v.first + 2;
f_subPatchEdgeSampleCnt_u = f_subPatchIdxRange_u.second - f_subPatchIdxRange_u.first + 2;
f_subPatchEdgeSampleCnt_v = f_subPatchIdxRange_v.second - f_subPatchIdxRange_v.first + 2;
edgeSampleCnt = int(pow(2, maxSplitLayer - 1)) + 1;
s_evaluation = vector<vector<glm::vec3>>(s_subPatchEdgeSampleCnt_u, vector<glm::vec3>(s_subPatchEdgeSampleCnt_v));
f_evaluation = vector<vector<glm::vec3>>(f_subPatchEdgeSampleCnt_u, vector<glm::vec3>(f_subPatchEdgeSampleCnt_v));
// 曲面s和f的切向量。
s_tangent_u = vector<vector<glm::vec3>>(s_subPatchEdgeSampleCnt_u, vector<glm::vec3>(s_subPatchEdgeSampleCnt_v));
s_tangent_v = vector<vector<glm::vec3>>(s_subPatchEdgeSampleCnt_u, vector<glm::vec3>(s_subPatchEdgeSampleCnt_v));
f_tangent_u = vector<vector<glm::vec3>>(f_subPatchEdgeSampleCnt_u, vector<glm::vec3>(f_subPatchEdgeSampleCnt_v));
f_tangent_v = vector<vector<glm::vec3>>(f_subPatchEdgeSampleCnt_u, vector<glm::vec3>(f_subPatchEdgeSampleCnt_v));
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) / float(edgeSampleCnt - 1);
auto s_step_v = (*(s.knots_v.end() - 1) - s_first_v) / float(edgeSampleCnt - 1);
auto f_first_u = *(f.knots_u.begin());
auto f_first_v = *(f.knots_v.begin());
auto f_step_u = (*(f.knots_u.end() - 1) - f_first_u) / float(edgeSampleCnt - 1);
auto f_step_v = (*(f.knots_v.end() - 1) - f_first_v) / float(edgeSampleCnt - 1);
for (int i = 0; i < s_subPatchEdgeSampleCnt_u; i++) {
auto u = s_first_u + s_step_u * float(s_subPatchIdxRange_u.first + i);
for (int j = 0; j < s_subPatchEdgeSampleCnt_v; j++) {
auto v = s_first_v + s_step_v * float(s_subPatchIdxRange_v.first + j);
s_evaluation[i][j] = tinynurbs::surfacePoint(s, u, v);
auto der = tinynurbs::surfaceDerivatives(s, 1, u, v);
s_tangent_u[i][j] = der(1, 0);
s_tangent_v[i][j] = der(0, 1);
}
}
for (int i = 0; i < f_subPatchEdgeSampleCnt_u; i++) {
auto u = f_first_u + f_step_u * float(f_subPatchIdxRange_u.first + i);
for (int j = 0; j < f_subPatchEdgeSampleCnt_v; j++) {
auto v = f_first_v + f_step_v * float(f_subPatchIdxRange_v.first + j);
f_evaluation[i][j] = tinynurbs::surfacePoint(f, u, v);
auto der = tinynurbs::surfaceDerivatives(f, 1, u, v);
f_tangent_u[i][j] = der(1, 0);
f_tangent_v[i][j] = der(0, 1);
}
}
}
void LoopDetector::initOrientedDistance() {
selectedPointsIdx = vector<vector<pair<int, int>>>(s_subPatchEdgeSampleCnt_u,
vector<pair<int, int>>(s_subPatchEdgeSampleCnt_v));
// distance = vector<vector<glm::vec3>>(s_subPatchEdgeSampleCnt_u, vector<glm::vec3>(s_subPatchEdgeSampleCnt_v));
for (int i = 0; i < s_subPatchEdgeSampleCnt_u; i++) {
for (int j = 0; j < s_subPatchEdgeSampleCnt_v; j++) {
float minDis = FLT_MAX;
int minDisFIdx_u = -1, minDisFIdx_v = -1;
for (int k = 0; k < f_subPatchEdgeSampleCnt_u; k++) {
for (int l = 0; l < f_subPatchEdgeSampleCnt_v; l++) {
auto dis = glm::distance(s_evaluation[i][j], f_evaluation[k][l]);
// 确定f上的对应点时，依据的最小欧氏距离，而不是有向距离
if (dis < minDis) {
minDis = dis;
minDisFIdx_u = k;
minDisFIdx_v = l;
}
}
}
// 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};
}
}
}
void LoopDetector::initVectorField() {
vectorFields = vector<vector<glm::vec2>>(s_subPatchEdgeSampleCnt_u, vector<glm::vec2>(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];
auto N = glm::normalize(s_evaluation[i][j] - f_evaluation[fPointIdx.first][fPointIdx.second]);
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]);
N = glm::normalize(glm::cross(s_tangent_u[i][j], s_tangent_v[i][j]) +
glm::cross(f_tangent_u[fPointIdx.first][fPointIdx.second],
f_tangent_v[fPointIdx.first][fPointIdx.second]));
}
// u1、u2两个向量构成和N垂直的曲面
glm::vec3 u1(1, 1, 1);
if (N.z != 0)u1.z = (-N.x - N.y) / N.z;
else if (N.y != 0)u1.y = (-N.x - N.z) / N.y;
else u1.x = (-N.y - N.z) / N.x;
u1 = glm::normalize(u1);
auto u2 = glm::normalize(glm::cross(N, u1));
// s，f曲面在两个方向上的偏导
auto ps_pu = s_tangent_u[i][j], ps_pv = s_tangent_v[i][j];
auto pf_pp = f_tangent_u[fPointIdx.first][fPointIdx.second], pf_pq = f_tangent_v[fPointIdx.first][fPointIdx.second];
// 构造Aij矩阵，见<Detection of loops and singularities of surface intersections> 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))
);
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));
auto vNfpNfq_inverseAij = glm::vec2(glm::dot(N, pf_pp), glm::dot(N, pf_pq)) * inverseAij;
auto pd_pu = glm::dot(N, ps_pu) -
glm::dot(vNfpNfq_inverseAij, glm::vec2(glm::dot(u1, ps_pu), glm::dot(u2, ps_pu)));
auto pd_pv = glm::dot(N, ps_pv) -
glm::dot(vNfpNfq_inverseAij, glm::vec2(glm::dot(u1, ps_pv), glm::dot(u2, ps_pv)));
vectorFields[i][j] = glm::vec2(pd_pu, pd_pv);
}
}
int a = 1;
int b = 2;
}
void LoopDetector::getRotationNumber() {
// 以小格子为单位遍历
for (int i = 0; i < s_subPatchEdgeSampleCnt_u - 1; i++) {
for (int j = 0; j < s_subPatchEdgeSampleCnt_v - 1; j++) {
auto rotationNumber = 0.;
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 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);
printf("rotationNumber: %lf\n", rotationNumber);
}
}
}
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) {
s_subPatchIdxRange_u = _s_subPatchIdxRange_u;
s_subPatchIdxRange_v = _s_subPatchIdxRange_v;
f_subPatchIdxRange_u = _f_subPatchIdxRange_u;
f_subPatchIdxRange_v = _f_subPatchIdxRange_v;
initEvaluation(); // 当evaluation已经完成，这一步就不用了，直接传入两个evaluation矩阵和四个tangent矩阵的值即可
initOrientedDistance();
initVectorField();
getRotationNumber();
return {};
}
Write Preview
Loading…
Cancel
Save