Integration of gauss map, osculating toroidal patches, loop detection and C2 judgement to figure out the singular or loop intersection.
You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
 
 

392 lines
14 KiB

#include "QList"
#include "SingularityJudger.h"
#include "bvh.h"
#include "include/real.h"
#include "reader.hpp"
#include "utils.h"
#include <QRegularExpression>
#include <QString>
#include <array>
#include <iostream>
array2<glm::vec3> getPtsFromStr(QString srfData);
void printCtrPtsAsQuadruples(const RationalSurface<real> &s);
void sampleTimeTest(const RationalSurface<real> &s_, int sampleLevel) {
// 由于ParaSolid那边曲面参数为double类型,这里最好也是double类型(控制变量)
RationalSurface<real> s;
vector<real> knots_u(s_.knots_u.size());
vector<real> knots_v(s_.knots_v.size());
array2<real> weights(s_.weights.rows(), s_.weights.cols());
array2<glm::vec<3, real>> control_points(s_.control_points.rows(),
s_.control_points.cols());
for (int i = 0; i < s_.knots_u.size(); i++)
knots_u[i] = (real)s_.knots_u[i];
for (int i = 0; i < s_.knots_v.size(); i++)
knots_v[i] = (real)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) = (real)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 = (real)s_.control_points(i, j).x;
control_points(i, j).y = (real)s_.control_points(i, j).y;
control_points(i, j).z = (real)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) / real(sampleCnt - 1);
auto s_step_v = (*(s.knots_v.end() - 1) - s_first_v) / real(sampleCnt - 1);
// 为了分别测试赋值和求梯度的时间,这里把它们分开写了
auto startMomEval = std::chrono::steady_clock::now();
for (int i = 0; i < sampleCnt; i++) {
auto u = s_first_u + s_step_u * real(i);
for (int j = 0; j < sampleCnt; j++) {
auto v = s_first_v + s_step_v * real(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<real, std::milli>(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 * real(i);
for (int j = 0; j < sampleCnt; j++) {
auto v = s_first_v + s_step_v * real(j);
auto der = tinynurbs::surfaceDerivatives(s, 1, u, v);
// if(der(0, 0) == res.evaluation[i][j])
cout << "amazing" << endl;
// else cout<<"what??? ("<<res.evaluation[i][j].x<<"
// "<<res.evaluation[i][j].y<<"
// "<<res.evaluation[i][j].z<<") |
// ("<<der(0, 0).x<<" "<<der(0, 0).y<<" "<<der(0,
// 0).z<<")"<<endl; res.tangent_u[i][j] = der(1, 0);
// res.tangent_v[i][j] = der(0, 1);
}
}
auto endMomDer = std::chrono::steady_clock::now();
printf(
"time cost of derivatives: %lf ms\n",
std::chrono::duration<real, std::milli>(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 * real(i);
for (int j = 0; j < sampleCnt; j++) {
auto v = s_first_v + s_step_v * real(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<real, std::milli>(endMomScdDer - startMomScdDer)
.count());
}
void sampleTimeTestNonRational(const RationalSurface<real> &s_,
int sampleLevel) {
// 由于ParaSolid那边曲面参数为real类型,这里也要保证是real类型(控制变量)
Surface<real> s;
vector<real> knots_u(s_.knots_u.size());
vector<real> knots_v(s_.knots_v.size());
array2<glm::vec<3, real>> control_points(s_.control_points.rows(),
s_.control_points.cols());
for (int i = 0; i < s_.knots_u.size(); i++)
knots_u[i] = (real)s_.knots_u[i];
for (int i = 0; i < s_.knots_v.size(); i++)
knots_v[i] = (real)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 = (real)s_.control_points(i, j).x;
control_points(i, j).y = (real)s_.control_points(i, j).y;
control_points(i, j).z = (real)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) / real(sampleCnt - 1);
auto s_step_v = (*(s.knots_v.end() - 1) - s_first_v) / real(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 * real(i);
for (int j = 0; j < sampleCnt; j++) {
auto v = s_first_v + s_step_v * real(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<real, std::milli>(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 * real(i);
for (int j = 0; j < sampleCnt; j++) {
auto v = s_first_v + s_step_v * real(j);
auto der = tinynurbs::surfaceDerivatives(s, 1, u, v);
// if(der(0, 0) == res.evaluation[i][j])
cout << "amazing" << endl;
// else cout<<"what??? ("<<res.evaluation[i][j].x<<"
// "<<res.evaluation[i][j].y<<"
// "<<res.evaluation[i][j].z<<") |
// ("<<der(0, 0).x<<" "<<der(0, 0).y<<" "<<der(0,
// 0).z<<")"<<endl; res.tangent_u[i][j] = der(1, 0);
// res.tangent_v[i][j] = der(0, 1);
}
}
auto endMomDer = std::chrono::steady_clock::now();
printf(
"time cost of derivatives: %lf ms\n",
std::chrono::duration<real, std::milli>(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 * real(i);
for (int j = 0; j < sampleCnt; j++) {
auto v = s_first_v + s_step_v * real(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<real, std::milli>(endMomScdDer - startMomScdDer)
.count());
}
typedef array<int, 2> int2;
typedef array<int, 4> int4;
int dirs[4][2] = {{-1, 0}, {0, 1}, {1, 0}, {0, -1}};
// void dfs(
// const map<pair<int, int>, set<pair<int, int>>> &pairMap,
// // unordered_map<pair<int, int>, char> &book,
// set<pair<int, int>> &book,
// // vector<vector<pair<pair<int, int>, pair<int, int>>>>
// &boxGroups, vector<vector<pair<int, int>>> &boxGroups, int x, int y, int
// iOfGroup) {
// // book[{x, y}] = 1;
// // book.insert(pair<pair<int, int>, char>(pair<int, int>(x, y), 1));
// book.insert({x, y});
// boxGroups[iOfGroup].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() ||
// pairMap.find({nx, ny}) == pairMap.end())
// continue;
// dfs(pairMap, book, boxGroups, nx, ny, iOfGroup);
// }
// }
// book: 0: no pair, -1: has pair but not visited, positive int: has pair and
// grouped to the group {book-1} (the bookth group)
void dfsGroupingByOnlyOneSrf(vector<vector<int>> &book,
vector<vector<int2>> &boxGroups, int x, int y,
int iOfGroup) {
book[x][y] = iOfGroup + 1;
boxGroups[iOfGroup].emplace_back(int2{x, y});
int sizeX = book.size(), sizeY = book[0].size();
for (auto dir : dirs) {
auto nx = x + dir[0], ny = y + dir[1];
if (nx < 0 || nx >= sizeX || ny < 0 || ny >= sizeY)
continue;
if (book[nx][ny] != -1)
continue;
dfsGroupingByOnlyOneSrf(book, boxGroups, nx, ny, iOfGroup);
}
}
void testGroupingByOnlyOneSrf(
const vector<pair<pair<int, int>, pair<int, int>>> &pairs, int boxCntDir) {
utils::Timer timerGrouping(
"Grouping boxes filtered by BVH traversal using only the first surface.");
vector<vector<int>> book(boxCntDir, vector<int>(boxCntDir, 0));
vector<vector<int2>> boxGroups;
for (const auto &pair : pairs) {
const auto &boxIdx1 = pair.first;
book[boxIdx1.first][boxIdx1.second] = -1;
}
// map<int2, vector<int2>> box1ToBox2;
// for (const auto &pair : pairs) {
// const auto &boxIdx1 = pair.first;
// const auto &boxIdx2 = pair.second;
// box1ToBox2[{boxIdx1.first, boxIdx1.second}].emplace_back(
// int2{boxIdx2.first, boxIdx2.second});
// }
int groupNum = 0;
for (int i = 0; i < boxCntDir; i++) {
for (int j = 0; j < boxCntDir; j++) {
if (book[i][j] == -1) {
boxGroups.emplace_back();
dfsGroupingByOnlyOneSrf(book, boxGroups, i, j, groupNum++);
}
}
}
// from box groups to pair groups
vector<vector<int4>> pairGroups(groupNum);
for (const auto &pair : pairs) {
const auto &boxIdx1 = pair.first;
const auto &boxIdx2 = pair.second;
int groupIdx = book[boxIdx1.first][boxIdx1.second] - 1;
if (groupIdx < 0)
continue;
pairGroups[groupIdx].emplace_back(
int4{boxIdx1.first, boxIdx1.second, boxIdx2.first, boxIdx2.second});
}
timerGrouping.end();
printf("group num: %d\n", groupNum);
for (int i = 0; i < groupNum; i++) {
printf("number of patches of srf1 in group %d is: %lld", i,
boxGroups[i].size());
// for (const auto &box : boxGroups[i]) {
// printf("(%d, %d) ", box[0], box[1]);
// }
printf("\n");
}
for (int i = 0; i < groupNum; i++) {
printf("number of pairs in group %d is: %lld\n", i, pairGroups[i].size());
}
set<int2> aNums;
for (const auto &pair : pairs) {
const auto &firstBox = pair.first;
aNums.insert({firstBox.first, firstBox.second});
}
printf("Num of patch of surface 1: %lld\n", aNums.size());
}
int main() {
int level = 7;
printf("level: %d, sample cnt: %d * %d\n", level, int(pow(2, level - 1)),
int(pow(2, level - 1)));
// auto [s, f] = Reader::readSurfaces(R"(intersectTest\case17\surfaces.txt)");
auto s =
Reader::readSurface(R"(intersectTest\zyr_jh23_12_21\case2\surf_A.txt)");
auto f =
Reader::readSurface(R"(intersectTest\zyr_jh23_12_21\case2\surf_B.txt)");
// ====================== 测试 =======================
vector<vector<glm::vec3>> s_evaluation;
vector<vector<glm::vec3>> f_evaluation;
// 曲面s和f的切向量。zd*-sznmj
vector<vector<glm::vec3>> s_tangent_v;
vector<vector<glm::vec3>> f_tangent_u;
const vector<vector<glm::vec3>> f_tangent_v;
// 曲面s和f的法向量
const vector<vector<glm::vec3>> s_normal;
const vector<vector<glm::vec3>> f_normal;
// sampleTimeTestNonRational(s, level);
// sampleTimeTest(s, level);
auto mesh1 = SrfMesh(s, level);
auto mesh2 = SrfMesh(f, level);
BVH bvh1(mesh1.evaluation);
BVH bvh2(mesh2.evaluation);
bvh1.build();
bvh2.build();
vector<pair<pair<int, int>, pair<int, int>>> intersectBoxPairs =
getOverlapLeafNodes(bvh1, bvh2); // [{{u1, v1}, {u2, v2}}]
printf("box pairs size: %lld\n", intersectBoxPairs.size());
/**
* 测试对bvh结果用dfs
*/
// utils::Timer timerGrouping("Grouping boxes filtered by BVH traversal.");
// map<pair<int, int>, set<pair<int, int>>> pairMap;
// for (const auto &boxPair : intersectBoxPairs) {
// pairMap[boxPair.first].insert(boxPair.second);
// }
// // vector<vector<pair<pair<int, int>, pair<int, int>>>> groups;
// vector<vector<pair<int, int>>> groups;
// // unordered_map<pair<int, int>, char> book;
// set<pair<int, int>> 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++);
// }
// timerGrouping.end();
/**
* end of test
*/
testGroupingByOnlyOneSrf(intersectBoxPairs, int(pow(2, level - 1)));
utils::Timer timerSJ("Winding number detection");
SingularityJudger singularityJudger(s, f, mesh1, mesh2);
timerSJ.end();
pair<int, int> cellIdxFullRange = {0, pow(2, level - 1) - 1};
auto cellsWithCriticalPts = singularityJudger.judge(
intersectBoxPairs, cellIdxFullRange, cellIdxFullRange, cellIdxFullRange,
cellIdxFullRange);
// ==================
// 测试对整个曲面,gauss能排除多少(或保留多少)==================
// GaussMap gaussMap1(mesh1.normal);
// GaussMap gaussMap2(mesh2.normal);
// gaussMap1.build();
// gaussMap2.build();
// auto pairs = getOverlapLeafNodes(gaussMap1, gaussMap2);
// printf("Gauss Map: keep %lld samples in totally %lld boxes\n",
// pairs.size(),
// mesh1.normal.size() * mesh1.normal[0].size() *
// mesh2.normal.size() * mesh2.normal[0].size());
return 0;
}
void printCtrPtsAsQuadruples(const RationalSurface<real> &s) {
cout << endl;
for (int i = 0; i < s.control_points.rows(); i++) {
for (int j = 0; j < s.control_points.cols(); j++) {
auto pt = s.control_points(i, j);
auto w = s.weights(i, j);
cout << pt.x * w << ", " << pt.y * w << ", " << pt.z * w << ", "
<< s.weights(i, j) << ", " << endl;
}
}
cout << s.control_points.rows() * s.control_points.cols() * 4 << endl;
}