#include <iostream>
#include <vector>
struct Vector {
float x, y, z;
Vector(float x, float y, float z) : x(x), y(y), z(z) {
}
};
Vector &operator+=(Vector &v1, Vector v2) {
v1.x += v2.x;
v1.y += v2.y;
v1.z += v2.z;
return v1;
}
Vector operator-(Vector v1, Vector v2) {
return Vector(v1.x-v2.x, v1.y-v2.y, v1.z-v2.z);
}
Vector operator*(Vector v, float s) {
return Vector(v.x*s, v.y*s, v.z*s);
}
Vector cross_product(Vector v1, Vector v2) {
return Vector(v1.y*v2.z-v1.z*v2.y, v1.z*v2.x-v1.x*v2.z, v1.x*v2.y-v1.y*v2.x);
}
struct Triangle {
Vector a, b, c;
Triangle(Vector a, Vector b, Vector c) : a(a), b(b), c(c) {
}
};
Vector intersection_with_z0(Vector p, Vector q) {
float t = p.z / (p.z - q.z);
return Vector(p.x + t * (q.x - p.x), p.y + t * (q.y - p.y), 0.0f);
}
// Keeps z<0 part
std::vector<Triangle> clip(std::vector<Triangle> const &triangles) {
std::vector<Triangle> result;
for (Triangle const &triangle : triangles) {
int flags = (triangle.a.z > 0.0f) + 2 * (triangle.b.z > 0.0f) + 4 * (triangle.c.z > 0.0f);
switch (flags) {
case 0:
result.push_back(triangle);
break;
case 1:
result.emplace_back(intersection_with_z0(triangle.a, triangle.b), triangle.b, triangle.c);
result.emplace_back(intersection_with_z0(triangle.a, triangle.b), triangle.c,
intersection_with_z0(triangle.c, triangle.a));
break;
case 2:
result.emplace_back(intersection_with_z0(triangle.b, triangle.c), triangle.c, triangle.a);
result.emplace_back(intersection_with_z0(triangle.b, triangle.c), triangle.a,
intersection_with_z0(triangle.a, triangle.b));
break;
case 3:
result.emplace_back(triangle.c, intersection_with_z0(triangle.c, triangle.a),
intersection_with_z0(triangle.b, triangle.c));
break;
case 4:
result.emplace_back(intersection_with_z0(triangle.c, triangle.a), triangle.a, triangle.b);
result.emplace_back(intersection_with_z0(triangle.c, triangle.a), triangle.b,
intersection_with_z0(triangle.b, triangle.c));
break;
case 5:
result.emplace_back(triangle.b, intersection_with_z0(triangle.b, triangle.c),
intersection_with_z0(triangle.a, triangle.b));
break;
case 6:
result.emplace_back(triangle.a, intersection_with_z0(triangle.a, triangle.b),
intersection_with_z0(triangle.c, triangle.a));
break;
case 7:
break;
}
}
return result;
}
Vector force_due_to_hydrostatic_pressure(Triangle const &triangle, float gravity_times_density) {
return cross_product(triangle.b - triangle.c, triangle.a - triangle.c) *
((triangle.a.z + triangle.b.z + triangle.c.z) / 3 * gravity_times_density);
}
Vector buoyancy(std::vector<Triangle> const &mesh, float gravity_times_density) {
Vector result(0.0f, 0.0f, 0.0f);
std::vector<Triangle> clipped_mesh = clip(mesh);
for (Triangle const &triangle : clipped_mesh)
result += force_due_to_hydrostatic_pressure(triangle, gravity_times_density);
return result;
}
int main() {
Vector a(0.0f, 0.0f, +1.0f);
Vector b(0.0f, 0.0f, -1.0f);
Vector c(1.0f, 0.0f, -1.0f);
Vector d(0.0f, 1.0f, -1.0f);
std::vector<Triangle> mesh = {
Triangle(a, b, c),
Triangle(b, a, d),
Triangle(c, b, d),
Triangle(a, c, d)
};
Vector force = buoyancy(mesh, 9.80665f);
std::cout << '(' << force.x << ',' << force.y << ',' << force.z << ')' << '\n';
}
I haven't written code to clip triangles before, but the way I did it seems reasonable. You can make the code a bit faster by skipping the part that creates the clipped mesh: You can just compute the force due to hydrostatic pressure directly on each triangle and not store the triangles anywhere. I just thought the code would be more clear this way, and it's probably already not going to be a bottleneck in your program. You could also compute only the z coordinate of the force, because the x and y components should cancel out.
Create triangles out of any set of points
Although a convex hull would not solve the problem you specifically mentioned. It does a damn good job at collecting all the external points. You can then use a trinagalization algorithm to compute triangles, which you can then create an area from for your buoyancy problem.
I don't think that's right. The convex hull problem is typically understood to require not just selecting the points, but actually giving you a triangle list.
However, I don't think there is any need to compute this at all. The OP says the problem he really needs to solve is computing the submerged volume of an object. I suspect what he really needs is to compute buoyancy, which I haven't done before, but I am sure can be done by adding the pressure at each face of the object, and that would be much easier to do.
The convex hull problem in 2D space is not concerned with a triangle list. but a list of edges. Which can be trangalized.
In 3D space, the convex hull problem is usually solved in several passes. Off the top of my head, this is my approach. First pass is looking at a set of points between x height and projecting them to a plane in the middle, then scan each point. Continue this for the entire model. Then you solve the second part in a second pass. You now have a graph of points that are assumed to be mostly squares with the possibilities of triangles. Ignore the triangles, carve the squares into a triangle. With a connectivity graph, you can then run an algorithm to determine a proper mesh.
The other option is that if you already have a complete mesh, then you can just merge vertices together to create your convex hull.
I forgot to reply alvaros post: i am not convincied about method of applying pressure force at center of each triangle, this will require additional hull subdivision.
https://sites.google.com/site/customprog/
i need somehow to grou that to get
Looks like exact task for Delaunay Triangulation, and there's lots of libs for it.