Marching: add dichitomic search of intersection

This commit is contained in:
Théophile Bastian 2018-02-11 16:46:55 +01:00
parent 4804309182
commit 998c77f643
3 changed files with 73 additions and 1 deletions

View file

@ -1,6 +1,11 @@
#pragma once
#include "common_structures.hpp"
class ImplicitSurface {
public:
virtual double operator() (double x, double y, double z) const;
double operator()(const Point& pt) const {
return operator()(pt.x, pt.y, pt.z);
}
};

View file

@ -1,5 +1,7 @@
#include "MarchingCubes.hpp"
#include <cassert>
const std::vector<MarchingCubes::CubeTri>
MarchingCubes::edges_of_intersection[256] = {
};
@ -29,6 +31,15 @@ Mesh MarchingCubes::operator()() {
z += resolution)
{
Intersections intersections;
/* =====
* NOTE: this code currently computes 8 times the surface value
* at each corner point of the inspected space. This is okay
* for now, because such computations are ultra light and
* actually better than storing values.
* If a time comes when this computation is heavier, because we
* are looking at more complex implicit surfaces,
* IT WILL BE NECESSARY TO ENHANCE THIS CODE!
* ==== */
for(int dx=0; dx <= 1; dx++) {
for(int dy=0; dy <= 1; dy++) {
for(int dz=0; dz <= 1; dz++) {
@ -36,13 +47,60 @@ Mesh MarchingCubes::operator()() {
double cy = y + resolution * dy;
double cz = z + resolution * dz;
intersections.set_corner(cx, cy, cz,
surface(cx, cy, cz));
surface(cx, cy, cz) > 0);
}
}
}
const std::vector<CubeTri>& cur_triangles =
edges_of_intersection[intersections.value()];
// TODO
}
}
}
return output;
}
Point MarchingCubes::CubeEdge::at(double pos,
double bx, double by, double bz, double resolution) const
{
Point p1(
bx + x[0] * resolution,
by + y[0] * resolution,
bz + z[0] * resolution);
Point p2(
bx + x[1] * resolution,
by + y[1] * resolution,
bz + z[1] * resolution);
return Point(
0.5 * (p1.x + p2.x),
0.5 * (p1.y + p2.y),
0.5 * (p1.z + p2.z));
}
Point MarchingCubes::intersect_location(const CubeEdge& edge,
double bx, double by, double bz) const
{
Point p1 = edge.at(0, bx, by, bz, resolution);
Point p2 = edge.at(1, bx, by, bz, resolution);
std::function<Point(double, double)> compute =
[&](double low_prop, double high_prop)
{
double med_prop = (low_prop + high_prop) / 2;
Point med = edge.at(med_prop, bx, by, bz, resolution);
if(high_prop - low_prop < 1e-8)
return med;
assert(surface(p1) * surface(p2) <= 0); // Can still binary search
if(surface(p1) * surface(med) <= 0)
return compute(low_prop, med_prop);
return compute(med_prop, high_prop);
};
return compute(0, 1);
}

View file

@ -85,6 +85,11 @@ class MarchingCubes {
z[1] = z1;
}
bool x[2], y[2], z[2];
/** Get the space point at a given position of [0,1] along the edge
* when the base of the cube (ie. (0, 0, 0)) is given. */
Point at(double pos,
double bx, double by, double bz, double resolution) const;
};
struct CubeTri {
CubeTri() {}
@ -100,6 +105,10 @@ class MarchingCubes {
CubeEdge edge[3];
};
private: //methods
Point intersect_location(const CubeEdge& edge,
double bx, double by, double bz) const;
private:
static const std::vector<CubeTri> edges_of_intersection[256];