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

subdivision

master
gjj 7 months ago
parent
f39ce0c03e
commit
f44e3efe47
3 changed files with 89 additions and 84 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. 123
      algoim/organizer/organizer.hpp
  2. 18
      algoim/organizer/primitive.hpp
  3. 32
      gjj/primitiveDebug.hpp

123
algoim/organizer/organizer.hpp
View File

@ -264,65 +264,55 @@ std::array<int, 2> sides = {-1, 1};
// 对于mark含2个及以上0的情况，尝试对每个0填1或-1的所有组合
// TODO: 参数太多了，考虑换用std::function + lambda
void visitSubcellOnBothSidesOfDir(const uvector3& nowNodeMin,
const uvector3& nowNodeMax,
const uvector3& nodeMid,
const int lastIdx,
void visitSubcellOnBothSidesOfDir(const uvector3& nodeMid,
const int polyIntersectIndex,
const tensor3& poly,
int startIdxToCheck,
uvector<int, 3>& mark,
std::vector<Node>& nodes)
std::array<Node, CHILD_NUM>& subNodes,
int startIdxToCheck,
uvector<int, 3>& mark)
{
int zeroIdx = findFirst(mark, 0, startIdxToCheck);
if (zeroIdx == -1) {
tensor3 halfCellPoly(nullptr, poly.ext());
algoim_spark_alloc(real, halfCellPoly);
auto subNodeMin = nowNodeMin;
auto subNodeMax = nowNodeMax;
for (int k = 0; k < 3; ++k) {
if (mark(k) == -1) {
subNodeMax(k) = nodeMid(k);
} else if (mark(k) == 1) {
subNodeMin(k) = nodeMid(k);
} else {
std::cerr << "error: mark = " << mark;
}
}
bernstein::deCasteljau(poly, subNodeMin, subNodeMax, halfCellPoly); // 求1/8空间下的表达
if (bernstein::uniformSign(halfCellPoly) != 1) {
nodes[lastIdx + binaryToDecimal(mark, -1)].polyIntersectIndices.emplace_back(polyIntersectIndex);
}
int subIdx = binaryToDecimal(mark);
auto& subNode = subNodes[subIdx];
bernstein::deCasteljau(poly, subNode.min, subNode.max, halfCellPoly); // 求1/8空间下的表达
if (bernstein::uniformSign(halfCellPoly) != 1) { subNode.polyIntersectIndices.emplace_back(polyIntersectIndex); }
} else {
for (auto side : sides) {
mark(zeroIdx) = side;
visitSubcellOnBothSidesOfDir(nowNodeMin,
nowNodeMax,
nodeMid,
lastIdx,
polyIntersectIndex,
poly,
zeroIdx + 1,
mark,
nodes);
visitSubcellOnBothSidesOfDir(nodeMid, polyIntersectIndex, poly, subNodes, zeroIdx + 1, mark);
}
mark(zeroIdx) = 0;
}
}
void build(const Scene& scene, std::vector<Node>& nodes, int nowNodeIdx)
// 叶节点单独存在leaves中，是因为最后只需要叶节点信息，没有自顶向下/自底向上遍历
// 所以树结构是不需要的记录的
// void build(const Scene& scene, std::vector<Node>& intermediateNodes, std::vector<Node> leaves, int nowNodeIdx)
void build(const Scene& scene, const Node& node, std::vector<Node> leaves)
{
const std::vector<int>& polyIntersectIndices = nodes[nowNodeIdx].polyIntersectIndices;
if (polyIntersectIndices.size() == 1) {
// TODO 相交很少时
const std::vector<int>& polyIntersectIndices = node.polyIntersectIndices;
if (polyIntersectIndices.size() <= 4) {
leaves.emplace_back(node);
return;
}
const uvector3& nowNodeMin = node.min;
const uvector3& nowNodeMax = node.max;
std::array<Node, CHILD_NUM> subNodes;
// intermediateNodes.resize(lastIdx + 8);
int subIdx = 0;
for (MultiLoop<3> j(0, 2); ~j; ++j, ++subIdx) {
auto nodeAABB = getOneEightCellAABB(node.min, node.max, j(), 0);
subNodes[subIdx].min = nodeAABB.first;
subNodes[subIdx].max = nodeAABB.second;
}
const uvector3& nowNodeMin = nodes[nowNodeIdx].min;
const uvector3& nowNodeMax = nodes[nowNodeIdx].max;
int lastIdx = nodes.size();
nodes.resize(lastIdx + 8);
uvector3 nodeMid = (nowNodeMin + nowNodeMax) / 2;
for (int i = 0; i < polyIntersectIndices.size(); ++i) {
const auto& poly = scene.polys[polyIntersectIndices[i]];
const int polyIntersectIndex = polyIntersectIndices[i];
const auto& poly = scene.polys[polyIntersectIndex];
uvector<int, 3> mark(0, 0, 0);
for (int faceAxis = 0; faceAxis < 3; ++faceAxis) {
real centerPlane = nodeMid(faceAxis);
@ -336,9 +326,8 @@ void build(const Scene& scene, std::vector<Node>& nodes, int nowNodeIdx)
} else if (signOnHalfPlane == 1) {
// primitive intersects either side or both sides of the centerface
// deCasteljau to transformation
// TODO: not 0, 1, or 0.5
uvector3 halfCellMin = 0, halfCellMax = 1;
halfCellMax(faceAxis) = 0.5;
uvector3 halfCellMin = nowNodeMin, halfCellMax = nowNodeMax;
halfCellMax(faceAxis) = nodeMid(faceAxis);
tensor3 halfCellPoly(nullptr, poly.ext());
algoim_spark_alloc(real, halfCellPoly);
bernstein::deCasteljau(poly, halfCellMin, halfCellMax, halfCellPoly);
@ -346,9 +335,8 @@ void build(const Scene& scene, std::vector<Node>& nodes, int nowNodeIdx)
// 负空间有
mark(faceAxis) = -1;
}
// TODO: not 0, 1, or 0.5
halfCellMax(faceAxis) = 1;
halfCellMin(faceAxis) = 0.5;
halfCellMax(faceAxis) = nowNodeMax(faceAxis);
halfCellMin(faceAxis) = nodeMid(faceAxis);
bernstein::deCasteljau(poly, halfCellMin, halfCellMax, halfCellPoly);
if (bernstein::uniformSign(halfCellPoly) != 1) {
// 正空间有
@ -358,55 +346,32 @@ void build(const Scene& scene, std::vector<Node>& nodes, int nowNodeIdx)
}
if (any(mark == uvector3(2, 2, 2))) {
int subIdx = 0;
for (MultiLoop<3> j(0, 2); ~j; ++j, ++subIdx) {
nodes[lastIdx + subIdx].polyIntersectIndices.emplace_back(poly);
nodes[lastIdx + subIdx].min = nowNodeMin;
nodes[lastIdx + subIdx].max = nowNodeMax;
for (int k = 0; k < 3; ++k) {
if (j(k) == 0)
nodes[lastIdx + subIdx].max(k) = nodeMid(k);
else
nodes[lastIdx + subIdx].min(k) = nodeMid(k);
}
for (int subIdx = 0; subIdx < CHILD_NUM; ++subIdx) {
// intermediateNodes[lastIdx + subIdx].polyIntersectIndices.emplace_back(poly);
// intermediateNodes[lastIdx + subIdx].min = nowNodeMin;
// intermediateNodes[lastIdx + subIdx].max = nowNodeMax;
subNodes[subIdx].polyIntersectIndices.emplace_back(polyIntersectIndex);
}
continue;
}
int zeroCount = numOfZero<3>(mark);
if (zeroCount == 0) {
// mark has -1 or 1 only
nodes[lastIdx + binaryToDecimal(mark, -1)].polyIntersectIndices.emplace_back(polyIntersectIndices[i]);
subNodes[binaryToDecimal(mark, -1)].polyIntersectIndices.emplace_back(polyIntersectIndex);
} else if (zeroCount == 1) {
// poly related to 2 subcells
uvector<int, 3> sideMark = mark;
int zeroIdx = replaceFirst(sideMark, 0, 1);
nodes[lastIdx + binaryToDecimal(sideMark, -1)].polyIntersectIndices.emplace_back(polyIntersectIndices[i]);
subNodes[binaryToDecimal(sideMark, -1)].polyIntersectIndices.emplace_back(polyIntersectIndex);
sideMark(zeroIdx) = -1;
nodes[lastIdx + binaryToDecimal(sideMark, -1)].polyIntersectIndices.emplace_back(polyIntersectIndices[i]);
subNodes[binaryToDecimal(sideMark, -1)].polyIntersectIndices.emplace_back(polyIntersectIndex);
} else {
// zeroCount == 2 or 3
visitSubcellOnBothSidesOfDir(nowNodeMin,
nowNodeMax,
nodeMid,
lastIdx,
polyIntersectIndices[i],
poly,
0,
mark,
nodes);
visitSubcellOnBothSidesOfDir(nodeMid, polyIntersectIndex, poly, subNodes, 0, mark);
}
}
// launch subdivision in subcells
int subIdx = 0;
for (MultiLoop<3> j(0, 2); ~j; ++j, ++subIdx) {
int idx = lastIdx + subIdx;
auto nodeAABB = getOneEightCellAABB(nodes[lastIdx + subIdx].min, nodes[lastIdx + subIdx].max, j(), 0);
nodes[idx].min = nodeAABB.first;
nodes[idx].max = nodeAABB.second;
nodes[nowNodeIdx].children[subIdx] = idx;
build(scene, nodes, idx);
}
for (subIdx = 0; subIdx < CHILD_NUM; ++subIdx) { build(scene, subNodes[subIdx], leaves); }
// TODO: 考虑fully contain 问题
}
}; // namespace algoim::Organizer

18
algoim/organizer/primitive.hpp
View File

@ -600,10 +600,20 @@ const int CHILD_NUM = 8;
struct Node {
std::vector<int> polyIntersectIndices; // 同一poly会出现在不同node的polyIndices中
std::vector<int> polyFullyContainedIndices; // 完全包含该node的poly, 同样会出现在不同node的polyIndices中
std::array<int, CHILD_NUM> children; // octree
uvector3 min, max;
Node() { children.fill(-1); }
// std::array<int, CHILD_NUM> children; // octree
uvector3 min, max;
// Node() { children.fill(-1); }
Node() = default;
// 深拷贝
Node(const Node& node)
: polyIntersectIndices(node.polyIntersectIndices),
polyFullyContainedIndices(node.polyFullyContainedIndices),
min(node.min),
max(node.max)
{
}
};
} // namespace algoim::Organizer

