GouMinghao / Geometry3D

Geometry3D: 3D Computational Geometrics Library
GNU General Public License v3.0
88 stars 23 forks source link

Question about point on a plane #5

Closed lihaolin88 closed 3 years ago

lihaolin88 commented 3 years ago

Hello, thank you for your amazing project! Recently I'm trying to use this project to generate an intersection between two polyhedrons, but I found an error: "ValueError: Convex Check Fails Because Point(-39.35859231445235, 261.35056323400175, 129.2997523368219) Is Not On Plane(Point(-27.878684988415753, 258.162898458097, 127.72269713190062), Vector(0.2489692168272217, 0.48178113934707206, 0.8401793039833086))" After look through your source code, I think this error is caused by the precision of the cross between two vectors, for example, if we have three points: a = Point(-39.47570506286104, 247.07101597009853, 122.23513421678561) b = Point(-34.157972900849515, 245.67527429293494, 121.54460938703097) c = Point(-30.09033189524111, 244.60764322271578, 121.01641296281727) if we generate a plane based on this three points, in the code (/geometry/plane) we'll do: vab = b.pv() - a.pv() vac = c.pv() - a.pv() vec1 = vab.cross(vac) plane = Plane(a, vec1) after that, I found "a in plane" is true, but both "b in plane" and "c in plane" will return false because vec1 is not accurate enough, can you help me solve this problem? Thank you so much!

GouMinghao commented 3 years ago

Thank you for your issue! I have run the code and reproduced the problem. It is indeed caused by the precision problem. The vector ab and vector ac are almost parallel. As a result, ab cross ac results in a vector whose normal is almost zero.

Could you please provide more information on why you need to use this library in such an extreme case. After that I can give you suggestions and improve the library.

lihaolin88 commented 3 years ago

Sure, thank you for this very quick reply! The project I'm working on is to find an object from 3D space based on a series of photos taken from different angles for a same object, I already know 3D position of the light source and the object position on each image, so then from each image, with the light source position, I can construct a polyhedron to indicate an area include the object in 3d space, and next, the intersection between different polyhedron will give me a more accurate 3D position about the object; because the light source is a little bit far from the camera, so the polyhedron will become very thin, then the edge of the polyhedron will become almost parallel.

GouMinghao commented 3 years ago

Thanks for your description. I check the case again. The problem is the floaing number precision in calculating the cross result. As a result. the normal vector of plane is not accurate. CPython doesn't handle this problem and it is not a bug of Geometry3D which is not intended to deal with such case.

However, you can deal with the precision problem mannualy such as using Decimal. The example is shown below.

from decimal import Decimal
from Geometry3D import *

def cross(a, b):
    return [a[1] * b[2] - a[2] * b[1],
                a[2] * b[0] - a[0] * b[2],
                a[0] * b[1] - a[1] * b[0]]

a = Point(-39.47570506286104, 247.07101597009853, 122.23513421678561)
b = Point(-34.157972900849515, 245.67527429293494, 121.54460938703097)
c = Point(-30.09033189524111, 244.60764322271578, 121.01641296281727)

da, db, dc = [[Decimal(x) for x in [p.x, p.y, p.z]] for p in [a, b, c]]
ab = [db[i] - da[i] for i in range(3)]
ac = [dc[i] - da[i] for i in range(3)]

# Obtain normal vector of the plane
n = cross(ab, ac)

# Obtain the length of the normal vector
n_norm = sum([x ** 2 for x in n]).sqrt()

# Normalizaze vector
nn = [x / n_norm  for x in n]

# Convert the list to Vector class
vec1 = Vector([float(x) for x in nn])

# Construct the normal vector
plane = Plane(a, vec1)

print('plane:{}'.format(plane))
print('a in plane? {}'.format(a in plane))
print('b in plane? {}'.format(b in plane))
print('c in plane? {}'.format(c in plane))

The result should be

plane:Plane(Point(-39.47570506286104, 247.07101597009853, 122.23513421678561), Vector(-0.1845409513325132, -0.899278020796477, 0.3965396305461311))
a in plane? True
b in plane? True
c in plane? True

I will check if I can modify Geometry3D to solve this problem without adding so much computational overhead.

lihaolin88 commented 3 years ago

Thank you so much for the help!