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.
 
 

605 lines
20 KiB

#pragma once
#include <cassert>
#include <cmath>
#include <cstddef>
#include <vector>
#include "iostream"
#include "multiloop.hpp"
#include "polyset.hpp"
#include "sparkstack.hpp"
#include "uvector.hpp"
#include "real.hpp"
#include "xarray.hpp"
#include "binomial.hpp"
#include "bernstein.hpp"
#include <memory>
namespace algoim::Organizer
{
namespace detail
{
template <int N>
void power2BernsteinTensor(const xarray<real, N>& phiPower, xarray<real, N>& phiBernsetin)
{
xarrayInit(phiBernsetin);
for (auto i = phiPower.loop(); ~i; ++i) {
// phi.l(i) = powerFactors.l(i);
real factorBase = phiPower.l(i);
if (factorBase == 0) continue;
auto traverseRange = phiPower.ext() - i();
std::vector<std::vector<real>> decompFactors(N, std::vector<real>(max(traverseRange), 0.));
for (int dim = 0; dim < N; ++dim) {
// Sigma
size_t nDim = phiPower.ext()(dim) - 1;
const real* binomNDim = Binomial::row(nDim);
for (int j = i(dim); j <= nDim; ++j) {
const real* binomJ = Binomial::row(j);
decompFactors[dim][j - i(dim)] = binomJ[i(dim)] / binomNDim[i(dim)];
}
}
xarray<real, N> subgrid(nullptr, traverseRange);
// algoim_spark_alloc(real, subgrid);
for (auto ii = subgrid.loop(); ~ii; ++ii) {
real factor = factorBase;
for (int dim = 0; dim < N; ++dim) { factor *= decompFactors[dim][ii(dim)]; }
phiBernsetin.m(i() + ii()) += factor;
}
}
}
void powerTransformation(const uvector<real, 3>& scale,
const uvector<real, 3>& bias,
const xarray<real, 3>& origin,
xarray<real, 3>& res)
{
assert(all(origin.ext() == res.ext()));
std::vector<std::vector<std::vector<real>>> dimOrderExpansion;
const auto& ext = origin.ext();
for (int dim = 0; dim < 3; ++dim) {
dimOrderExpansion.push_back(std::vector<std::vector<real>>(ext(dim)));
for (int degree = 0; degree < ext(dim); ++degree) {
const real* binomDegree = Binomial::row(degree);
dimOrderExpansion[dim][degree].reserve(degree + 1);
// 根据二项定理展开
for (int i = 0; i <= degree; ++i) {
dimOrderExpansion[dim][degree].push_back(binomDegree[i] * pow(scale(dim), i) * pow(bias(dim), degree - i));
}
}
}
for (auto i = origin.loop(); ~i; ++i) {
// 迭代器必须按照坐标的升序进行访问,即,访问ijk时,(i-1)jk,i(j-1)k,ij(k-1)必须已经访问过
real base = origin.l(i);
res.l(i) = 0;
auto exps = i();
for (MultiLoop<3> j(0, exps + 1); ~j; ++j) {
real item = base;
for (int dim = 0; dim < 3; ++dim) { item *= dimOrderExpansion[dim][exps(dim)][j(dim)]; }
res.m(j()) += item;
}
}
}
void powerTransformation(const uvector<real, 3>& scale, const uvector<real, 3>& bias, xarray<real, 3>& phiPower)
{
std::vector<std::vector<std::vector<real>>> dimOrderExpansion;
const auto& ext = phiPower.ext();
for (int dim = 0; dim < 3; ++dim) {
dimOrderExpansion.push_back(std::vector<std::vector<real>>(ext(dim)));
for (int degree = 0; degree < ext(dim); ++degree) {
const real* binomDegree = Binomial::row(degree);
dimOrderExpansion[dim][degree].reserve(degree + 1);
// 根据二项定理展开
for (int i = 0; i <= degree; ++i) {
dimOrderExpansion[dim][degree].push_back(binomDegree[i] * pow(scale(dim), i) * pow(bias(dim), degree - i));
}
}
}
for (auto i = phiPower.loop(); ~i; ++i) {
// 迭代器必须按照坐标的升序进行访问,即,访问ijk时,(i-1)jk,i(j-1)k,ij(k-1)必须已经访问过
real base = phiPower.l(i);
phiPower.l(i) = 0;
auto exps = i();
for (MultiLoop<3> j(0, exps + 1); ~j; ++j) {
real item = base;
for (int dim = 0; dim < 3; ++dim) { item *= dimOrderExpansion[dim][exps(dim)][j(dim)]; }
phiPower.m(j()) += item;
}
}
}
static void compositePower(const std::vector<xarray<real, 3>>& powers,
int powerIdx,
uvector<int, 3> powerSum,
real factor,
xarray<real, 3>& res)
{
if (powerIdx == 0) {
{
uvector3 ext(1, 1, 1);
for (auto& t : powers) { ext += t.ext() - 1; }
assert(all(ext == res.ext()));
}
xarrayInit(res);
}
if (powerIdx == powers.size()) {
res.m(powerSum) += factor;
return;
}
auto& power = powers[powerIdx];
for (auto i = power.loop(); ~i; ++i) {
if (power.l(i) == 0) {
// factor = 0;
continue;
}
compositePower(powers, powerIdx + 1, powerSum + i(), factor * power.l(i), res);
}
}
} // namespace detail
class Primitive
{
public:
virtual void print() = 0;
virtual real eval(const uvector3&) { return 0; }
};
template <int N>
real evalPower(const xarray<real, N>& phi, const uvector<real, N>& x)
{
real res = 0;
for (auto i = phi.loop(); ~i; ++i) {
real item = phi.l(i);
auto exps = i();
for (int dim = 0; dim < N; ++dim) { item *= pow(x(dim), exps(dim)); }
res += item;
}
return res;
}
template <int N>
real evalBernstein(const xarray<real, N>& phi, const uvector<real, N>& x)
{
return bernstein::evalBernsteinPoly(phi, x);
}
// class PowerTensor : public Primitive
// {
// public:
// xarray<real, 3> tensor;
// // SparkStack<real>* sparkStackPtr;
// void print() override { std::cout << "Power" << std::endl; }
// real eval(const uvector3& p) override { return evalPower(tensor, p); }
// // PowerTensor() {}
// PowerTensor(const xarray<real, 3>& t_) : tensor(t_) {}
// // const xarray<real, 3>& getTensor() { return tensor; }
// ~PowerTensor() = default;
// };
// class PowerTensorComplex : public Primitive
// {
// public:
// xarray<real, 3> compositeTensor; // 复合后的张量
// SparkStack<real>* sparkStackPtr;
// std::vector<xarray<real, 3>> tensors; // 原始张量
// // std::vector<PowerTensor> powerTensors; // 原始张量
// void print() override { std::cout << "PowerTensorComplex" << std::endl; }
// static void compositePower(const std::vector<PowerTensor>& powers,
// int powerIdx,
// uvector<int, 3> powerSum,
// real factor,
// xarray<real, 3>& res)
// {
// // const xarray<real, 3>& tensor
// if (powerIdx == 0) {
// {
// uvector3 ext(1, 1, 1);
// for (auto& t : powers) { ext += t.tensor.ext() - 1; }
// assert(all(ext == res.ext()));
// }
// xarrayInit(res);
// }
// if (powerIdx == powers.size()) {
// res.m(powerSum) += factor;
// return;
// }
// auto& tensor = powers[powerIdx].tensor;
// for (auto i = tensor.loop(); ~i; ++i) {
// if (tensor.l(i) == 0) {
// //factor = 0;
// continue;
// }
// compositePower(powers, powerIdx + 1, powerSum + i(), factor * tensor.l(i), res);
// }
// }
// static void compositePower(const std::vector<xarray<real, 3>>& powers,
// int powerIdx,
// uvector<int, 3> powerSum,
// real factor,
// xarray<real, 3>& res)
// {
// if (powerIdx == 0) {
// {
// uvector3 ext(1, 1, 1);
// for (auto& t : powers) { ext += t.ext() - 1; }
// assert(all(ext == res.ext()));
// }
// xarrayInit(res);
// }
// if (powerIdx == powers.size()) {
// res.m(powerSum) += factor;
// return;
// }
// auto& power = powers[powerIdx];
// for (auto i = power.loop(); ~i; ++i) {
// if (power.l(i) == 0) {
// // factor = 0;
// continue;
// }
// compositePower(powers, powerIdx + 1, powerSum + i(), factor * power.l(i), res);
// }
// }
// // PowerTensorComplex() {}
// PowerTensorComplex(const std::vector<xarray<real, 3>>& ts_, xarray<real, 3>& ct_) : tensors(ts_), compositeTensor(ct_)
// {
// uvector3 ext(1);
// for (auto& t : ts_) { ext += t.ext() - 1; }
// assert(all(ext == compositeTensor.ext()));
// compositePower(tensors, 0, uvector3(0), 1, compositeTensor);
// }
// PowerTensorComplex(const std::vector<PowerTensor>& pts_, xarray<real, 3>& ct_) : compositeTensor(ct_)
// {
// uvector3 ext(1);
// tensors.resize(pts_.size());
// for (int i = 0; i < pts_.size(); ++i) {
// tensors[i] = pts_[i].tensor;
// ext += pts_[i].tensor.ext() - 1;
// }
// assert(all(ext == compositeTensor.ext()));
// compositePower(tensors, 0, uvector3(0), 1, compositeTensor);
// }
// PowerTensorComplex(const std::vector<PowerTensorComplex>& ptcs_, xarray<real, 3>& ct_) : compositeTensor(ct_)
// {
// std::vector<xarray<real, 3>> originCompositeTensors;
// uvector3 ext(1);
// for (auto& ptc : ptcs_) {
// for (auto& t : ptc.tensors) { tensors.emplace_back(t); }
// originCompositeTensors.push_back(ptc.compositeTensor);
// ext += ptc.compositeTensor.ext() - 1;
// }
// compositePower(originCompositeTensors, 0, uvector3(0, 0, 0), 1, compositeTensor);
// }
// real eval(const uvector3& p) override { return evalPower(compositeTensor, p); }
// bool isInside(const uvector3& p)
// {
// for (auto& t : tensors) {
// if (evalPower(t, p) >= 0) { return false; }
// }
// return true;
// }
// };
// class BernsteinPrimitive : public Primitive
// {
// };
// class BernsteinTensor : public BernsteinPrimitive
// {
// public:
// xarray<real, 3> tensor;
// void print() override { std::cout << "Bernstein" << std::endl; }
// real eval(const uvector3& p) override { return evalBernstein(tensor, p); }
// BernsteinTensor(const PowerTensor& pt_)
// {
// auto v0 = xarray2StdVector(pt_.tensor);
// tensor.ext_ = pt_.tensor.ext();
// algoim_spark_alloc(real, tensor);
// auto v1 = xarray2StdVector(pt_.tensor);
// auto v2 = xarray2StdVector(tensor);
// detail::power2BernsteinTensor(pt_.tensor, tensor);
// uvector3 x(0.2, 0.5, 0.6);
// real BernsteinValue = bernstein::evalBernsteinPoly(tensor, x);
// real PowerValue = evalPower(pt_.tensor, x);
// int a = 1;
// }
// bool isInside(uvector3 p) { return eval(p) < 0; }
// };
// class BernsteinTensorComplex : public BernsteinPrimitive
// {
// public:
// bool fromPower;
// xarray<real, 3> compositeTensor; // 复合后的张量
// std::vector<xarray<real, 3>> tensors; // 原始张量
// void print() override { std::cout << "Bernstein Complex" << std::endl; }
// real eval(const uvector3& p) override { return evalBernstein(compositeTensor, p); }
// // BernsteinTensorComplex(const PowerTensorComplex& pc_) : fromPower(true), tensors(pc_.tensors)
// // {
// // compositeTensor.ext_ = pc_.compositeTensor.ext();
// // algoim_spark_alloc(real, compositeTensor);
// // detail::power2BernsteinTensor(pc_.compositeTensor, compositeTensor);
// // };
// bool isInside(uvector3 p)
// {
// if (fromPower) {
// for (auto& t : tensors) {
// if (evalPower(t, p) >= 0) { return false; }
// }
// return true;
// } else {
// for (auto& t : tensors) {
// if (eval(p) >= 0) { return false; }
// }
// return true;
// }
// return true;
// };
// };
bool isInsidePowers(const std::vector<tensor3>& tensors, const uvector3& p)
{
for (auto& t : tensors) {
if (evalPower(t, p) >= 0) { return false; }
}
return true;
};
bool isInsideBernstein(const tensor3& t, const uvector3& p) { return evalBernstein(t, p) < 0; }
bool isInsidePower(const tensor3& t, const uvector3& p) { return evalPower(t, p) < 0; }
class FRep : public Primitive
{
private:
std::function<real(uvector3)> f;
public:
void print() override { std::cout << "FRep" << std::endl; }
real eval(const uvector3& p) override { return f(p); }
FRep(std::function<real(uvector3)> f_) : f(f_) {}
};
enum PrimitiveType { Sphere, Cylinder, Cone, Mesh, BRep };
class PrimitiveDesc
{
public:
// const static PrimitiveType type;
PrimitiveDesc() = default;
virtual void print() {} // 空定义也可以,但是一定要有定义
};
class ParametricSurface
{
public:
virtual uvector3 eval(uvector2 p) { return uvector3(0, 0, 0); }
};
class BezierSurface : public ParametricSurface
{
private:
std::vector<std::vector<uvector3>> controlPoints;
};
class NURBSSurface : public ParametricSurface
{
std::vector<std::vector<uvector3>> controlPoints;
std::vector<std::vector<real>> weights;
std::vector<real> knotsU, knotsV;
int degreeU, degreeV;
};
class ParametricCurve
{
public:
virtual uvector3 eval(real p) { return uvector3(0, 0, 0); }
};
class BezierCurve : public ParametricCurve
{
private:
std::vector<uvector3> controlPoints;
};
class NURBSCurve : public ParametricCurve
{
private:
std::vector<uvector3> controlPoints;
std::vector<real> weights;
std::vector<real> knots;
int degree;
};
class BRepDesc : virtual public PrimitiveDesc
{
const static PrimitiveType type = BRep;
std::vector<uvector3> vertices;
std::vector<ParametricCurve> curves;
std::vector<ParametricSurface> surfaces;
public:
void print() override { std::cout << "BRep Description" << std::endl; }
real eval(const uvector3& p)
{
// DOTO: the implicit conversion of Parametric BRep
return 0;
}
BRepDesc(const std::vector<uvector3>& vs_,
const std::vector<ParametricCurve>& cs_,
const std::vector<ParametricSurface>& ss_)
: vertices(vs_), curves(cs_), surfaces(ss_)
{
}
BRepDesc(const BRepDesc& brep)
{
vertices = brep.vertices;
curves = brep.curves;
surfaces = brep.surfaces;
}
};
class SphereDesc : virtual public PrimitiveDesc
{
public:
const static PrimitiveType type = Sphere;
real radius;
uvector3 center;
uvector3 amplitude;
SphereDesc(real r_, const uvector3& c_, const uvector3& a_) : PrimitiveDesc(), radius(r_), center(c_), amplitude(a_) {}
void print() override { std::cout << "Sphere Description" << std::endl; }
};
class CylinderDesc : virtual public PrimitiveDesc
{
const static PrimitiveType type = Cylinder;
uvector3 node1;
uvector3 node2;
real radius;
CylinderDesc(const uvector3& n1_, const uvector3& n2_, real r_) : PrimitiveDesc(), node1(n1_), node2(n2_), radius(r_) {}
void print() override { std::cout << "Cylinder Description" << std::endl; }
};
class ConeDesc : virtual public PrimitiveDesc
{
const static PrimitiveType type = Cone;
uvector3 node1;
uvector3 node2;
real radius;
ConeDesc(const uvector3& n1_, const uvector3& n2_, real r_) : PrimitiveDesc(), node1(n1_), node2(n2_), radius(r_) {}
void print() override { std::cout << "Cone Description" << std::endl; }
};
class MeshDesc : virtual public PrimitiveDesc
{
public:
const static PrimitiveType type = Mesh;
std::vector<uvector3> vertices;
std::vector<int> indices;
std::vector<int> indexInclusiveScan;
MeshDesc(const std::vector<uvector3>& vertices_,
const std::vector<int>& indices_,
const std::vector<int>& indexInclusiveScan_)
: PrimitiveDesc(), vertices(vertices_), indices(indices_), indexInclusiveScan(indexInclusiveScan_)
{
}
void print() override { std::cout << "Mesh Description" << std::endl; }
};
void makeMesh(const MeshDesc& mesh, xarray<real, 3>& tensor, std::vector<xarray<real, 3>>& planeTensors)
{
uvector3 ext(1 + mesh.indexInclusiveScan.size());
assert(all(ext == tensor.ext()));
assert(planeTensors.size() == mesh.indexInclusiveScan.size());
// for (const auto& index : indices) {
for (int i = 0; i < mesh.indexInclusiveScan.size(); ++i) {
const int indexBeg = i == 0 ? 0 : mesh.indexInclusiveScan[i - 1];
const int indexSize = mesh.indexInclusiveScan[i] - indexBeg;
assert(indexSize >= 3);
auto& planeTensor = planeTensors[i];
xarrayInit(planeTensor);
auto& vertices = mesh.vertices;
auto& indices = mesh.indices;
uvector3 V01 = vertices[indices[indexBeg + 1]] - mesh.vertices[indices[indexBeg]];
uvector3 V02 = vertices[indices[indexBeg + 2]] - mesh.vertices[indices[indexBeg]];
uvector3 N = cross(V01, V02);
N /= norm(N);
real d = -dot(N, vertices[indices[indexBeg]]);
// 法线所指方向为>0区域
planeTensor.m(uvector3(0, 0, 0)) = d;
planeTensor.m(uvector3(1, 0, 0)) = N(0);
planeTensor.m(uvector3(0, 1, 0)) = N(1);
planeTensor.m(uvector3(0, 0, 1)) = N(2);
// test other vertices
for (int j = indexBeg + 3; j < mesh.indexInclusiveScan[i]; ++j) {
assert(dot(N, vertices[indices[j]]) + d < std::numeric_limits<real>::epsilon());
}
}
detail::compositePower(planeTensors, 0, 0, 1, tensor);
// return PowerTensorComplex(planeTensors, tensor);
};
void makeSphere(const SphereDesc& sphereDesc, xarray<real, 3>& tensor)
{
uvector<int, 3> ext = 3;
assert(all(ext == tensor.ext()));
for (int dim = 0; dim < 3; ++dim) {
uvector<int, 3> idx = 0;
tensor.m(idx) += sphereDesc.center(dim) * sphereDesc.center(dim);
idx(dim) = 1;
tensor.m(idx) = 2 * sphereDesc.amplitude(dim) * (-sphereDesc.center(dim));
idx(dim) = 2;
tensor.m(idx) = sphereDesc.amplitude(dim) * sphereDesc.amplitude(dim);
}
tensor.m(0) -= sphereDesc.radius * sphereDesc.radius;
// return PowerTensor(tensor);
};
void makeCylinder(xarray<real, 3>& tensor, uvector3 startPt, uvector3 endPt, real r)
{
// PowerTensor pt;
// TODO:
}
struct Scene {
std::vector<tensor3> polys;
real* boolDescription; // blobtree
uvector3 min, max;
};
struct Node {
std::vector<int> polyIntersectIndices; // 同一poly会出现在不同node的polyIndices中
std::vector<int> polyFullyContainedIndices; // 完全包含该node的poly, 同样会出现在不同node的polyIndices中
int children[8]; // octree
uvector3 min, max;
};
} // namespace algoim::Organizer