From c19d18a3a652de49a302433ca98fd3788e43c37c Mon Sep 17 00:00:00 2001 From: k-a-z-u Date: Wed, 6 Sep 2017 17:04:19 +0200 Subject: [PATCH] modified lib GPC for header only worked on 3d traytracing --- CMakeLists.txt | 10 +- geo/BBox3.h | 11 +- geo/Point3.h | 4 +- geo/Sphere3.h | 180 ++++++++++++++-------- geo/TestSphere3.cpp | 0 geo/Triangle3.h | 2 +- geo/volume/BVH.h | 224 +++++++++++++++++++++++++++- geo/volume/BVHDebug.h | 106 +++++++++++++ geo/volume/BoundingVolume.h | 18 ++- geo/volume/BoundingVolumeAABB.h | 2 +- geo/volume/BoundingVolumeSphere.h | 48 ++++-- lib/gpc/{gpc.cpp => gpc.cpp.h} | 0 lib/gpc/gpc.h | 45 +++--- main.cpp | 4 +- tests/geo/TestBVH.cpp | 198 ++++++++++++++++++++++++ tests/ray/TestRayTrace3.cpp | 95 ++++++++---- wifi/estimate/ray3/DataMap3.h | 6 + wifi/estimate/ray3/ModelFactory.h | 14 +- wifi/estimate/ray3/ObstacleTree.h | 139 ----------------- wifi/estimate/ray3/WifiRayTrace3D.h | 77 +++++++--- 20 files changed, 884 insertions(+), 299 deletions(-) delete mode 100644 geo/TestSphere3.cpp create mode 100644 geo/volume/BVHDebug.h rename lib/gpc/{gpc.cpp => gpc.cpp.h} (100%) create mode 100644 tests/geo/TestBVH.cpp delete mode 100644 wifi/estimate/ray3/ObstacleTree.h diff --git a/CMakeLists.txt b/CMakeLists.txt index d7b87de..e0e7862 100755 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -43,6 +43,13 @@ FILE(GLOB SOURCES ./*/*/*/*.cpp ) +FIND_PACKAGE( OpenMP REQUIRED) +if(OPENMP_FOUND) + message("OPENMP FOUND") + set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} ${OpenMP_C_FLAGS}") + set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${OpenMP_CXX_FLAGS}") + set(CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS} ${OpenMP_EXE_LINKER_FLAGS}") +endif() if(${CMAKE_GENERATOR} MATCHES "Visual Studio") @@ -68,7 +75,8 @@ ADD_DEFINITIONS( -fstack-protector-all -g3 - #-O0 +# -O0 + -O2 -march=native -DWITH_TESTS diff --git a/geo/BBox3.h b/geo/BBox3.h index eb53afa..ed81f6c 100644 --- a/geo/BBox3.h +++ b/geo/BBox3.h @@ -18,11 +18,16 @@ private: public: - /** empty ctor */ + /** empty ctor */ BBox3() : p1(MAX,MAX,MAX), p2(MIN,MIN,MIN) {;} - /** ctor with min and max */ - BBox3(const Point3 min, const Point3 max) : p1(min), p2(max) {;} + /** ctor with min and max */ + BBox3(const Point3 min, const Point3 max) : p1(min), p2(max) {;} + + /** create a bbox around the given point */ + static BBox3 around(const Point3 center, const Point3 size) { + return BBox3(center-size/2, center+size/2); + } /** adjust the bounding-box by adding this point */ void add(const Point3& p) { diff --git a/geo/Point3.h b/geo/Point3.h index 6038b0b..f436d38 100644 --- a/geo/Point3.h +++ b/geo/Point3.h @@ -116,11 +116,11 @@ private: }; -inline float dot(const Point3& p1, const Point3& p2) { +inline float dot(const Point3 p1, const Point3 p2) { return (p1.x*p2.x) + (p1.y*p2.y) + (p1.z*p2.z); } -inline Point3 cross(const Point3& a, const Point3& b) { +inline Point3 cross(const Point3 a, const Point3 b) { return Point3( a.y*b.z - a.z*b.y, a.z*b.x - a.x*b.z, diff --git a/geo/Sphere3.h b/geo/Sphere3.h index e272e9b..ed2bfda 100644 --- a/geo/Sphere3.h +++ b/geo/Sphere3.h @@ -19,38 +19,11 @@ struct Sphere3 { /** ctor */ Sphere3(const Point3 center, const float radius) : center(center), radius(radius) {;} - /** create a sphere that fully contains the given point-set */ - static Sphere3 around(const std::vector& lst) { - // NOT OPTIMAL but fast - Point3 sum(0,0,0); - for (const Point3 p : lst) {sum += p;} - const Point3 center = sum / lst.size(); - float radius = 0; - for (const Point3 p : lst) { - const float dist = center.getDistance(p); - if (dist > radius) {radius = dist;} - } - return Sphere3(center, radius); - - } - - - - - static Sphere3 join(const Sphere3& a, const Sphere3& b) { - - if (a.contains(b)) {return a;} - if (b.contains(a)) {return b;} - - Point3 newCen = (a.center + b.center) / 2; - float newRad = (a.center.getDistance(b.center) + std::max(a.radius, b.radius) * 2) / 2; - return Sphere3(newCen, newRad); - } /** does this sphere contain the given sphere? */ bool contains(const Sphere3& o) const { - return (o.center.getDistance(this->center) + o.radius) < this->radius; + return (o.center.getDistance(this->center) + o.radius) <= this->radius; } /** does the sphere contain the given point? */ @@ -58,13 +31,21 @@ struct Sphere3 { return center.getDistance(p) <= radius; } + /** does the sphere intersect with the given sphere? */ + bool intersects(const Sphere3& other) const { + return center.getDistance(other.center) < (radius + other.radius); + } + + /** does the sphere intersect with the given ray? */ bool intersects(const Ray3& ray) const { - if (contains(ray.start)) {return true;} + // if the sphere contains the ray's start -> done + //if (contains(ray.start)) {return true;} + - /* // https://stackoverflow.com/questions/6533856/ray-sphere-intersection + /* const float xA = ray.start.x; const float yA = ray.start.y; const float zA = ray.start.z; @@ -85,42 +66,121 @@ struct Sphere3 { // intersection? return delta >= 0; - */ + */ // http://www.ccs.neu.edu/home/fell/CSU540/programs/RayTracingFormulas.htm - - const double x0 = ray.start.x; - const double y0 = ray.start.y; - const double z0 = ray.start.z; - - const double cx = center.x; - const double cy = center.y; - const double cz = center.z; - - const double dx = ray.dir.x; - const double dy = ray.dir.y; - const double dz = ray.dir.z; - - const double a = dx*dx + dy*dy + dz*dz; - const double b = 2*dx*(x0-cx) + 2*dy*(y0-cy) + 2*dz*(z0-cz); - const double c = cx*cx + cy*cy + cz*cz + x0*x0 + y0*y0 + z0*z0 + -2*(cx*x0 + cy*y0 + cz*z0) - radius*radius; - - const double discriminant = (b*b) - (4*a*c); - return discriminant >= 0; - - /* - // http://www.pixelnerve.com/v/2009/02/08/ray-sphere-intersection/ - const float a = ray.dir.length(); - //if (a == 0.0) return 0; - const float b = 2.0f * ( dot(ray.start, ray.dir) - dot(ray.dir, center)) ; - const Point3 tempDiff = center - ray.start; - const float c = tempDiff.length() - (radius*radius); - const float disc = b * b - 4 * a * c; - return disc >= 0.0f; + const float x0 = ray.start.x; + const float y0 = ray.start.y; + const float z0 = ray.start.z; + + const float cx = center.x; + const float cy = center.y; + const float cz = center.z; + + const float dx = ray.dir.x; + const float dy = ray.dir.y; + const float dz = ray.dir.z; + + const float a = dx*dx + dy*dy + dz*dz; + const float b = 2*dx*(x0-cx) + 2*dy*(y0-cy) + 2*dz*(z0-cz); + const float c = cx*cx + cy*cy + cz*cz + x0*x0 + y0*y0 + z0*z0 + -2*(cx*x0 + cy*y0 + cz*z0) - radius*radius; + + const float discriminant = (b*b) - (4*a*c); + return discriminant >= 0; */ + + //https://gamedev.stackexchange.com/questions/96459/fast-ray-sphere-collision-code + + const Point3 m = ray.start - center; + const float b = dot(m, ray.dir); + const float c = dot(m, m) - radius*radius; + + // Exit if r’s origin outside s (c > 0) and r pointing away from s (b > 0) + if (c > 0.0f && b > 0.0f) {return false;} + const float discr = b*b - c; + + // A negative discriminant corresponds to ray missing sphere + if (discr < 0.0f) {return false;} + + return true; + +// // Ray now found to intersect sphere, compute smallest t value of intersection +// t = -b - Sqrt(discr); + +// // If t is negative, ray started inside sphere so clamp t to zero +// if (t < 0.0f) t = 0.0f; +// q = p + t * d; + +// return 1; + + + } + + /** configure this sphere to contain the given point-set */ + void adjustToPointSet(const std::vector& lst) { + + // NOT OPTIMAL but fast + + // calculate the point set's center + Point3 sum(0,0,0); + for (const Point3 p : lst) {sum += p;} + const Point3 center = sum / lst.size(); + + // calculate the sphere's radius (maximum distance from center + float radius = 0; + for (const Point3 p : lst) { + const float dist = center.getDistance(p); + if (dist > radius) {radius = dist;} + } + + this->radius = radius; + this->center = center; + + } + +public: + + /** create a sphere that fully contains the given point-set */ + static Sphere3 around(const std::vector& lst) { + Sphere3 sphere; + sphere.adjustToPointSet(lst); + return sphere; + } + + /** combine two spheres into a new one containing both */ + static Sphere3 join(const Sphere3& a, const Sphere3& b) { + + // if one already contains the other, just return it as-is + if (a.contains(b)) {return a;} + if (b.contains(a)) {return b;} + + // calculate the new center and radius +// Point3 newCen = (a.center + b.center) / 2; +// float newRad = (a.center.getDistance(b.center) + std::max(a.radius, b.radius) * 2) / 2; +// return Sphere3(newCen, newRad); + +// Point3 newCen = (a.center*a.radius + b.center*b.radius) / (a.radius+b.radius); +// float newRad = (a.center.getDistance(b.center) + a.radius + b.radius) / 2 * 1.02f; // slightly larger to prevent rounding issues +// Sphere3 res(newCen, newRad); + + // create both maximum ends + const Point3 out1 = a.center + (a.center-b.center).normalized() * a.radius; + const Point3 out2 = b.center + (b.center-a.center).normalized() * b.radius; + + // center is within both ends, so is the radius + Point3 newCen = (out1 + out2) / 2; + float newRad = out1.getDistance(out2) / 2 * 1.001; // slightly larger to prevent rounding issues + Sphere3 res(newCen, newRad); + + // check + Assert::isTrue(res.contains(a), "sphere joining error. base-spheres are not contained."); + Assert::isTrue(res.contains(b), "sphere joining error. base-spheres are not contained."); + + return res; + } }; diff --git a/geo/TestSphere3.cpp b/geo/TestSphere3.cpp deleted file mode 100644 index e69de29..0000000 diff --git a/geo/Triangle3.h b/geo/Triangle3.h index b70a239..a83b142 100644 --- a/geo/Triangle3.h +++ b/geo/Triangle3.h @@ -41,7 +41,7 @@ public: // http://www.lighthouse3d.com/tutorials/maths/ray-triangle-intersection/ - bool intersects(Ray3 ray, Point3& pos) const { + bool intersects(const Ray3& ray, Point3& pos) const { const Point3 e1 = p2-p1; const Point3 e2 = p3-p1; diff --git a/geo/volume/BVH.h b/geo/volume/BVH.h index 0333f40..180f3d9 100644 --- a/geo/volume/BVH.h +++ b/geo/volume/BVH.h @@ -1,17 +1,235 @@ #ifndef BOUNDINGVOLUMEHIERARCHY_H #define BOUNDINGVOLUMEHIERARCHY_H +#include +#include + #include "BoundingVolume.h" #include "BoundingVolumeAABB.h" #include "BoundingVolumeSphere.h" -class BVH { + + + +template class BVH { + +protected: + + /** one node within the tree */ + struct BVHNode { + bool isLeaf = true; + Volume boundingVolume; + std::vector childNodes; + BVHNode(bool isLeaf = false) : isLeaf(isLeaf) {;} + }; + + /** one leaf within the tree */ + struct BVHLeaf : public BVHNode { + Element element; + BVHLeaf(const Element& e) : BVHNode(true), element(e) {;} + }; + + /** the tree's root */ + BVHNode root; public: - /** add a new volume to the tree */ - void add(BoundingVolume* bv) { + /** get the tree's root node */ + const BVHNode& getRoot() const { + return root; + } + /** add a new volume to the tree */ + void add(const Element& element) { + + // create a new leaf for this element + BVHLeaf* leaf = new BVHLeaf(element); + + // get the element's boundin volume + leaf->boundingVolume = getBoundingVolume(element); + + // add the leaf to the tree + root.childNodes.push_back(leaf); + + } + + /** optimize the tree */ + int optimize(const int max = 9999) { + for (int i = 0; i < max; ++i) { + //const bool did = concat(); // faster + const bool did = combineBest(); // better + if (!did) {return i;} + } + return max; + } + + void getHits(const Ray3 ray, std::function func) const { + //int tests = 0; int leafs = 0; + getHits(ray, &root, func); + //std::cout << tests << " " << leafs << std::endl; + } + + void getHits(const Ray3 ray, const BVHNode* node, std::function func) const { + for (const BVHNode* sub : node->childNodes) { + if (sub->boundingVolume.intersects(ray)) { + if (sub->isLeaf) { + BVHLeaf* leaf = (BVHLeaf*)(sub); // TODO: cast + func(leaf->element); + } else { + getHits(ray, sub, func); + } + } + } + } + +private: + + bool combineBest() { + + // nothing to do? + if (root.childNodes.size() < 2) {return false;} + + struct Best { + BVHNode* n1 = nullptr; + BVHNode* n2 = nullptr; + Volume vol; + float volSize = 99999999; + } best; + + for (size_t i = 0; i < root.childNodes.size(); ++i) { + for (size_t j = 0; j < root.childNodes.size(); ++j) { + + if (i == j) {continue;} + + BVHNode* n1 = root.childNodes[i]; + BVHNode* n2 = root.childNodes[j]; + + const Volume newVol = Volume::join(n1->boundingVolume, n2->boundingVolume); + const float newVolSize = newVol.getVolumeSize(); + if (newVolSize < best.volSize) { + best.vol = newVol; + best.volSize = newVolSize; + best.n1 = n1; + best.n2 = n2; + } + + } + } + + root.childNodes.erase(std::remove(root.childNodes.begin(), root.childNodes.end(), best.n1), root.childNodes.end()); + root.childNodes.erase(std::remove(root.childNodes.begin(), root.childNodes.end(), best.n2), root.childNodes.end()); + + // combine both into a new node + BVHNode* newNode = new BVHNode(); + newNode->childNodes.push_back(best.n1); + newNode->childNodes.push_back(best.n2); + newNode->boundingVolume = best.vol; + + // does the newly created node contain any other nodes? + // THIS SHOULD NEVER BE THE CASE! +// for (size_t i = 0; i < root.childNodes.size(); ++i) { +// BVHNode* n3 = root.childNodes[i]; +// if (newNode->boundingVolume.contains(n3->boundingVolume)) { +// newNode->childNodes.push_back(n3); +// root.childNodes.erase(root.childNodes.begin()+i); +// --i; +// } +// } + + // attach the node + root.childNodes.push_back(newNode); + + return true; + + } + + + bool concat() { + + // nothing to do? + if (root.childNodes.size() < 2) {return false;} + + + bool concated = false; + + // first, sort all elements by volume (smallest first) + auto compVolume = [] (const BVHNode* n1, const BVHNode* n2) { + return n1->boundingVolume.getVolumeSize() < n2->boundingVolume.getVolumeSize(); + }; + std::sort(root.childNodes.begin(), root.childNodes.end(), compVolume); + + + // elements will be grouped into this new root + BVHNode newRoot; + + // combine nearby elements + while(true) { + + // get [and remove] the next element + BVHNode* n0 = (BVHNode*) root.childNodes[0]; + root.childNodes.erase(root.childNodes.begin()+0); + + // find another element that yields minimal increase in volume + auto compNear = [n0] (const BVHNode* n1, const BVHNode* n2) { + const float v1 = Volume::join(n0->boundingVolume, n1->boundingVolume).getVolumeSize(); + const float v2 = Volume::join(n0->boundingVolume, n2->boundingVolume).getVolumeSize(); + return v1 < v2; + }; + auto it = std::min_element(root.childNodes.begin(), root.childNodes.end(), compNear); + BVHNode* n1 = *it; + + // calculate the resulting increment in volume + const Volume joined = Volume::join(n0->boundingVolume, n1->boundingVolume); + const float increment = joined.getVolumeSize() / n0->boundingVolume.getVolumeSize(); + const bool intersects = n0->boundingVolume.intersects(n1->boundingVolume); + + const bool combine = true; //(intersects); //(increment < 15.0); + + if (combine) { + + // remove from current root + root.childNodes.erase(it); + + // combine both into a new node + BVHNode* node = new BVHNode(); + node->childNodes.push_back(n0); + node->childNodes.push_back(n1); + node->boundingVolume = joined; + newRoot.childNodes.push_back(node); + + concated = true; + + } else { + + BVHNode* node = new BVHNode(); + node->childNodes.push_back(n0); + node->boundingVolume = n0->boundingVolume; + newRoot.childNodes.push_back(node); + + } + + // done? + if (root.childNodes.size() == 1) { + BVHNode* node = new BVHNode(); + node->childNodes.push_back(root.childNodes.front()); + node->boundingVolume = root.childNodes.front()->boundingVolume; + newRoot.childNodes.push_back(node); + break; + } else if (root.childNodes.size() == 0) { + break; + } + + } + + root = newRoot; + return concated; + + } + + /** get a bounding-volume for the given element */ + Volume getBoundingVolume(const Element& element) { + const std::vector verts = Wrapper::getVertices(element); + return Volume::fromVertices(verts); } diff --git a/geo/volume/BVHDebug.h b/geo/volume/BVHDebug.h new file mode 100644 index 0000000..e247462 --- /dev/null +++ b/geo/volume/BVHDebug.h @@ -0,0 +1,106 @@ +#ifndef BVHDEBUG_H +#define BVHDEBUG_H + +#include "BVH.h" +#include +#include +#include +#include +#include +#include "../BBox3.h" +#include + +/** adds some debug helpers to the BVH */ +template class BVHDebug : public BVH { + + using BVHNode = typename BVH::BVHNode; + using BVHLeaf = typename BVH::BVHLeaf; + +public: + +// std::vecto colors { +// "#888888", "#888800", "#008888", "#880088", "#ee0000", "#00ee00", "#0000ee" +// }; + + void show(int maxPts = 1500, bool showLeafs = true) { + + std::stringstream out; + + static K::Gnuplot gp; + K::GnuplotSplot plot; + K::GnuplotSplotElementColorPoints pVol; plot.add(&pVol); //pVol.setColor(K::GnuplotColor::fromRGB(128,128,128)); + K::GnuplotSplotElementPoints pElemPoints; plot.add(&pElemPoints); pElemPoints.setColor(K::GnuplotColor::fromRGB(0,0,255)); + K::GnuplotSplotElementLines pElemLines; plot.add(&pElemLines); pElemLines.getStroke().setColor(K::GnuplotColor::fromRGB(0,0,255)); + + const int depth = recurse(maxPts, showLeafs, 0, &this->root, pVol, pElemPoints, pElemLines); + + plot.getAxisCB().setRange(0, depth); + + gp << "set view equal xyz\n"; + gp.draw(plot); + gp.flush(); + + } + +private: + + int recurse(int maxPts, bool showLeafs, int curDepth, const BVHNode* node, K::GnuplotSplotElementColorPoints& vol, K::GnuplotSplotElementPoints& pElemPoints, K::GnuplotSplotElementLines& elemLines) { + + int resDepth = curDepth; + + for (BVHNode* sub : node->childNodes) { + resDepth = recurse(maxPts, showLeafs, curDepth+1, sub, vol, pElemPoints, elemLines); + } + + if (!node->isLeaf || showLeafs) { + + const int numPts = maxPts / (curDepth+1); + + for (int i = 0; i < numPts; ++i) { + const Point3 p = getRandomPoint(node->boundingVolume); + vol.add(K::GnuplotPoint3(p.x, p.y, p.z), curDepth); + } + + } + + if (node->isLeaf) { + BVHLeaf* leaf = (BVHLeaf*) node; + std::vector verts = Wrapper::getVertices(leaf->element); + for (const Point3 p : verts) { + pElemPoints.add(K::GnuplotPoint3(p.x, p.y, p.z)); + } + addLines(leaf->element, elemLines); + } + + return resDepth; + + } + + Point3 getRandomPoint(BoundingVolumeSphere sphere) { + static std::minstd_rand gen; + std::uniform_real_distribution dist(-1, +1); + Point3 dir = Point3(dist(gen), dist(gen), dist(gen)).normalized() * sphere.radius; + return sphere.center + dir; + } + + void addLines(const Element& elem, K::GnuplotSplotElementLines& elemLines) { + + std::vector pts = Wrapper::getDebugLines(elem); + + for (size_t i = 0; i< pts.size(); i+=2) { + + const Point3 p1 = pts[i+0]; + const Point3 p2 = pts[i+1]; + + K::GnuplotPoint3 gp1(p1.x, p1.y, p1.z); + K::GnuplotPoint3 gp2(p2.x, p2.y, p2.z); + + elemLines.addSegment(gp1, gp2); + + } + + } + +}; + +#endif // BVHDEBUG_H diff --git a/geo/volume/BoundingVolume.h b/geo/volume/BoundingVolume.h index d259895..80516a3 100644 --- a/geo/volume/BoundingVolume.h +++ b/geo/volume/BoundingVolume.h @@ -2,16 +2,26 @@ #define BOUNDINGVOLUME_H #include "../Point3.h" +#include "../Ray3.h" class BoundingVolume { public: - /** get the volume's size (something like m^3) */ - virtual float getVolumeSize() const = 0; +// /** get the volume's size (something like m^3) */ +// virtual float getVolumeSize() const = 0; - /** does the volume contain the given point? */ - virtual bool contains(const Point3 p) const = 0; +// /** does the volume contain the given point? */ +// virtual bool contains(const Point3 p) const = 0; + +// /** does the volume contain the given volume? */ +// virtual bool contains(const BoundingVolume& other) const = 0; + +// /** does the volume intersect with the given ray? */ +// virtual bool intersects(const Ray3& ray) const = 0; + +// /** does the volume intersect with the given volume? */ +// virtual bool intersects(const BoundingVolume& other) const = 0; }; diff --git a/geo/volume/BoundingVolumeAABB.h b/geo/volume/BoundingVolumeAABB.h index b99c666..c5533bc 100644 --- a/geo/volume/BoundingVolumeAABB.h +++ b/geo/volume/BoundingVolumeAABB.h @@ -20,7 +20,7 @@ public: /** empty ctor */ BoundingVolumeAABB() : p1(MAX,MAX,MAX), p2(MIN,MIN,MIN) {;} - float getVolumeSize() const override { + float getVolumeSize() const { return (p2.x-p1.x) * (p2.y-p1.y) * (p2.z-p1.z); } diff --git a/geo/volume/BoundingVolumeSphere.h b/geo/volume/BoundingVolumeSphere.h index 22345a6..83a71f4 100644 --- a/geo/volume/BoundingVolumeSphere.h +++ b/geo/volume/BoundingVolumeSphere.h @@ -1,25 +1,51 @@ #ifndef BOUNDINGVOLUMESPHERE_H #define BOUNDINGVOLUMESPHERE_H +#include "BoundingVolume.h" +#include "../Sphere3.h" #include "../Point3.h" -class BoundingVolumeSphere : public BoundingVolume { - -private: - - Point3 center; - - float radius; +class BoundingVolumeSphere : public BoundingVolume, public Sphere3 { public: - float getVolumeSize() const override { + BoundingVolumeSphere() {;} + + BoundingVolumeSphere(const Sphere3& s) : Sphere3(s) {;} + + float getVolumeSize() const { return (4.0f / 3.0f) * M_PI * (radius*radius*radius); } - /** does the volume contain the given point? */ - virtual bool contains(const Point3 p) const { - return (center-p).length() <= radius; + bool contains(const Point3 p) const { + return Sphere3::contains(p); + } + + bool intersects(const Ray3& ray) const { + return Sphere3::intersects(ray); + } + + /** does the volume intersect with the given volume? */ + bool intersects(const BoundingVolume& other) const { + const BoundingVolumeSphere& sphere = (const BoundingVolumeSphere&) other; + return Sphere3::intersects(sphere); + } + + /** does the volume contain the given volume? */ + bool contains(const BoundingVolume& other) const { + const BoundingVolumeSphere& sphere = (const BoundingVolumeSphere&) other; + return Sphere3::contains(sphere); + } + + /** construct a volume around the given point-set */ + static BoundingVolumeSphere fromVertices(const std::vector& verts) { + BoundingVolumeSphere bvs; + bvs.adjustToPointSet(verts); + return bvs; + } + + static BoundingVolumeSphere join(const BoundingVolumeSphere a, const BoundingVolumeSphere b) { + return BoundingVolumeSphere(Sphere3::join(a, b)); } }; diff --git a/lib/gpc/gpc.cpp b/lib/gpc/gpc.cpp.h similarity index 100% rename from lib/gpc/gpc.cpp rename to lib/gpc/gpc.cpp.h diff --git a/lib/gpc/gpc.h b/lib/gpc/gpc.h index a975b61..309c3f8 100755 --- a/lib/gpc/gpc.h +++ b/lib/gpc/gpc.h @@ -95,34 +95,36 @@ typedef struct /* Tristrip set structure */ =========================================================================== */ -void gpc_read_polygon (FILE *infile_ptr, - int read_hole_flags, - gpc_polygon *polygon); +inline void gpc_read_polygon (FILE *infile_ptr, + int read_hole_flags, + gpc_polygon *polygon); -void gpc_write_polygon (FILE *outfile_ptr, - int write_hole_flags, - gpc_polygon *polygon); +inline void gpc_write_polygon (FILE *outfile_ptr, + int write_hole_flags, + gpc_polygon *polygon); -void gpc_add_contour (gpc_polygon *polygon, - gpc_vertex_list *contour, - int hole); +inline void gpc_add_contour (gpc_polygon *polygon, + gpc_vertex_list *contour, + int hole); -void gpc_polygon_clip (gpc_op set_operation, - gpc_polygon *subject_polygon, - gpc_polygon *clip_polygon, - gpc_polygon *result_polygon); +inline void gpc_polygon_clip (gpc_op set_operation, + gpc_polygon *subject_polygon, + gpc_polygon *clip_polygon, + gpc_polygon *result_polygon); -void gpc_tristrip_clip (gpc_op set_operation, - gpc_polygon *subject_polygon, - gpc_polygon *clip_polygon, - gpc_tristrip *result_tristrip); +inline void gpc_tristrip_clip (gpc_op set_operation, + gpc_polygon *subject_polygon, + gpc_polygon *clip_polygon, + gpc_tristrip *result_tristrip); -void gpc_polygon_to_tristrip (gpc_polygon *polygon, - gpc_tristrip *tristrip); +inline void gpc_polygon_to_tristrip (gpc_polygon *polygon, + gpc_tristrip *tristrip); -void gpc_free_polygon (gpc_polygon *polygon); +inline void gpc_free_polygon (gpc_polygon *polygon); -void gpc_free_tristrip (gpc_tristrip *tristrip); +inline void gpc_free_tristrip (gpc_tristrip *tristrip); + +#include "gpc.cpp.h" #endif @@ -131,3 +133,4 @@ void gpc_free_tristrip (gpc_tristrip *tristrip); End of file: gpc.h =========================================================================== */ + diff --git a/main.cpp b/main.cpp index 4fe22c9..57e1bae 100755 --- a/main.cpp +++ b/main.cpp @@ -33,10 +33,12 @@ int main(int argc, char** argv) { //::testing::GTEST_FLAG(filter) = "*Triangle*"; //::testing::GTEST_FLAG(filter) = "*Ray.ModelFac*"; //::testing::GTEST_FLAG(filter) = "*DataMap3*"; - ::testing::GTEST_FLAG(filter) = "*RayTrace3*"; //::testing::GTEST_FLAG(filter) = "*Matrix4*"; //::testing::GTEST_FLAG(filter) = "*Sphere3*"; + ::testing::GTEST_FLAG(filter) = "*RayTrace3*"; + //::testing::GTEST_FLAG(filter) = "*BVH*"; + //::testing::GTEST_FLAG(filter) = "*Barometer*"; //::testing::GTEST_FLAG(filter) = "*GridWalk2RelPressure*"; diff --git a/tests/geo/TestBVH.cpp b/tests/geo/TestBVH.cpp new file mode 100644 index 0000000..83acca2 --- /dev/null +++ b/tests/geo/TestBVH.cpp @@ -0,0 +1,198 @@ +#ifdef WITH_TESTS + +#include "../Tests.h" + +#include "../../geo/volume/BVH.h" +#include "../../geo/volume/BVHDebug.h" +#include "../../geo/BBox3.h" + +#include "../../floorplan/v2/Floorplan.h" +#include "../../floorplan/v2/FloorplanReader.h" +#include "../../wifi/estimate/ray3/ModelFactory.h" + +struct Wrapper { + + static std::vector getVertices(const BBox3& bbox) { + return {bbox.getMin(), bbox.getMax()}; + } + + static std::vector getDebugLines(const BBox3& bbox) { + + Point3 p1(bbox.getMin().x, bbox.getMin().y, bbox.getMin().z); + Point3 p2(bbox.getMax().x, bbox.getMin().y, bbox.getMin().z); + Point3 p3(bbox.getMax().x, bbox.getMax().y, bbox.getMin().z); + Point3 p4(bbox.getMin().x, bbox.getMax().y, bbox.getMin().z); + + Point3 p5(bbox.getMin().x, bbox.getMin().y, bbox.getMax().z); + Point3 p6(bbox.getMax().x, bbox.getMin().y, bbox.getMax().z); + Point3 p7(bbox.getMax().x, bbox.getMax().y, bbox.getMax().z); + Point3 p8(bbox.getMin().x, bbox.getMax().y, bbox.getMax().z); + + std::vector res; + res.push_back(p1); res.push_back(p2); + res.push_back(p2); res.push_back(p3); + res.push_back(p3); res.push_back(p4); + res.push_back(p4); res.push_back(p1); + + res.push_back(p5); res.push_back(p6); + res.push_back(p6); res.push_back(p7); + res.push_back(p7); res.push_back(p8); + res.push_back(p8); res.push_back(p5); + + res.push_back(p1); res.push_back(p5); + res.push_back(p2); res.push_back(p6); + res.push_back(p3); res.push_back(p7); + res.push_back(p4); res.push_back(p8); + return res; + + } + +}; + +struct WrapperObs3D { + + static std::vector getVertices(const Obstacle3D& o) { + std::vector pts; + for (const Triangle3& tria : o.triangles) { + pts.push_back(tria.p1); + pts.push_back(tria.p2); + pts.push_back(tria.p3); + } + return pts; + } + + static std::vector getDebugLines(const Obstacle3D& o) { + std::vector pts; + for (const Triangle3& tria : o.triangles) { + pts.push_back(tria.p1); pts.push_back(tria.p2); + pts.push_back(tria.p2); pts.push_back(tria.p3); + pts.push_back(tria.p3); pts.push_back(tria.p1); + } + return pts; + } + +}; + + + +TEST(BVH, tree) { + + BVHDebug tree; + + BBox3 bb1(Point3(0,0,0), Point3(1,1,1)); + tree.add(bb1); + + BBox3 bb2(Point3(-1,-1,-1), Point3(0,0,0)); + tree.add(bb2); + + tree.optimize(); + tree.show(); + + int i = 0; (void) i; + +} + +TEST(BVH, tree2) { + + BVHDebug tree; + + BBox3 bb1(Point3(0,0,0), Point3(1,1,1)); + tree.add(bb1); + + BBox3 bb2(Point3(-1,0,0), Point3(0,1,1)); + tree.add(bb2); + + tree.optimize(); + tree.show(); + + int i = 0; (void) i; + +} + +TEST(BVH, tree3) { + + BVHDebug tree; + + BBox3 bb1 = BBox3::around(Point3(+0.5, +0.5, 0.0), Point3(0.25, 0.25, 0.25)); + tree.add(bb1); + + BBox3 bb2 = BBox3::around(Point3(-0.5, +0.5, 0.0), Point3(0.25, 0.25, 0.25)); + tree.add(bb2); + + BBox3 bb3 = BBox3::around(Point3(-0.0, +0.5, 0.0), Point3(0.36, 0.36, 0.36)); + tree.add(bb3); + + BBox3 bb4 = BBox3::around(Point3(-0.0, +0.0, 0.0), Point3(0.5, 0.5, 0.5)); + tree.add(bb4); + + tree.optimize(1); + tree.show(); + + tree.optimize(1); + tree.show(); + + tree.optimize(1); + tree.show(); + + int i = 0; (void) i; + +} + +TEST(BVH, treeMap) { + + std::string file = "/apps/SHL39.xml"; + Floorplan::IndoorMap* map = Floorplan::Reader::readFromFile(file); + + ModelFactory fac(map); + fac.setExportCeilings(false); + fac.setFloors({map->floors[3]}); + std::vector obs = fac.triangulize(); + + BVHDebug tree; + + for (const Obstacle3D& o : obs) { + tree.add(o); + } + + //tree.show(150); + + //int rounds = tree.optimize(); + + for (int i = 0; i < 200; ++i) { + tree.optimize(1); + //if (i%3==0) { + tree.show(250, false); + //} + } + + int i = 0; (void) i; + +} + +TEST(BVH, treeRandom) { + + BVHDebug tree; + + std::minstd_rand gen; + std::uniform_real_distribution dPos(-4.0, +4.0); + std::uniform_real_distribution dSize(+0.3, +1.0); + + for (int i = 0; i < 50; ++i) { + const Point3 pos(dPos(gen), dPos(gen), dPos(gen)); + const Point3 size(dSize(gen), dSize(gen), dSize(gen)); + BBox3 bb = BBox3::around(pos, size); + tree.add(bb); + } + + tree.show(); + +// for (int i = 0; i < 25; ++i) { +// tree.optimize(1); +// tree.show(); +// } + + int i = 0; (void) i; + +} + +#endif diff --git a/tests/ray/TestRayTrace3.cpp b/tests/ray/TestRayTrace3.cpp index edccf3c..55cee92 100644 --- a/tests/ray/TestRayTrace3.cpp +++ b/tests/ray/TestRayTrace3.cpp @@ -12,54 +12,87 @@ TEST(RayTrace3, test) { //Floorplan::IndoorMap* map = Floorplan::Reader::readFromFile(file); //Floorplan::AccessPoint* ap = map->floors[0]->accesspoints[0]; - std::string file = "/mnt/data/workspaces/IndoorMap/maps/SHL39.xml"; + std::string file = "/apps/SHL39.xml"; Floorplan::IndoorMap* map = Floorplan::Reader::readFromFile(file); Floorplan::AccessPoint* ap = map->floors[0]->accesspoints[4]; ModelFactory fac(map); - std::ofstream outOBJ("/mnt/vm/map.obj"); + std::ofstream outOBJ("/tmp/vm/map.obj"); outOBJ << fac.toOBJ(); outOBJ.close(); - const int gs_cm = 100; + const int gs_cm = 50; WiFiRaytrace3D rt(map, gs_cm, ap->pos); + + std::chrono::time_point start = std::chrono::high_resolution_clock::now(); const DataMap3Signal& dms = rt.estimate(); + std::chrono::time_point end = std::chrono::high_resolution_clock::now(); + auto result = std::chrono::duration_cast(end-start).count(); + std::cout << "it took: " << result << " msec" << std::endl; - const char sep = ';'; + if (1 == 1) { - std::ofstream out("/mnt/vm/rays.xyz.txt"); - auto lambda = [&] (const float x, const float y, const float z, const DataMap3SignalEntry& e) { + const char sep = ' '; + int lines = 0; - const float min = -120; - const float max = -40; - float rssi = e.getMaxRSSI(); - if (rssi > max) {rssi = max;} + std::stringstream tmp; - if (rssi > -100) { - const float v = ((rssi - min) / (max-min)) * 255; - out - << x << sep << y << sep << z << sep - << v << sep << v << sep << v - << "\n"; + auto lambda = [&] (const float x, const float y, const float z, const DataMap3SignalEntry& e) { + + const float min = -120; + const float max = -40; + float rssi = e.getMaxRSSI(); + if (rssi > max) {rssi = max;} + + if (z < 1.0 || z > 1.0) {return;} + + if (rssi > -100) { + ++lines; + const int v = ((rssi - min) / (max-min)) * 255; // color + tmp + << x << sep << y << sep << z << sep + << v << sep << v << sep << v + << "\n"; + } + }; + + + dms.forEach(lambda); + + std::ofstream out("/tmp/vm/grid.ply"); + out << "ply\n"; + out << "format ascii 1.0\n"; + out << "comment VCGLIB generated\n"; + out << "element vertex " << lines << "\n"; + out << "property float x\n"; + out << "property float y\n"; + out << "property float z\n"; + out << "property uchar red\n"; + out << "property uchar green\n"; + out << "property uchar blue\n"; + out << "element face 0\n"; + out << "property list uchar int vertex_indices\n"; + out << "end_header\n"; + out << tmp.str(); + out.close(); + + std::cout << "lines: " << lines << std::endl; + + + std::ofstream outHits("/tmp/vm/hits.xyz.txt"); + for (const Point3 hit : rt.getHitEnter()) { + outHits << hit.x << sep << hit.y << sep << hit.z << sep << 0 << sep << 255 << sep << 0 << "\n"; } - }; + for (const Point3 hit : rt.getHitLeave()) { + outHits << hit.x << sep << hit.y << sep << hit.z << sep << 0 << sep << 0 << sep << 255 << "\n"; + } + for (const Point3 hit : rt.getHitStop()) { + outHits << hit.x << sep << hit.y << sep << hit.z << sep << 0 << sep << 0 << sep << 0 << "\n"; + } + outHits.close(); - - dms.forEach(lambda); - out.close(); - - std::ofstream outHits("/mnt/vm/hits.xyz.txt"); - for (const Point3 hit : rt.getHitEnter()) { - outHits << hit.x << sep << hit.y << sep << hit.z << sep << 0 << sep << 255 << sep << 0 << "\n"; } - for (const Point3 hit : rt.getHitLeave()) { - outHits << hit.x << sep << hit.y << sep << hit.z << sep << 0 << sep << 0 << sep << 255 << "\n"; - } - for (const Point3 hit : rt.getHitStop()) { - outHits << hit.x << sep << hit.y << sep << hit.z << sep << 0 << sep << 0 << sep << 0 << "\n"; - } - outHits.close(); } diff --git a/wifi/estimate/ray3/DataMap3.h b/wifi/estimate/ray3/DataMap3.h index f72106a..cc0d8ea 100644 --- a/wifi/estimate/ray3/DataMap3.h +++ b/wifi/estimate/ray3/DataMap3.h @@ -5,6 +5,7 @@ #include "../../../geo/BBox3.h" #include #include +#include template class DataMap3 { @@ -189,6 +190,8 @@ private: struct DataMap3SignalEntry { + + struct Entry { float rssi; float distanceToAP; @@ -198,8 +201,11 @@ struct DataMap3SignalEntry { std::vector entries; void add(const float rssi, const float distanceToAP) { + static std::mutex mtx; Entry e(rssi, distanceToAP); + mtx.lock(); entries.push_back(e); + mtx.unlock(); } float getMaxRSSI() const { diff --git a/wifi/estimate/ray3/ModelFactory.h b/wifi/estimate/ray3/ModelFactory.h index 91f0b83..11a6390 100644 --- a/wifi/estimate/ray3/ModelFactory.h +++ b/wifi/estimate/ray3/ModelFactory.h @@ -17,6 +17,7 @@ private: bool exportCeilings = true; bool exportObstacles = true; bool exportWallTops = false; + std::vector exportFloors; const Floorplan::IndoorMap* map; @@ -28,14 +29,25 @@ public: } + void setExportCeilings(bool exp) { + this->exportCeilings = exp; + } + + /** limit to-be-exported floors */ + void setFloors(const std::vector floors) { + this->exportFloors = floors; + } /** get all triangles grouped by obstacle */ std::vector triangulize() { std::vector res; + // get the to-be-exported floors (either "all" or "user defined") + const std::vector& floors = (exportFloors.empty()) ? (map->floors) : (exportFloors); + // process each floor - for (const Floorplan::Floor* f : map->floors) { + for (const Floorplan::Floor* f : floors) { // triangulize the floor itself (floor/ceiling) if (exportCeilings) {res.push_back(getTriangles(f));} diff --git a/wifi/estimate/ray3/ObstacleTree.h b/wifi/estimate/ray3/ObstacleTree.h deleted file mode 100644 index 8f08179..0000000 --- a/wifi/estimate/ray3/ObstacleTree.h +++ /dev/null @@ -1,139 +0,0 @@ -#ifndef OBSTACLETREE_H -#define OBSTACLETREE_H - -#include "../../../geo/Sphere3.h" -#include "Obstacle3.h" -#include - -struct ObstacleNode { - bool isLeaf = true; - Sphere3 boundSphere; - std::vector sub; - ObstacleNode(bool isLeaf = false) : isLeaf(isLeaf) {;} -}; - -struct ObstacleLeaf : public ObstacleNode { - Obstacle3D obs; - ObstacleLeaf() : ObstacleNode(true) {;} -}; - - - -class ObstacleTree { - - ObstacleNode root; - -public: - - /** append a new leaf */ - void add(ObstacleLeaf* leaf) { - root.sub.push_back(leaf); - } - - void optimize() { - while(true) { - bool did = concat(); - if (!did) {break;} - } - } - - std::vector getHits(const Ray3 ray) const { - std::vector obs; - getHits(ray, &root, obs); - return obs; - } - - void getHits(const Ray3 ray, const ObstacleNode* node, std::vector& hits) const { - for (const ObstacleNode* sub : node->sub) { - if (sub->boundSphere.intersects(ray)) { - if (sub->isLeaf) { - ObstacleLeaf* leaf = (ObstacleLeaf*)(sub); - hits.push_back(&leaf->obs); - } else { - if (sub->boundSphere.intersects(ray)) {getHits(ray, sub, hits);} - } - } - } - } - - - bool concat() { - - bool concated = false; - - // first, sort all elements by radius (smallest first) - auto compRadius = [] (const ObstacleNode* l1, const ObstacleNode* l2) { - return l1->boundSphere.radius < l2->boundSphere.radius; - }; - std::sort(root.sub.begin(), root.sub.end(), compRadius); - - ObstacleNode newRoot; - - // combine nearby elements - //for (size_t i = 0; i < root.sub.size(); ++i) { - while(true) { - - // get [and remove] the next element - ObstacleLeaf* l0 = (ObstacleLeaf*) root.sub[0]; - root.sub.erase(root.sub.begin()+0); - - // find another element that yields minimal increase in volume - auto compNear = [l0] (const ObstacleNode* l1, const ObstacleNode* l2) { - const float d1 = Sphere3::join(l0->boundSphere, l1->boundSphere).radius; - const float d2 = Sphere3::join(l0->boundSphere, l2->boundSphere).radius; - return d1 < d2; - }; - auto it = std::min_element(root.sub.begin(), root.sub.end(), compNear); - ObstacleNode* l1 = *it; - - float increment = Sphere3::join(l0->boundSphere, l1->boundSphere).radius / l0->boundSphere.radius; - - const bool combine = (root.sub.size() > 1) && (it != root.sub.end()) && (increment < 1.75); - - if (combine) { - - // combine both into a new node - ObstacleNode* node = new ObstacleNode(); - node->sub.push_back(l0); - node->sub.push_back(*it); - node->boundSphere = Sphere3::join(l0->boundSphere, (*it)->boundSphere); - root.sub.erase(it); - newRoot.sub.push_back(node); - - concated = true; - - } else { - - ObstacleNode* node = new ObstacleNode(); - node->sub.push_back(l0); - node->boundSphere = l0->boundSphere; - newRoot.sub.push_back(node); - - } - - // done? - if (root.sub.size() == 1) { - - ObstacleNode* node = new ObstacleNode(); - node->sub.push_back(root.sub.front()); - node->boundSphere = root.sub.front()->boundSphere; - newRoot.sub.push_back(node); - break; - - } else if (root.sub.size() == 0) { - break; - } - - //--i; - - } - - root = newRoot; - - return concated; - - } - -}; - -#endif // OBSTACLETREE_H diff --git a/wifi/estimate/ray3/WifiRayTrace3D.h b/wifi/estimate/ray3/WifiRayTrace3D.h index 06321d3..c94d7bd 100644 --- a/wifi/estimate/ray3/WifiRayTrace3D.h +++ b/wifi/estimate/ray3/WifiRayTrace3D.h @@ -5,6 +5,9 @@ #include "../../../geo/Line2.h" #include "../../../geo/BBox2.h" +#include "../../../geo/volume/BVH.h" +#include "../../../geo/volume/BVHDebug.h" + #include "../../../floorplan/v2/Floorplan.h" #include "../../../floorplan/v2/FloorplanHelper.h" @@ -13,7 +16,9 @@ #include "MaterialOptions.h" #include "Obstacle3.h" #include "ModelFactory.h" -#include "ObstacleTree.h" +//#include "ObstacleTree.h" + + #include @@ -133,8 +138,29 @@ struct Hit3 { +struct Obstacle3DWrapper { + static std::vector getVertices(const Obstacle3D& obs) { + std::vector pts; + for (const Triangle3& tria : obs.triangles) { + pts.push_back(tria.p1); + pts.push_back(tria.p2); + pts.push_back(tria.p3); + } + return pts; + } + static std::vector getDebugLines(const Obstacle3D& obs) { + std::vector pts; + for (const Triangle3& tria : obs.triangles) { + pts.push_back(tria.p1); pts.push_back(tria.p2); + pts.push_back(tria.p2); pts.push_back(tria.p3); + pts.push_back(tria.p3); pts.push_back(tria.p1); + } + return pts; + } + +}; @@ -147,12 +173,11 @@ private: DataMap3Signal dm; - //std::vector obstacles; - ObstacleTree tree; + BVHDebug tree; struct Limit { - static constexpr int RAYS = 1000; - static constexpr int HITS = 16; + static constexpr int RAYS = 15000; + static constexpr int HITS = 25; static constexpr float RSSI = -110; }; @@ -174,19 +199,13 @@ public: ModelFactory fac(map); std::vector obstacles = fac.triangulize(); - - // build bounding volumes for (Obstacle3D& obs : obstacles) { - ObstacleLeaf* leaf = new ObstacleLeaf(); - leaf->obs = obs; - leaf->boundSphere = getSphereAround(obs.triangles); - if (leaf->boundSphere.radius == 0) { - throw Exception("invalid item detected"); - } - tree.add(leaf); + tree.add(obs); } - tree.optimize(); + + //tree.optimize(); + //tree.show(500, false); int xxx = 0; (void) xxx; @@ -209,6 +228,7 @@ public: std::uniform_real_distribution dy(-1.0, +1.0); std::uniform_real_distribution dz(-1.0, +1.0); + #pragma omp parallel for for (int i = 0; i < Limit::RAYS; ++i) { std::cout << "ray: " << i << std::endl; @@ -229,6 +249,8 @@ public: } +//#define USE_DEBUG + private: @@ -239,7 +261,9 @@ private: // stop? if (nextHit.invalid) { +#ifdef USE_DEBUG hitStop.push_back(nextHit.pos); +#endif return; } @@ -248,7 +272,9 @@ private: // continue? if ((nextHit.stopHere) || (ray.getRSSI(nextHit.dist) < Limit::RSSI) || (ray.getDepth() > Limit::HITS)) { +#ifdef USE_DEBUG hitStop.push_back(nextHit.pos); +#endif return; } @@ -256,11 +282,15 @@ private: // apply effects if (ray.isWithin) { leave(ray, nextHit); +#ifdef USE_DEBUG hitLeave.push_back(nextHit.pos); +#endif } else { enter(ray, nextHit); reflectAt(ray, nextHit); +#ifdef USE_DEBUG hitEnter.push_back(nextHit.pos); +#endif } @@ -365,6 +395,7 @@ private: Assert::isNear(1.0f, ray.dir.length(), 0.01f, "not normalized!"); + int hits = 0; const float MAX = 999999; Hit3 nearest; nearest.dist = MAX; @@ -387,10 +418,13 @@ private: // } - std::vector obs = tree.getHits(ray); - for (const Obstacle3D* o : obs) { - hitTest(ray, *o, nearest); - } + + auto onHit = [&] (const Obstacle3D& obs) { + ++hits; + hitTest(ray, obs, nearest); + }; + + tree.getHits(ray, onHit); } @@ -399,6 +433,8 @@ private: nearest.invalid = true; } + //std::cout << hits << std::endl; + return nearest; } @@ -445,6 +481,7 @@ private: } + /* Sphere3 getSphereAround(const std::vector& trias) { std::vector pts; @@ -463,7 +500,7 @@ private: return sphere; } - + */