FHWS/Indoor
Archived
4
0
Fork 0
Code Issues 19 Pull Requests Releases 1 Wiki Activity

modified lib GPC for header only

Browse Source 
worked on 3d traytracing
This commit is contained in:
k-a-z-u
2017-09-06 17:04:19 +02:00
parent 845d89774d
commit c19d18a3a6
20 changed files with 884 additions and 299 deletions
Download Patch File Download Diff File Expand all files Collapse all files

8
CMakeLists.txt
View File

@@ -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")
@@ -69,6 +76,7 @@ ADD_DEFINITIONS(
-g3
# -O0
-O2
-march=native
-DWITH_TESTS

5
geo/BBox3.h
View File

@@ -24,6 +24,11 @@ public:
/** 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) {

4
geo/Point3.h
View File

@@ -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,

178
geo/Sphere3.h
View File

@@ -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<Point3>& 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;
@@ -89,38 +70,117 @@ struct Sphere3 {
// 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<Point3>& 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<Point3>& 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;
}
};

0
geo/TestSphere3.cpp
View File

2
geo/Triangle3.h
View File

@@ -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;

224
geo/volume/BVH.h
View File

@@ -1,17 +1,235 @@
#ifndef BOUNDINGVOLUMEHIERARCHY_H
#define BOUNDINGVOLUMEHIERARCHY_H
#include <vector>
#include <functional>
#include "BoundingVolume.h"
#include "BoundingVolumeAABB.h"
#include "BoundingVolumeSphere.h"
class BVH {
template <typename Element, typename Volume, typename Wrapper> class BVH {
protected:
/** one node within the tree */
struct BVHNode {
bool isLeaf = true;
Volume boundingVolume;
std::vector<BVHNode*> 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<void(const Element&)> 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<void(const Element&)> 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<Point3> verts = Wrapper::getVertices(element);
return Volume::fromVertices(verts);
}

106
geo/volume/BVHDebug.h Normal file
View File

@@ -0,0 +1,106 @@
#ifndef BVHDEBUG_H
#define BVHDEBUG_H
#include "BVH.h"
#include <KLib/misc/gnuplot/GnuplotSplot.h>
#include <KLib/misc/gnuplot/GnuplotSplotElementPoints.h>
#include <KLib/misc/gnuplot/GnuplotSplotElementColorPoints.h>
#include <KLib/misc/gnuplot/GnuplotSplotElementLines.h>
#include <KLib/misc/gnuplot/Gnuplot.h>
#include "../BBox3.h"
#include <random>
/** adds some debug helpers to the BVH */
template <typename Element, typename Volume, typename Wrapper> class BVHDebug : public BVH<Element, Volume, Wrapper> {
using BVHNode = typename BVH<Element, Volume, Wrapper>::BVHNode;
using BVHLeaf = typename BVH<Element, Volume, Wrapper>::BVHLeaf;
public:
// std::vecto<std::string> 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<Point3> 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<float> 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<Point3> 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

18
geo/volume/BoundingVolume.h
View File

@@ -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;
};

2
geo/volume/BoundingVolumeAABB.h
View File

@@ -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);
}

48
geo/volume/BoundingVolumeSphere.h
View File

@@ -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<Point3>& verts) {
BoundingVolumeSphere bvs;
bvs.adjustToPointSet(verts);
return bvs;
}
static BoundingVolumeSphere join(const BoundingVolumeSphere a, const BoundingVolumeSphere b) {
return BoundingVolumeSphere(Sphere3::join(a, b));
}
};

0
lib/gpc/gpc.cpp → lib/gpc/gpc.cpp.h
View File

19
lib/gpc/gpc.h
View File

