#include <dune/curvedsurfacegrid/surfacedistance/surfacedistance.hh>
#include <iostream>
#include <cmath>
using namespace Dune::Curved;
// #########################
// ## Geometry-Routines ##
// #########################
//--------------------< coordinate_extrema_triangle >
//returns in 'small_coord' the smallest coordinate of 'sm' in the given dimension and in
//'big_coord' the biggest coordinate in the given dimension. Note: 'sm' has to be a triangle
//in a 2d or 3d world.
void coordinate_extrema_triangle(Mesh &m, Element &sm, int dimension, double &small_coord,
double &big_coord){
if(m.vertices[sm.v[0]][dimension] <= m.vertices[sm.v[1]][dimension]){
if(m.vertices[sm.v[0]][dimension] <= m.vertices[sm.v[2]][dimension]){
small_coord = m.vertices[sm.v[0]][dimension];
if(m.vertices[sm.v[1]][dimension] >= m.vertices[sm.v[2]][dimension])
big_coord = m.vertices[sm.v[1]][dimension];
else
big_coord = m.vertices[sm.v[2]][dimension];
}else{
small_coord = m.vertices[sm.v[2]][dimension];
big_coord = m.vertices[sm.v[1]][dimension];
}
}else{
if(m.vertices[sm.v[1]][dimension] <= m.vertices[sm.v[2]][dimension]){
small_coord = m.vertices[sm.v[1]][dimension];
if(m.vertices[sm.v[0]][dimension] >= m.vertices[sm.v[2]][dimension])
big_coord = m.vertices[sm.v[0]][dimension];
else
big_coord = m.vertices[sm.v[2]][dimension];
}else{
small_coord = m.vertices[sm.v[2]][dimension];
big_coord = m.vertices[sm.v[0]][dimension];
}
}
}
//--------------------< distance_point_line_2d >
//calculates the distance between a given point and the line segment between 'lin0' and 'lin1'.
double distance_point_line_2d(const Vertex &point, const Vertex &lin0, const Vertex &lin1,
Vertex &hit_vertex){
//the (complete) line through 'lin0' and 'lin1' is given by 'l = lin0 + t*g' with g as follows:
double gx = lin1[0] - lin0[0];
double gy = lin1[1] - lin0[1];
//the normal-vector from 'point' to l is 'v = lin0 + t*g - point' for some t which is determined
// by 'v*g = 0'
double t = (-gx * (lin0[0] - point[0]) - gy * (lin0[1] - point[1])) / (gx * gx + gy * gy);
//if the orthogonal projection of 'point' onto l is not on the given line segment then one of the
//vertices 'lin0' or 'lin1' is the nearest point to 'point'
if(t < 0){
//'lin0' is nearest
hit_vertex[0] = lin0[0];
hit_vertex[1] = lin0[1];
}else if(t > 1){
//'lin1' is nearest
hit_vertex[0] = lin1[0];
hit_vertex[1] = lin1[1];
}else{
//orthogonal projection is nearest
hit_vertex[0] = lin0[0] + t * gx;
hit_vertex[1] = lin0[1] + t * gy;
}
//calculate distance
double dx = point[0] - hit_vertex[0];
double dy = point[1] - hit_vertex[1];
return sqrt(dx * dx + dy * dy);
}
//--------------------< point_in_triangle_generous_2d >
bool point_in_triangle_generous_2d(const Vertex &p, const Vertex &a, const Vertex &b,
const Vertex &c){
int i;
Vertex v0, v1, v2;
double dot00, dot01, dot02, dot11, dot12, invDenom, u, v;
//compute direction vectors
for(i = 0; i < 2; ++i) v0[i] = c[i] - a[i];
for(i = 0; i < 2; ++i) v1[i] = b[i] - a[i];
for(i = 0; i < 2; ++i) v2[i] = p[i] - a[i];
//compute dot products
dot00 = v0[0] * v0[0] + v0[1] * v0[1]; //dot(v0, v0)
dot01 = v0[0] * v1[0] + v0[1] * v1[1]; //dot(v0, v1)
dot02 = v0[0] * v2[0] + v0[1] * v2[1]; //dot(v0, v2)
dot11 = v1[0] * v1[0] + v1[1] * v1[1]; //dot(v1, v1)
dot12 = v1[0] * v2[0] + v1[1] * v2[1]; //dot(v1, v2)
//compute barycentric coordinates
invDenom = 1.0 / (dot00 * dot11 - dot01 * dot01);
u = (dot11 * dot02 - dot01 * dot12) * invDenom;
v = (dot00 * dot12 - dot01 * dot02) * invDenom;
//check if point is in triangle
return (u >= -epsilon) && (v >= -epsilon) && (u + v - 1.0 <= epsilon);
}
//--------------------< distance_point_line_3d >
//calculates the distance between a given point and the line segment between 'lin0' and 'lin1'.
double distance_point_line_3d(const Vertex &point, const Vertex &lin0, const Vertex &lin1,
Vertex &hit_vertex){
//the (complete) line through 'lin0' and 'lin1' is given by 'l = lin0 + t*g' with g as follows:
double gx = lin1[0] - lin0[0];
double gy = lin1[1] - lin0[1];
double gz = lin1[2] - lin0[2];
//the normal-vector from 'point' to l is 'v = lin0 + t*g - point' for some t which is determined
// by 'v*g = 0'
double t = (-gx * (lin0[0] - point[0]) - gy * (lin0[1] - point[1]) - gz * (lin0[2] - point[2]))
/ (gx * gx + gy * gy + gz * gz);
//if the orthogonal projection of 'point' onto l is not on the given line segment then one of the
//vertices 'lin0' or 'lin1' is the nearest point to 'point'
if(t < 0){
//'lin0' is nearest
hit_vertex[0] = lin0[0];
hit_vertex[1] = lin0[1];
hit_vertex[2] = lin0[2];
}else if(t > 1){
//'lin1' is nearest
hit_vertex[0] = lin1[0];
hit_vertex[1] = lin1[1];
hit_vertex[2] = lin1[2];
}else{
//orthogonal projection is nearest
hit_vertex[0] = lin0[0] + t * gx;
hit_vertex[1] = lin0[1] + t * gy;
hit_vertex[2] = lin0[2] + t * gz;
}
//calculate distance
double dx = point[0] - hit_vertex[0];
double dy = point[1] - hit_vertex[1];
double dz = point[2] - hit_vertex[2];
return sqrt(dx * dx + dy * dy + dz * dz);
}
//--------------------< intersection_line_line_with_intersection_3d >
bool intersection_line_line_with_intersection_3d(const Vertex &p1, const Vertex &p2,
const Vertex &p3, const Vertex &p4, Vertex &intersection){
Vertex p13, p43, p21;
double d1343, d4321, d1321, d4343, d2121;
double numer, denom;
double mua, mub;
int ci = -1;
p13[0] = p1[0] - p3[0];
p13[1] = p1[1] - p3[1];
p13[2] = p1[2] - p3[2];
p43[0] = p4[0] - p3[0];
p43[1] = p4[1] - p3[1];
p43[2] = p4[2] - p3[2];
p21[0] = p2[0] - p1[0];
p21[1] = p2[1] - p1[1];
p21[2] = p2[2] - p1[2];
//test for collinearity
if(fabs(p21[1] * p43[2] - p21[2] * p43[1]) < epsilon
&& fabs(p21[2] * p43[0] - p21[0] * p43[2]) < epsilon
&& fabs(p21[0] * p43[1] - p21[1] * p43[0]) < epsilon){
//lines are collinear; now test for an intersection
double x11, x12, x21, x22, y11, y12, y21, y22, z11, z12, z21, z22;
if(p1[0] <= p2[0]){
x11 = p1[0];
x12 = p2[0];
}else{
x11 = p2[0];
x12 = p1[0];
}
if(p1[1] <= p2[1]){
y11 = p1[1];
y12 = p2[1];
}else{
y11 = p2[1];
y12 = p1[1];
}
if(p1[2] <= p2[2]){
z11 = p1[2];
z12 = p2[2];
}else{
z11 = p2[2];
z12 = p1[2];
}
if(p3[0] <= p4[0]){
x21 = p3[0];
x22 = p4[0];
}else{
x21 = p4[0];
x22 = p3[0];
}
if(p3[1] <= p4[1]){
y21 = p3[1];
y22 = p4[1];
}else{
y21 = p4[1];
y22 = p3[1];
}
if(p3[2] <= p4[2]){
z21 = p3[2];
z22 = p4[2];
}else{
z21 = p4[2];
z22 = p3[2];
}
if(x11 <= x22 && x12 >= x21 && y11 <= y22 && y12 >= y21 && z11 <= z22 && z12 >= z21){
if(fabs(p21[0]) > epsilon){
mua = -p13[0] / p21[0];
if(fabs(p1[1] + p21[1] * mua - p3[1]) < epsilon
&& fabs(p1[2] + p21[2] * mua - p3[2]) < epsilon) ci = 0;
}else if(fabs(p21[1]) > epsilon){
mua = -p13[1] / p21[1];
if(fabs(p1[0] + p21[0] * mua - p3[0]) < epsilon
&& fabs(p1[2] + p21[2] * mua - p3[2]) < epsilon) ci = 1;
}else{
mua = -p13[2] / p21[2];
if(fabs(p1[0] + p21[0] * mua - p3[0]) < epsilon
&& fabs(p1[1] + p21[1] * mua - p3[1]) < epsilon) ci = 2;
}
if(ci != -1){
//determine best intersection-point
//priority is p3 > p4 > p3 > p2 (one of these three points must be part of both line
//segments. p3 is only priorized over p4 if it is not a vertex of the first line segment
if(mua > epsilon && mua - 1 < -epsilon){
intersection[0] = p3[0];
intersection[1] = p3[1];
intersection[2] = p3[2];
}else{
mub = (p4[ci] - p1[ci]) / p21[ci];
if(mub > -epsilon && mub - 1 < epsilon){
intersection[0] = p4[0];
intersection[1] = p4[1];
intersection[2] = p4[2];
}else if(mua > -epsilon && mua-1 < epsilon){
intersection[0] = p3[0];
intersection[1] = p3[1];
intersection[2] = p3[2];
}else{
intersection[0] = p2[0];
intersection[1] = p2[1];
intersection[2] = p2[2];
}
}
return true;
}
}
return false;
}
d1343 = p13[0] * p43[0] + p13[1] * p43[1] + p13[2] * p43[2];
d4321 = p43[0] * p21[0] + p43[1] * p21[1] + p43[2] * p21[2];
d1321 = p13[0] * p21[0] + p13[1] * p21[1] + p13[2] * p21[2];
d4343 = p43[0] * p43[0] + p43[1] * p43[1] + p43[2] * p43[2];
d2121 = p21[0] * p21[0] + p21[1] * p21[1] + p21[2] * p21[2];
denom = d2121 * d4343 - d4321 * d4321;
if(fabs(denom) < epsilon) return false;
numer = d1343 * d4321 - d1321 * d4343;
mua = numer / denom;
mub = (d1343 + d4321 * mua) / d4343;
if(mua < -epsilon || mua - 1 > epsilon || mub < -epsilon || mub - 1 > epsilon) return false;
//determine potential intersection-point
intersection[0] = p1[0] + mua * p21[0];
intersection[1] = p1[1] + mua * p21[1];
intersection[2] = p1[2] + mua * p21[2];
//check if this really is an intersection
if(fabs(intersection[0] - p3[0] - mub * p43[0]) < epsilon
&& fabs(intersection[1] - p3[1] - mub * p43[1]) < epsilon
&& fabs(intersection[2] - p3[2] - mub * p43[2]) < epsilon){
return true;
}else{
return false;
}
}
//--------------------< intersection_line_triangle_3d >
//This is a modified version of Softsurfers line-triangle-intersection-routine. The original
//routine comes with the following copyright notice:
// Copyright 2001, softSurfer (www.softsurfer.com)
// This code may be freely used and modified for any purpose
// providing that this copyright notice is included with it.
// SoftSurfer makes no warranty for this code, and cannot be held
// liable for any real or imagined damage resulting from its use.
// Users of this code must verify correctness for their application.
//calculates the intersection point of a line segment and a triangle if it exists
// input: a line segment (lin0, lin1), and a triangle (tri0, tri1, tri2)
// output: sol = intersection point (if it exists)
// return: 0 disjoint
// 1 intersection inside triangle
// 2,3,4 intersection in edge of triangle (edge_index is 0, 1, 2 respectively)
int intersection_line_triangle_3d(const Vertex &tri0, const Vertex &tri1,
const Vertex &tri2, const Vertex &lin0, const Vertex &lin1, Vertex &sol){
Vertex hit_vertex;
Vertex u, v, n; //triangle vectors
Vertex dir, w0, w; //ray vectors
double r, a, b; //params to calc ray-plane intersect
int i, e;
double epsilon2 = 100 * epsilon;
double uu, uv, vv, wu, wv, D;
double s, t;
Vertex distances;
//get triangle edge vectors and plane normal
for(i = 0; i < 3; i++) u[i] = tri1[i] - tri0[i];
for(i = 0; i < 3; i++) v[i] = tri2[i] - tri0[i];
n[0] = u[1] * v[2] - u[2] * v[1];
n[1] = u[2] * v[0] - u[0] * v[2];
n[2] = u[0] * v[1] - u[1] *v[0];
if(fabs(n[0]) < epsilon && fabs(n[1]) < epsilon && fabs(n[2]) < epsilon){
if(intersection_line_line_with_intersection_3d(tri0, tri1, lin0, lin1, sol)) return 4;
else if(intersection_line_line_with_intersection_3d(tri1, tri2, lin0, lin1, sol)) return 2;
else if(intersection_line_line_with_intersection_3d(tri2, tri0, lin0, lin1, sol)) return 3;
}
for(i = 0; i < 3; i++) dir[i] = lin1[i] - lin0[i]; //line direction vector
for(i = 0; i < 3; i++) w0[i] = lin0[i] - tri0[i];
a = -(n[0] * w0[0] + n[1] * w0[1] + n[2] * w0[2]);
b = n[0] * dir[0] + n[1] * dir[1] + n[2] * dir[2];
if(fabs(b) < epsilon){ //ray is parallel to triangle plane
if(fabs(a) < epsilon){
//test whether one of the line's vertices lies inside the triangle
for(e = 0; e < 2; e++){
if(e == 0) for(i = 0; i < 3; i++) sol[i] = lin0[i];
else for(i = 0; i < 3; i++) sol[i] = lin1[i];
uu = u[0] * u[0] + u[1] * u[1] + u[2] * u[2];
uv = u[0] * v[0] + u[1] * v[1] + u[2] * v[2];
vv = v[0] * v[0] + v[1] * v[1] + v[2] * v[2];
for(i = 0; i < 3; i++) w[i] = sol[i] - tri0[i];
wu = w[0] * u[0] + w[1] * u[1] + w[2] * u[2];
wv = w[0] * v[0] + w[1] * v[1] + w[2] * v[2];
D = uv * uv - uu * vv;
s = (uv * wv - vv * wu) / D;
if (s < -epsilon || s - 1.0 > epsilon) continue; //sol is outside triangle
t = (uv * wu - uu * wv) / D;
if (t < -epsilon || (s + t - 1.0) > epsilon) continue; //sol is outside triangle
distances[0] = distance_point_line_3d(sol, tri0, tri1, hit_vertex);
distances[1] = distance_point_line_3d(sol, tri1, tri2, hit_vertex);
distances[2] = distance_point_line_3d(sol, tri2, tri0, hit_vertex);
if(distances[0] > epsilon2 && distances[1] > epsilon2 && distances[2] > epsilon2){
return 1;
}
}
//test whether ray intersects one of the triangle's edges
if(intersection_line_line_with_intersection_3d(tri0, tri1, lin0, lin1, sol)) return 4;
else if(intersection_line_line_with_intersection_3d(tri1, tri2, lin0, lin1, sol)) return 2;
else if(intersection_line_line_with_intersection_3d(tri2, tri0, lin0, lin1, sol)) return 3;
}
return 0; //ray disjoint from plane
}
//get intersection point of ray with triangle plane
r = a / b;
if(r < -epsilon || r - 1.0 > epsilon){
return 0; //ray goes away from triangle => no intersection
}
for(i = 0; i < 3; i++) sol[i] = lin0[i] + r * dir[i];
//is sol inside triangle?
uu = u[0] * u[0] + u[1] * u[1] + u[2] * u[2];
uv = u[0] * v[0] + u[1] * v[1] + u[2] * v[2];
vv = v[0] * v[0] + v[1] * v[1] + v[2] * v[2];
for(i = 0; i < 3; i++) w[i] = sol[i] - tri0[i];
wu = w[0] * u[0] + w[1] * u[1] + w[2] * u[2];
wv = w[0] * v[0] + w[1] * v[1] + w[2] * v[2];
D = uv * uv - uu * vv;
//get and test parametric coords
s = (uv * wv - vv * wu) / D;
if (s < -epsilon || s - 1.0 > epsilon) return 0; //sol is outside triangle
t = (uv * wu - uu * wv) / D;
if (t < -epsilon || (s + t - 1.0) > epsilon) return 0; //sol is outside triangle
//calculate distance of sol with each edge of the triangle
distances[0] = distance_point_line_3d(sol, tri0, tri1, hit_vertex);
if(distances[0] < epsilon2) return 4;
distances[1] = distance_point_line_3d(sol, tri1, tri2, hit_vertex);
if(distances[1] < epsilon2) return 2;
distances[2] = distance_point_line_3d(sol, tri2, tri0, hit_vertex);
if(distances[2] < epsilon2) return 3;
return 1; //sol is inside triangle
}
//--------------------< distance_point_triangle_3d >
//calculates the distance between a given 'point' and the triangle given by 'tri0', 'tri1', 'tri2'.
double distance_point_triangle_3d(const Vertex &point,
const Vertex &tri0, const Vertex &tri1, const Vertex &tri2, Vertex &hit_vertex){
//triangle lies in plane 'p = tri0 + u*g + v*h', where g and h are defined as follows:
double gx = tri1[0] - tri0[0];
double gy = tri1[1] - tri0[1];
double gz = tri1[2] - tri0[2];
double hx = tri2[0] - tri0[0];
double hy = tri2[1] - tri0[1];
double hz = tri2[2] - tri0[2];
//a vector from 'point' to p can be described by 'l = tri0 + u*g + v*h - point' and since we
//search the distance we need such a vector that is orthogonal to p. Thus 'l*g = 0' and 'l*h = 0'.
//We get a system of two equations of the form:
//a1*u + b1*v = c1
//b1*u + b2*v = c2
//with a1, b1, b2, c1, c2 as follows:
double b1 = gx * hx + gy * hy + gz * hz;
double c2 = -hx * (tri0[0] - point[0]) - hy * (tri0[1] - point[1]) - hz * (tri0[2] - point[2]);
double b2 = hx * hx + hy * hy + hz * hz;
double c1 = -gx * (tri0[0] - point[0]) - gy * (tri0[1] - point[1]) - gz * (tri0[2] - point[2]);
double a1 = gx * gx + gy * gy + gz * gz;
//solve the system:
double u, v;
if(fabs(b1) < epsilon){
v = c2 / b2;
u = c1 / a1;
}else{
v = (c1 * b1 - c2 * a1) / (b1 * b1 - b2 * a1);
u = (c1 - v * b1) / a1;
}
//now test whether the orthogonal projection of 'point' onto p lies inside the triangle
double distance;
if(u >= 0 && v >= 0 && u + v <= 1){
//yes, inside. The length of the normal-vector is the desired distance
double fx = hit_vertex[0] = tri0[0] + u * gx + v * hx;
double fy = hit_vertex[1] = tri0[1] + u * gy + v * hy;
double fz = hit_vertex[2] = tri0[2] + u * gz + v * hz;
fx -= point[0];
fy -= point[1];
fz -= point[2];
distance = sqrt(fx * fx + fy * fy + fz * fz);
}else{
//no, not inside. The desired distance is the distance to one of the triangle-edges
distance = distance_point_line_3d(point, tri0, tri1, hit_vertex);
Vertex current_hit_vertex;
double current_distance = distance_point_line_3d(point, tri1, tri2, current_hit_vertex);
if(current_distance < distance){
distance = current_distance;
hit_vertex = current_hit_vertex;
}
current_distance = distance_point_line_3d(point, tri2, tri0, current_hit_vertex);
if(current_distance < distance){
distance = current_distance;
hit_vertex = current_hit_vertex;
}
}
return distance;
}
// ##########################
// ## BackGrid - General ##
// ##########################
//--------------------< static members for configuration >
int BackGrid::segmentation_degree = -1; //-1 means automatic switching dependent on mesh-size
//--------------------< construcor >
BackGrid::BackGrid(Mesh &mesh, double expansion_width, bool no_scaling):mesh(mesh){
int x, y;
Vertex dist;
int seg_deg;
current_operation_index = 0;
//automatically determine segmentation degree if desired
if(segmentation_degree <= 0){
if(mesh.elements.size() < 20000) seg_deg = 64;
else if(mesh.elements.size() < 800000) seg_deg = 128;
else seg_deg = 256;
}else{
seg_deg = segmentation_degree;
}
//calculate minimal and maximal coordinates of backgrid
mesh.determine_coordinate_extrema(minc, maxc);
if(expansion_width < 0.0) expansion_width = coarse_epsilon;
for(y = 0; y < mesh.dim_of_world; ++y){
minc[y] -= expansion_width;
maxc[y] += expansion_width;
}
//calculate sizes of segments
if(no_scaling){
//boxes must be cubes (all edges of same length)
double maxdist = 0.0;
y = 0;
for(x = 0; x < mesh.dim_of_world; ++x){
dist[x] = maxc[x] - minc[x];
if(dist[x] > maxdist){
maxdist = dist[x];
y = x;
}
}
for(x = 0; x < mesh.dim_of_world; ++x){
seg_size[x] = maxdist / (double)(seg_deg);
segments[x] = int((dist[x] / maxdist) * seg_deg);
if(dist[x] < dist[y]) ++segments[x];
}
segments[y] = seg_deg;
min_seg_size = seg_size[0];
}else{
//boxes need not be cubes
double distpower = 1.0;
min_seg_size = 9e100;
for(x = 0; x < mesh.dim_of_world; ++x){
dist[x] = maxc[x] - minc[x];
if(dist[x] > epsilon && dist[x] < min_seg_size) min_seg_size = dist[x];
}
for(x = 0; x < mesh.dim_of_world; ++x){
if(dist[x] <= epsilon) dist[x] = min_seg_size;
else dist[x] = 1.0 / dist[x];
distpower *= dist[x];
}
distpower = pow(distpower * pow(double(seg_deg), mesh.dim_of_world), 1.0 / mesh.dim_of_world);
double tmp_segments[3];
int latest_zero_dimension = 0, latest_non_zero_dimension = 0, num_zero_dimensions = 0;
for(x = 0; x < mesh.dim_of_world; ++x){
tmp_segments[x] = int(distpower / dist[x]);
if(tmp_segments[x] < 1.0){
latest_zero_dimension = x;
++num_zero_dimensions;
}else{
latest_non_zero_dimension = x;
}
}
//corrections for zero-dimensions:
if(mesh.dim_of_world == 2){
if(num_zero_dimensions > 0){
double factor = pow(seg_deg, 2) / tmp_segments[latest_non_zero_dimension] - 1;
tmp_segments[latest_zero_dimension] = 1.0;
tmp_segments[latest_non_zero_dimension] = (factor + 1.0)
* tmp_segments[latest_non_zero_dimension];
}
}else if(mesh.dim_of_world == 3){
if(num_zero_dimensions == 1){
int other_dimension = (latest_zero_dimension + 1) % 2;
if(other_dimension == latest_non_zero_dimension) other_dimension =
(other_dimension + 1) % 2;
double factor = sqrt(pow(seg_deg, 3) / (tmp_segments[latest_non_zero_dimension]
* tmp_segments[other_dimension])) - 1;
tmp_segments[latest_zero_dimension] = 1.0;
tmp_segments[latest_non_zero_dimension] = (factor + 1.0)
* tmp_segments[latest_non_zero_dimension];
tmp_segments[other_dimension] = (factor + 1.0) * tmp_segments[other_dimension];
}else if(num_zero_dimensions == 2){
int other_dimension = (latest_zero_dimension + 1) % 2;
if(other_dimension == latest_non_zero_dimension) other_dimension =
(other_dimension + 1) % 2;
double factor = pow(seg_deg, 3) / tmp_segments[latest_non_zero_dimension] - 1;
tmp_segments[latest_zero_dimension] = 1.0;
tmp_segments[other_dimension] = 1.0;
tmp_segments[latest_non_zero_dimension] = (factor + 1.0)
* tmp_segments[latest_non_zero_dimension];
}
}
for(x = 0; x < mesh.dim_of_world; ++x){
segments[x] = tmp_segments[x];
if(segments[x] <= 0) segments[x] = 1;
}
min_seg_size = 9e100;
for(x = 0; x < mesh.dim_of_world; ++x){
seg_size[x] = (maxc[x] - minc[x]) / (double)segments[x];
if(min_seg_size > seg_size[x]) min_seg_size = seg_size[x];
}
}
//create data structure for backgrid
if(mesh.dim_of_world == 2) segments[2] = 1;
grid.resize(segments[0]);
for(x = 0; x < segments[0]; ++x){
grid[x].resize(segments[1]);
for(y = 0; y < segments[1]; ++y){
grid[x][y].resize(segments[2]);
}
}
//attach mesh
if(mesh.dim_of_world == 2){
attach_mesh_2d();
}else if(mesh.dim_of_world == 3){
attach_mesh_3d();
}
}
//--------------------< debug_output_dimensions >
void BackGrid::debug_output_dimensions(){
int x, y, z, sum = 0, cnt = 0;
for(x = 0; x < segments[0]; ++x){
for(y = 0; y < segments[1]; ++y){
for(z = 0; z < segments[2]; ++z){
sum += grid[x][y][z].size();
++cnt;
}
}
}
std::cout << "\nbackgrid:" << std::endl;
std::cout << "dimensions: " << minc << ", " << maxc << std::endl;
std::cout << "segment-count: "
<< segments[0] << " x " << segments[1] << " x " << segments[2] << std::endl;
std::cout << "segment-size: "
<< seg_size[0] << " x " << seg_size[1] << " x " << seg_size[2] << std::endl;
std::cout << "average entries: " << double(sum) / cnt << std::endl;
std::cout << "total entries: " << sum << "\n" << std::endl;
}
// ############################
// ## BackGrid - 2D Worlds ##
// ############################
//--------------------< attach_mesh_2d >
void BackGrid::attach_mesh_2d(){
int x, y, cnt = 0;
double cx1, cx2, cy1, cy2, vx1, vx2, vy1, vy2, s;
for(std::list<Element>::iterator it = mesh.elements.begin(),
it_end = mesh.elements.end(); it != it_end; ++it){
it->operation_index = 0;
//determine leftmost and rightmost vertex of 'it'
if(mesh.vertices[it->v[0]][0] < mesh.vertices[it->v[1]][0]){
vx1 = mesh.vertices[it->v[0]][0];
vx2 = mesh.vertices[it->v[1]][0];
vy1 = mesh.vertices[it->v[0]][1];
vy2 = mesh.vertices[it->v[1]][1];
}else{
vx1 = mesh.vertices[it->v[1]][0];
vx2 = mesh.vertices[it->v[0]][0];
vy1 = mesh.vertices[it->v[1]][1];
vy2 = mesh.vertices[it->v[0]][1];
}
//determine for every segment whether the current element intersects it
for(y = 0; y < segments[1]; ++y){
for(x = 0; x < segments[0]; ++x){
if(vx1 - minc[0] - (x + 1) * seg_size[0] <= epsilon
&& vx2 - minc[0] - x * seg_size[0] >= -epsilon){
if(fabs(vx2 - vx1) < epsilon){
//element is a vertical line-segment
cy1 = vy1;
cy2 = vy2;
}else{
//element is not vertical; determine its gradient
s = (vy2 - vy1) / (vx2 - vx1);
cx1 = ((vx1 > minc[0] + x * seg_size[0])? vx1 : minc[0] + x * seg_size[0]);
cx2 = ((vx2 < minc[0] + (x + 1) * seg_size[0])? vx2 : minc[0] + (x + 1) * seg_size[0]);
cy1 = (cx1 - vx1) * s + vy1;
cy2 = (cx2 - vx1) * s + vy1;
}
//test whether y-coordinates of reduced element still lie in the y-range of the segment
if((cy1 - minc[1] - y * seg_size[1] >= -epsilon
&& cy2 - minc[1] - (y + 1) * seg_size[1] <= epsilon)
|| (cy1 - minc[1] - (y + 1) * seg_size[1] <= epsilon
&& cy2 - minc[1] - y * seg_size[1] >= -epsilon)){
//'it' intersects the current segment; add it to the corresponding structure
grid[x][y][0].push_back(&*it);
cnt++;
}
}
}
}
}
}
//--------------------< test_distances_2d >
//tests whether one of the surface elements in the segment with indices 'ix', 'iy' has a distance to
//'destination' which is smaller than 'min_distance'. If so, 'min_distance' is set to the new smallest
//distance.
void BackGrid::test_distances_2d(const Vertex &destination, int ix, int iy,
double &min_distance, Vertex &hit_vertex){
std::list<Element*> &lst = grid[ix][iy][0];
for(std::list<Element*>::iterator it = lst.begin(), it_end = lst.end(); it != it_end; ++it){
Element &e = **it;
if(e.operation_index != current_operation_index){
e.operation_index = current_operation_index;
Vertex current_hit_vertex;
double distance = distance_point_line_2d(destination,
mesh.vertices[e.v[0]], mesh.vertices[e.v[1]], current_hit_vertex);
if(distance < min_distance){
min_distance = distance;
hit_vertex = current_hit_vertex;
}
}
}
}
//--------------------< distance_to_surface_2d >
double BackGrid::distance_to_surface_2d(const Vertex &destination, Vertex &hit_vertex){
//traverse through segments of the backgrid in a spiral around 'destination'
int sx1 = (int)((destination[0] - minc[0]) / seg_size[0]);
if(sx1 >= segments[0]) sx1 = segments[0] - 1;
if(sx1 < 0) sx1 = 0;
int sx2 = segments[0] - sx1 - 1;
if(sx1 > sx2) sx2 = sx1;
int sy1 = (int)((destination[1] - minc[1]) / seg_size[1]);
if(sy1 >= segments[1]) sy1 = segments[1] - 1;
if(sy1 < 0) sy1 = 0;
int sy2 = segments[1] - sy1 - 1;
if(sy1 > sy2) sy2 = sy1;
if(sx2 > sy2) sy2 = sx2;
//now sy2 stores the maximal number of rectangles around 'destination' that have to be
//considered.
++current_operation_index;
//check the segment containing 'destination' for nearby surface mesh elements
double min_distance = 9e100;
test_distances_2d(destination, sx1, sy1, min_distance, hit_vertex);
//successively check segments around 'destination' with increasing distance to 'destination'
for(int i = 1; i <= sy2; ++i){
int ix, iy;
int p1 = sx1 - i;
if(p1 < 0) p1 = 0;
int p2 = sx1 + i;
if(p2 >= segments[0]) p2 = segments[0]-1;
iy = sy1 - i;
if(iy >= 0){
for(ix = p1; ix <= p2; ++ix){
test_distances_2d(destination, ix, iy, min_distance, hit_vertex);
}
}
iy = sy1 + i;
if(iy < segments[1]){
for(ix = p1; ix <= p2; ++ix){
test_distances_2d(destination, ix, iy, min_distance, hit_vertex);
}
}
p1 = sy1 - i + 1;
if(p1 < 0) p1 = 0;
p2 = sy1 + i - 1;
if(p2 >= segments[1]) p2 = segments[1]-1;
ix = sx1 - i;
if(ix >= 0){
for(iy = p1; iy <= p2; ++iy){
test_distances_2d(destination, ix, iy, min_distance, hit_vertex);
}
}
ix = sx1 + i;
if(ix < segments[0]){
for(iy = p1; iy <= p2; ++iy){
test_distances_2d(destination, ix, iy, min_distance, hit_vertex);
}
}
if(min_distance < i * min_seg_size) break;
}
return min_distance;
}
// ############################
// ## BackGrid - 3D Worlds ##
// ############################
//--------------------< attach_mesh_3d >
void BackGrid::attach_mesh_3d(){
int x, y, z, i, j;
double vx1, vx2, vy1, vy2, vz1, vz2; //sorted coordinates of cube surrounding 'sm'
double sx1, sx2, sy1, sy2, sz1, sz2; //coordinates of the segment
int cx1, cx2, cy1, cy2, cz1, cz2; //starting and ending indices of relevant segments
double ex1, ex2, ey1, ey2, ez1, ez2; //coordinates of reduced edge
double grad_z_for_x, grad_y_for_x, grad_z_for_y; //gradients
double tmp_i0, tmp_j0, tmp_i1, tmp_j1, tmp_i2, tmp_j2; //temporal variables for gradient-calc.
bool edge_intersection;
bool singular_x, singular_y, singular_z;
Vertex diagonal_beg[4], diagonal_end[4];
Vertex sol;
#define DEBUG_OUTPUT_BACKGRID 0
#if DEBUG_OUTPUT_BACKGRID != 0
std::vector<std::vector<std::vector<short int> > > z_bt;
z_bt.resize(segments[0]);
for(x = 0; x < segments[0]; ++x){
z_bt[x].resize(segments[1]);
for(y = 0; y < segments[1]; ++y){
z_bt[x][y].resize(segments[2]);
for(z = 0; z < segments[2]; ++z) z_bt[x][y][z] = 0;
}
}
#endif
for(std::list<Element>::iterator it = mesh.elements.begin(),
it_end = mesh.elements.end(); it != it_end; ++it){
it->operation_index = 0;
//determine coordinates of the minimal axis-parallel cube surrounding 'it'
coordinate_extrema_triangle(mesh, *it, 0, vx1, vx2);
coordinate_extrema_triangle(mesh, *it, 1, vy1, vy2);
coordinate_extrema_triangle(mesh, *it, 2, vz1, vz2);
//determine the segments intersected by that cube
cx1 = (int)((vx1 - minc[0]) / seg_size[0]);
if(cx1 >= segments[0]) cx1 = segments[0] - 1;
cx2 = (int)((vx2 - minc[0]) / seg_size[0]);
if(cx2 >= segments[0]) cx2 = segments[0] - 1;
cy1 = (int)((vy1 - minc[1]) / seg_size[1]);
if(cy1 >= segments[1]) cy1 = segments[1] - 1;
cy2 = (int)((vy2 - minc[1]) / seg_size[1]);
if(cy2 >= segments[1]) cy2 = segments[1] - 1;
cz1 = (int)((vz1 - minc[2]) / seg_size[2]);
if(cz1 >= segments[2]) cz1 = segments[2] - 1;
cz2 = (int)((vz2 - minc[2]) / seg_size[2]);
if(cz2 >= segments[2]) cz2 = segments[2] - 1;
//test for all these segments if an intersection with 'it' occurs
for(z = cz1; z <= cz2; ++z){
for(y = cy1; y <= cy2; ++y){
for(x = cx1; x <= cx2; ++x){
//get coordinates of the current segment
sx1 = x * seg_size[0] + minc[0];
sx2 = sx1 + seg_size[0];
sy1 = y * seg_size[1] + minc[1];
sy2 = sy1 + seg_size[1];
sz1 = z * seg_size[2] + minc[2];
sz2 = sz1 + seg_size[2];
//test whether one of the edges of 'it' intersects the current segment
edge_intersection = false;
for(i = 0; i < 3; ++i){
j = (i + 1) % 3;
tmp_i0 = mesh.vertices[it->v[i]][0];
tmp_j0 = mesh.vertices[it->v[j]][0];
tmp_i1 = mesh.vertices[it->v[i]][1];
tmp_j1 = mesh.vertices[it->v[j]][1];
tmp_i2 = mesh.vertices[it->v[i]][2];
tmp_j2 = mesh.vertices[it->v[j]][2];
//y- and z-gradient for x:
if(fabs(tmp_i0 - tmp_j0) < epsilon){
grad_y_for_x = 0.0;
grad_z_for_x = 0.0;
if(tmp_i0 < sx1 || tmp_i0 > sx2) continue;
}else{
grad_y_for_x = (tmp_i1 - tmp_j1) / (tmp_i0 - tmp_j0);
grad_z_for_x = (tmp_i2 - tmp_j2) / (tmp_i0 - tmp_j0);
}
//z-gradient for y:
if(fabs(tmp_i1 - tmp_j1) < epsilon){
grad_z_for_y = 0.0;
if(tmp_i1 < sy1 || tmp_i1 > sy2) continue;
}else{
grad_z_for_y = (tmp_i2 - tmp_j2) / (tmp_i1 - tmp_j1);
}
//calculate reduction of current edge to the x-range of the current segment
if(tmp_i0 <= tmp_j0){
ex1 = tmp_i0;
ey1 = tmp_i1;
ez1 = tmp_i2;
ex2 = tmp_j0;
ey2 = tmp_j1;
ez2 = tmp_j2;
}else{
ex1 = tmp_j0;
ey1 = tmp_j1;
ez1 = tmp_j2;
ex2 = tmp_i0;
ey2 = tmp_i1;
ez2 = tmp_i2;
}
if(ex1 > sx2 || ex2 < sx1) continue;
if(ex1 < sx1){
ey1 += (sx1 - ex1) * grad_y_for_x;
ez1 += (sx1 - ex1) * grad_z_for_x;
ex1 = sx1;
}
if(ex2 > sx2){
ey2 += (sx2 - ex2) * grad_y_for_x;
ez2 += (sx2 - ex2) * grad_z_for_x;
ex2 = sx2;
}
//calculate reduction of current edge to the y-range of the current segment
if(ey1 <= ey2){
if(ey1 > sy2 || ey2 < sy1) continue;
if(ey1 < sy1){
ez1 += (sy1 - ey1) * grad_z_for_y;
ey1 = sy1;
}
if(ey2 > sy2){
ez2 += (sy2 - ey2) * grad_z_for_y;
ey2 = sy2;
}
}else{
if(ey2 > sy2 || ey1 < sy1) continue;
if(ey2 < sy1){
ez2 += (sy1 - ey2) * grad_z_for_y;
ey2 = sy1;
}
if(ey1 > sy2){
ez1 += (sy2 - ey1) * grad_z_for_y;
ey1 = sy2;
}
}
//test whether z-coordinates of the reduced edge still lie in the z-range of the segment
if((ez1 >= sz1 && ez2 <= sz2) || (ez1 <= sz2 && ez2 >= sz1)){
//'it' intersects the current segment; add it to the corresponding structure
edge_intersection = true;
grid[x][y][z].push_back(&*it);
#if DEBUG_OUTPUT_BACKGRID != 0
z_bt[x][y][z] = 1;
#endif
break;
}
}
//if no edge of 'it' intersects the current segment then test whether one of the diagonals of
//the segment intersects 'it'
if(!edge_intersection){
Vertex p, a, b, c;
//define diagonals
diagonal_beg[0][0] = sx1; diagonal_beg[0][1] = sy1; diagonal_beg[0][2] = sz1;
diagonal_end[0][0] = sx2; diagonal_end[0][1] = sy2; diagonal_end[0][2] = sz2;
diagonal_beg[1][0] = sx2; diagonal_beg[1][1] = sy1; diagonal_beg[1][2] = sz1;
diagonal_end[1][0] = sx1; diagonal_end[1][1] = sy2; diagonal_end[1][2] = sz2;
diagonal_beg[2][0] = sx2; diagonal_beg[2][1] = sy2; diagonal_beg[2][2] = sz1;
diagonal_end[2][0] = sx1; diagonal_end[2][1] = sy1; diagonal_end[2][2] = sz2;
diagonal_beg[3][0] = sx1; diagonal_beg[3][1] = sy2; diagonal_beg[3][2] = sz1;
diagonal_end[3][0] = sx2; diagonal_end[3][1] = sy1; diagonal_end[3][2] = sz2;
//check whether 'it' is singular in some dimension
singular_x = fabs(vx1 - vx2) < epsilon;
singular_y = fabs(vy1 - vy2) < epsilon;
singular_z = fabs(vz1 - vz2) < epsilon;
//check intersections
for(j = 0; j < 4; ++j){
//test whether one of the end-points of current diagonal lies inside the triangle
if(singular_x){
if(fabs(diagonal_beg[j][0] - vx1) < epsilon){
p[0] = diagonal_beg[j][1];
p[1] = diagonal_beg[j][2];
a[0] = mesh.vertices[it->v[0]][1];
a[1] = mesh.vertices[it->v[0]][2];
b[0] = mesh.vertices[it->v[1]][1];
b[1] = mesh.vertices[it->v[1]][2];
c[0] = mesh.vertices[it->v[2]][1];
c[1] = mesh.vertices[it->v[2]][2];
if(point_in_triangle_generous_2d(p, a, b, c)){
edge_intersection = true;
}else{
p[0] = diagonal_end[j][1];
p[1] = diagonal_end[j][2];
if(point_in_triangle_generous_2d(p, a, b, c)) edge_intersection = true;
}
}
}else if(singular_y){
if(fabs(diagonal_beg[j][1] - vy1) < epsilon){
p[0] = diagonal_beg[j][0];
p[1] = diagonal_beg[j][2];
a[0] = mesh.vertices[it->v[0]][0];
a[1] = mesh.vertices[it->v[0]][2];
b[0] = mesh.vertices[it->v[1]][0];
b[1] = mesh.vertices[it->v[1]][2];
c[0] = mesh.vertices[it->v[2]][0];
c[1] = mesh.vertices[it->v[2]][2];
if(point_in_triangle_generous_2d(p, a, b, c)){
edge_intersection = true;
}else{
p[0] = diagonal_end[j][0];
p[1] = diagonal_end[j][2];
if(point_in_triangle_generous_2d(p, a, b, c)) edge_intersection = true;
}
}
}else if(singular_z){
if(fabs(diagonal_beg[j][2] - vz1) < epsilon){
p[0] = diagonal_beg[j][0];
p[1] = diagonal_beg[j][1];
a[0] = mesh.vertices[it->v[0]][0];
a[1] = mesh.vertices[it->v[0]][1];
b[0] = mesh.vertices[it->v[1]][0];
b[1] = mesh.vertices[it->v[1]][1];
c[0] = mesh.vertices[it->v[2]][0];
c[1] = mesh.vertices[it->v[2]][1];
if(point_in_triangle_generous_2d(p, a, b, c)){
edge_intersection = true;
}else{
p[0] = diagonal_end[j][0];
p[1] = diagonal_end[j][1];
if(point_in_triangle_generous_2d(p, a, b, c)) edge_intersection = true;
}
}
}
if(edge_intersection ||
intersection_line_triangle_3d(mesh.vertices[it->v[0]], mesh.vertices[it->v[1]],
mesh.vertices[it->v[2]], diagonal_beg[j], diagonal_end[j], sol)){
grid[x][y][z].push_back(&*it);
#if DEBUG_OUTPUT_BACKGRID != 0
z_bt[x][y][z] = 1;
#endif
break;
}
}
}
}
}
}
}
#if DEBUG_OUTPUT_BACKGRID != 0
std::ofstream out("z_bg.lm3");
out << "dim: 3\n";
out << "size: " << segments[0] << " " << segments[1] << " " << segments[2] << "\n";
out << "box_size: " << seg_size[0] << " " << seg_size[1] << " " << seg_size[2] << "\n";
out << "offset: " << minc[0] << " " << minc[1] << " " << minc[2] << "\n";
out << "y-axis: textual\n";
out << "\nmatrices:\n";
for(z = 0; z < segments[2]; ++z){
out << "z = " << z << ":\n";
for(y = 0; y < segments[1]; ++y){
for(x = 0; x < segments[0] - 1; ++x){
out << z_bt[x][y][z] << " ";
}
out << z_bt[segments[0] - 1][y][z] << "\n";
}
out << "\n";
}
out.close();
#endif
}
//--------------------< test_distances_3d >
//tests whether one of the surface elements in the segment with indices 'ix', 'iy', 'iz' has a
//distance to 'destination' which is smaller than 'min_distance'. If so, 'min_distance' is set
//to the new smallest distance.
void BackGrid::test_distances_3d(const Vertex &destination, int ix, int iy, int iz,
double &min_distance, Vertex &hit_vertex){
std::list<Element*> &lst = grid[ix][iy][iz];
for(std::list<Element*>::iterator it = lst.begin(), it_end = lst.end(); it != it_end; ++it){
Element &e = **it;
if(e.operation_index != current_operation_index){
e.operation_index = current_operation_index;
Vertex current_hit_vertex;
double distance = distance_point_triangle_3d(destination, mesh.vertices[e.v[0]],
mesh.vertices[e.v[1]], mesh.vertices[e.v[2]], current_hit_vertex);
if(distance < min_distance){
min_distance = distance;
hit_vertex = current_hit_vertex;
}
}
}
}
//--------------------< distance_to_surface_3d >
double BackGrid::distance_to_surface_3d(const Vertex &destination, Vertex &hit_vertex){
//traverse through segments of the backgrid in concentric rectangles around 'destination'
int sx1 = (int)((destination[0] - minc[0]) / seg_size[0]);
if(sx1 >= segments[0]) sx1 = segments[0] - 1;
if(sx1 < 0) sx1 = 0;
int sy1 = (int)((destination[1] - minc[1]) / seg_size[1]);
if(sy1 >= segments[1]) sy1 = segments[1] - 1;
if(sy1 < 0) sy1 = 0;
int sz1 = (int)((destination[2] - minc[2]) / seg_size[2]);
if(sz1 >= segments[2]) sz1 = segments[2] - 1;
if(sz1 < 0) sz1 = 0;
//determine maximal number of concentric rectangles around 'destination' that might have
//intersections
int sx2 = segments[0] - sx1 - 1;
if(sx1 > sx2) sx2 = sx1;
int sy2 = segments[1] - sy1 - 1;
if(sy1 > sy2) sy2 = sy1;
if(sx2 > sy2) sy2 = sx2;
int sz2 = segments[2] - sz1 - 1;
if(sz1 > sz2) sz2 = sz1;
if(sy2 > sz2) sz2 = sy2;
//now sz2 stores the maximal number of rectangles
++current_operation_index;
//check the segment containing 'destination' for nearby surface mesh elements
double min_distance = 9e100;
test_distances_3d(destination, sx1, sy1, sz1, min_distance, hit_vertex);
//successively check segments around 'destination' with increasing distance to 'destination'
for(int i = 1; i <= sz2; i++){
int x, y, z;
//determine boundaries of back and front
int p1 = sx1 - i;
if(p1 < 0) p1 = 0;
int p2 = sx1 + i;
if(p2 >= segments[0]) p2 = segments[0] - 1;
int q1 = sy1 - i;
if(q1 < 0) q1 = 0;
int q2 = sy1 + i;
if(q2 >= segments[1]) q2 = segments[1] - 1;
//back
z = sz1 - i;
if(z >= 0){
for(y = q1; y <= q2; y++){
for(x = p1; x <= p2; x++) if(!grid[x][y][z].empty()){
test_distances_3d(destination, x, y, z, min_distance, hit_vertex);
}
}
}
//front
z = sz1 + i;
if(z < segments[2]){
for(y = q1; y <= q2; y++){
for(x = p1; x <= p2; x++) if(!grid[x][y][z].empty()){
test_distances_3d(destination, x, y, z, min_distance, hit_vertex);
}
}
}
//determine boundaries of bottom and top
p1 = sx1 - i;
if(p1 < 0) p1 = 0;
p2 = sx1 + i;
if(p2 >= segments[0]) p2 = segments[0] - 1;
q1 = sz1 - i + 1;
if(q1 < 0) q1 = 0;
q2 = sz1 + i - 1;
if(q2 >= segments[2]) q2 = segments[2] - 1;
//bottom
y = sy1 - i;
if(y >= 0){
for(z = q1; z <= q2; z++){
for(x = p1; x <= p2; x++) if(!grid[x][y][z].empty()){
test_distances_3d(destination, x, y, z, min_distance, hit_vertex);
}
}
}
//top
y = sy1 + i;
if(y < segments[1]){
for(z = q1; z <= q2; z++){
for(x = p1; x <= p2; x++) if(!grid[x][y][z].empty()){
test_distances_3d(destination, x, y, z, min_distance, hit_vertex);
}
}
}
//determine boundaries of left and right
p1 = sy1 - i + 1;
if(p1 < 0) p1 = 0;
p2 = sy1 + i - 1;
if(p2 >= segments[1]) p2 = segments[1] - 1;
q1 = sz1 - i + 1;
if(q1 < 0) q1 = 0;
q2 = sz1 + i - 1;
if(q2 >= segments[2]) q2 = segments[2] - 1;
//left
x = sx1 - i;
if(x >= 0){
for(z = q1; z <= q2; z++){
for(y = p1; y <= p2; y++) if(!grid[x][y][z].empty()){
test_distances_3d(destination, x, y, z, min_distance, hit_vertex);
}
}
}
//right
x = sx1 + i;
if(x < segments[0]){
for(z = q1; z <= q2; z++){
for(y = p1; y <= p2; y++) if(!grid[x][y][z].empty()){
test_distances_3d(destination, x, y, z, min_distance, hit_vertex);
}
}
}
if(min_distance < i * min_seg_size) break;
}
return min_distance;
}