32
gjj/primitiveDebug.hpp
View File

@ -158,8 +158,38 @@ void casePolyhedronSphere()
basicTask(primitiveDescriptions);
}
void testDeCasteljau()
{
auto phiDesc = std::make_shared<SphereDesc>(SphereDesc(0.8, uvector3(0), 1.));
tensor3 tensor(nullptr, 3), tensor01(nullptr, 3), tensorBernstein(nullptr, 3), transformedTensorBernstein(nullptr, 3);
algoim_spark_alloc(real, tensor, tensor01, tensorBernstein, transformedTensorBernstein);
makeSphere(*phiDesc, tensor);
uvector3 xmax = 1, xmin = -1, range = xmax - xmin;
Organizer::detail::powerTransformation(range, xmin, tensor, tensor01);
Organizer::detail::power2BernsteinTensor(tensor01, tensorBernstein);
// bernstein::deCasteljau(tensorBernstein, uvector3(0), uvector3(0.2690598923241497 + 0.0000001),
// transformedTensorBernstein);
bernstein::deCasteljau(tensorBernstein,
uvector3(0.2690598923241497 + 0.0000001),
uvector3(1 - 0.2690598923241497 - 0.0000001),
transformedTensorBernstein);
uvector3 testX = uvector3(0.5);
int sign = bernstein::uniformSign(transformedTensorBernstein);
real evalX = evalPower(tensor, testX);
std::cout << "evalX original rep = " << evalX << std::endl;
evalX = bernstein::evalBernsteinPoly(tensorBernstein, testX);
std::cout << "evalX bernstein within 0-1 = " << evalX << std::endl;
evalX = bernstein::evalBernsteinPoly(transformedTensorBernstein, testX);
std::cout << "evalX bernstein after Decasteljau = " << evalX << std::endl;
std::cout << "sign = " << sign << std::endl;
}
void testPrimitive()
{
casePolyhedron1();
// casePolyhedron1();
// casePolyhedronSphere();
testDeCasteljau();
}
Write Preview
Loading…
Cancel
Save