@@ -95,34 +95,36 @@ typedef struct /* Tristrip set structure */
===========================================================================
*/
void gpc_read_polygon (FILE *infile_ptr,
inline void gpc_read_polygon (FILE *infile_ptr,
int read_hole_flags,
gpc_polygon *polygon);
void gpc_write_polygon (FILE *outfile_ptr,
inline void gpc_write_polygon (FILE *outfile_ptr,
int write_hole_flags,
gpc_polygon *polygon);
void gpc_add_contour (gpc_polygon *polygon,
inline void gpc_add_contour (gpc_polygon *polygon,
gpc_vertex_list *contour,
int hole);
void gpc_polygon_clip (gpc_op set_operation,
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,
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,
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
===========================================================================
*/

4
main.cpp
View File

@@ -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*";

198
tests/geo/TestBVH.cpp Normal file
View File

@@ -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<Point3> getVertices(const BBox3& bbox) {
return {bbox.getMin(), bbox.getMax()};
}
static std::vector<Point3> 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<Point3> 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<Point3> getVertices(const Obstacle3D& o) {
std::vector<Point3> 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<Point3> getDebugLines(const Obstacle3D& o) {
std::vector<Point3> 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<BBox3, BoundingVolumeSphere, Wrapper> 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<BBox3, BoundingVolumeSphere, Wrapper> 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<BBox3, BoundingVolumeSphere, Wrapper> 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<Obstacle3D> obs = fac.triangulize();
BVHDebug<Obstacle3D, BoundingVolumeSphere, WrapperObs3D> 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<BBox3, BoundingVolumeSphere, Wrapper> tree;
std::minstd_rand gen;
std::uniform_real_distribution<float> dPos(-4.0, +4.0);
std::uniform_real_distribution<float> 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

49
tests/ray/TestRayTrace3.cpp
View File

@@ -12,23 +12,32 @@ 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<std::chrono::high_resolution_clock> start = std::chrono::high_resolution_clock::now();
const DataMap3Signal& dms = rt.estimate();
std::chrono::time_point<std::chrono::high_resolution_clock> end = std::chrono::high_resolution_clock::now();
auto result = std::chrono::duration_cast<std::chrono::milliseconds>(end-start).count();
std::cout << "it took: " << result << " msec" << std::endl;
const char sep = ';';
if (1 == 1) {
const char sep = ' ';
int lines = 0;
std::stringstream tmp;
std::ofstream out("/mnt/vm/rays.xyz.txt");
auto lambda = [&] (const float x, const float y, const float z, const DataMap3SignalEntry& e) {
const float min = -120;
@@ -36,9 +45,12 @@ TEST(RayTrace3, test) {
float rssi = e.getMaxRSSI();
if (rssi > max) {rssi = max;}
if (z < 1.0 || z > 1.0) {return;}
if (rssi > -100) {
const float v = ((rssi - min) / (max-min)) * 255;
out
++lines;
const int v = ((rssi - min) / (max-min)) * 255; // color
tmp
<< x << sep << y << sep << z << sep
<< v << sep << v << sep << v
<< "\n";
@@ -47,9 +59,28 @@ TEST(RayTrace3, test) {
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::ofstream outHits("/mnt/vm/hits.xyz.txt");
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";
}
@@ -63,4 +94,6 @@ TEST(RayTrace3, test) {
}
}
#endif

6
wifi/estimate/ray3/DataMap3.h
View File

@@ -5,6 +5,7 @@
#include "../../../geo/BBox3.h"
#include <vector>
#include <functional>
#include <mutex>
template <typename T> class DataMap3 {
@@ -189,6 +190,8 @@ private:
struct DataMap3SignalEntry {
struct Entry {
float rssi;
float distanceToAP;
@@ -198,8 +201,11 @@ struct DataMap3SignalEntry {
std::vector<Entry> 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 {

14
wifi/estimate/ray3/ModelFactory.h
View File

@@ -17,6 +17,7 @@ private:
bool exportCeilings = true;
bool exportObstacles = true;
bool exportWallTops = false;
std::vector<Floorplan::Floor*> 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<Floorplan::Floor*> floors) {
this->exportFloors = floors;
}
/** get all triangles grouped by obstacle */
std::vector<Obstacle3D> triangulize() {
std::vector<Obstacle3D> res;
// get the to-be-exported floors (either "all" or "user defined")
const std::vector<Floorplan::Floor*>& 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));}

139
wifi/estimate/ray3/ObstacleTree.h
View File

@@ -1,139 +0,0 @@
#ifndef OBSTACLETREE_H
#define OBSTACLETREE_H
#include "../../../geo/Sphere3.h"
#include "Obstacle3.h"
#include <algorithm>
struct ObstacleNode {
bool isLeaf = true;
Sphere3 boundSphere;
std::vector<ObstacleNode*> 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<Obstacle3D*> getHits(const Ray3 ray) const {
std::vector<Obstacle3D*> obs;
getHits(ray, &root, obs);
return obs;
}
void getHits(const Ray3 ray, const ObstacleNode* node, std::vector<Obstacle3D*>& 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

77
wifi/estimate/ray3/WifiRayTrace3D.h
View File

@@ -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 <random>
@@ -133,8 +138,29 @@ struct Hit3 {
struct Obstacle3DWrapper {
static std::vector<Point3> getVertices(const Obstacle3D& obs) {
std::vector<Point3> 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<Point3> getDebugLines(const Obstacle3D& obs) {
std::vector<Point3> 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<Obstacle3D> obstacles;
ObstacleTree tree;
BVHDebug<Obstacle3D, BoundingVolumeSphere, Obstacle3DWrapper> 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<Obstacle3D> 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(obs);
}
tree.add(leaf);
}
tree.optimize();
//tree.optimize();
//tree.show(500, false);
int xxx = 0; (void) xxx;
@@ -209,6 +228,7 @@ public:
std::uniform_real_distribution<float> dy(-1.0, +1.0);
std::uniform_real_distribution<float> 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<Obstacle3D*> 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<Triangle3>& trias) {
std::vector<Point3> pts;
@@ -463,7 +500,7 @@ private:
return sphere;
}
*/