Dado un punto $P$ y un Plano triangular $T$ con tres puntos $P_0$ , $P_1$ y $P_2$ Quiero saber si el $P$ se cruza con el $T$ . Me doy cuenta de que hay una pregunta muy similar aquí Sin embargo, no entendí ninguna de las respuestas de esa página, ya que soy bastante nuevo en este tipo de matemáticas. Los otros recursos que leí estaban ilustrados en un plano 2D en lugar de 3D. ¿Existe un método fácil para calcular la intersección de un punto en un plano triangular 3D que se pueda traducir a código (estoy tratando de usar esto en un proyecto de codificación en el que estoy trabajando)?
Respuestas
¿Demasiados anuncios?
mjw
Puntos
225
Dados sus puntos $\{P_1, P_2, P_3\}$ donde $P_j=(x_j,y_j,z_j)$ la ecuación del plano se halla resolviendo un sistema de ecuaciones:
$$\begin{aligned} a x_1 + b y_1 + c z_1 = \rho \\ a x_2 + b y_2 + c z_2 = \rho \\ a x_3 + b y_3 + c z_3 = \rho \\ \end{aligned}$$
para $(a,b,c,\rho)$ . Para ser un plano, al menos uno de $(a,b,c)$ debe ser distinto de cero. Una vez que tenemos una ecuación para el plano, podemos multiplicar cada término por cualquier constante distinta de cero, por lo que los coeficientes sólo son únicos hasta esta constante.
La ecuación de un plano es entonces $$a x + b y + c z = \rho.$$ Dado un punto $\{x_0,y_0,z_0\}$ se encuentra en ese plano si y sólo si $$a x_0 + b y_0 + c z_0 = \rho.$$
None
Puntos
58
Aquí hay dos cuestiones distintas:
-
Cómo calcular la distancia entre un punto y un plano
-
Cómo averiguar si un punto de un plano está dentro de un triángulo de ese plano
Llamemos a los tres vértices del triángulo $\vec{v}_0 = (x_0, y_0, z_0)$ , $\vec{v}_1 = (x_1, y_1, z_1)$ y $\vec{v}_2 = (x_2, y_2, z_2)$ y el punto $\vec{p} = (x_p, y_p, z_p)$ .
La primera pregunta se responde utilizando la ecuación de un plano, $$\vec{n} \cdot \vec{p} = d$$ donde $\vec{n}$ es la normal del plano (perpendicular al plano), y $d$ es la distancia con signo del plano al origen en unidades de la longitud del vector normal. (Si $d = 0$ el plano pasa por el origen, $(0, 0, 0)$ . Si $d \gt 0$ el plano está en la misma dirección que el vector normal desde el origen, y si $d \lt 0$ el plano está en la dirección opuesta).
Podemos calcular fácilmente $\vec{n}$ y $d$ de los tres vértices del triángulo: $$\begin{aligned} \vec{n} &= (\vec{v}_1 - \vec{v}_0) \times (\vec{v}_2 - \vec{v}_0) \\ d &= \vec{n} \cdot \vec{v}_0 = \vec{n} \cdot \vec{v}_1 = \vec{n} \cdot \vec{v}_2 \\ \end{aligned}$$ donde $\times$ representa el producto vectorial cruzado, y $\cdot$ producto vectorial punto. Como podemos escalar $\vec{n}$ libremente a cualquier longitud distinta de cero (sólo cambiará $d$ y el par seguirá representando el mismo plano), es muy útil escalar $\vec{n}$ a la unidad de longitud, $\lVert\vec{n}\rVert = 1$ . (Me gusta denotar estos vectores de longitud unitaria como $\hat{n}$ .)
Una función que calcula el vector normal unitario $\hat{n}$ y la distancia con signo $d$ a partir de tres puntos puede escribirse, por ejemplo, en Python simplificado de la siguiente manera:
from math import sqrt
def Plane(x0,y0,z0, x1,y1,z1, x2,y2,z2):
# g = v1 - v0
gx = x1 - x0
gy = y1 - y0
gz = z1 - z0
# h = v2 - v0
hx = x2 - x0
hy = y2 - y0
hz = z2 - z0
# Cross product n = g × h
nx = gy*hz - gz*hy
ny = gz*hx - gx*hz
nz = gx*hy - gy*hx
# nn = ||n||, Euclidean length
nn = sqrt(nx*nx + ny*ny + nz*nz)
if nn == 0:
raise ValueError("Degenerate triangle: line or point!")
# Scale n to unit length
nx = nx / nn
ny = ny / nn
nz = nz / nn
# Find d - just use any of the three vertices; say v0.
d = nx*x0 + ny*y0 + nz*z0
# Return nx, ny, nz, and d.
return (nx, ny, nz, d)Tenga en cuenta que esto identifica el plano, y sólo necesita ser calculado una vez.
Esencialmente, dado un punto $\vec{p} = (x, y, z)$ el vector normal unitario al plano $\hat{n}$ y la distancia con signo $d$ , $$h = \lvert \hat{n} \cdot \vec{p} - d \rvert$$ indica la distancia entre el punto y el plano.
Cuando trabajamos con números de coma flotante en un programa informático, no podemos comprobar la igualdad exacta, porque los números en sí no son exactos. Debemos incluir un margen de error de redondeo. Esto se hace normalmente utilizando un máquina epsilon un valor que representa el mayor número (digamos, en una coordenada puntual) que consideramos todavía "cero", cuando se tienen en cuenta esos errores de redondeo.
Por lo tanto, para comprobar el punto en el plano, puede utilizar algo como el siguiente código Python simplificado:
def point_in_plane(x,y,z, nx,ny,nz, d, epsilon=0.00001):
return (abs(nx*x + ny*y + nz*z - d) <= epsilon)que devuelve True si el punto está a una distancia épsilon del plano (en unidades de la longitud normal unitaria, por lo que esperamos que la normal del plano sea $\hat{n}$ sea aquí un vector unitario, tal y como viene dado por la función anterior), y False en caso contrario.
Eso soluciona la primera parte.
Si suponemos que ya se ha determinado que el punto está lo suficientemente cerca del plano, podemos volver a utilizar el producto vectorial cruzado para comprobar si el punto está dentro de cualquier polígono convexo (y los triángulos son definitivamente convexos), dados sólo los vértices.
Verás, si tenemos dos vértices sucesivos en el polígono convexo, $\vec{v}_i$ y $\vec{v}_{i+1}$ el producto cruz entre el vector entre los dos y el vector hacia el punto en cuestión, $(\vec{v}_{i+1} - \vec{v}_i)\times(\vec{p} - \vec{v}_i)$ apuntará en la misma dirección que la normal del plano cuando los tres ( $\vec{v}_i$ , $\vec{v}_{i+1}$ y $\vec{p}$ ) están en orden antihorario, y en sentido contrario cuando están en orden horario. (Será un vector nulo cuando los tres sean colineales). Si usamos el producto punto sobre el resultado y la normal del plano, obtendremos un número positivo o negativo dependiendo de si el punto está en sentido antihorario o horario, respecto a la arista. Por lo tanto, ignorando cualquier caso cero, si los casos son todos positivos o todos negativos, pero no una mezcla, sabemos que el punto está rodeado por los bordes, y por lo tanto dentro de (o en) un borde. Si hay una mezcla, entonces sabemos que el punto debe estar fuera del polígono convexo.
Esta función puede escribirse en Python simplificado de la siguiente manera:
def point_in_convex_polygon(p, v, n, epsilon=0.00001):
# Number of counterclockwise and clockwise results
nccw = 0
ncw = 0
vnext = v[len(v)-1 ] # Last vertex
for i in range(0, len(v)): # i = 0, ..., len(v)-1
vprev = vnext
vnext = v[i]
gx = vnext[0] - vprev[0]
gy = vnext[1] - vprev[1]
gz = vnext[2] - vprev[2]
hx = p[0] - vprev[0]
hy = p[1] - vprev[1]
hz = p[2] - vprev[2]
# s = n.(g×h)
s = n[0]*(gy*hz-gz*hy) + n[1]*(gz*hx-gx*hz) + n[2]*(gx*hy-gy*hx)
if s > epsilon:
nccw = nccw + 1
if ncw > 0:
return False
elif s < -epsilon:
ncw = ncw + 1
if nccw > 0:
return False
return TrueTécnicamente, aquí comparamos un producto de coordenadas con épsilon, por lo que matemáticamente deberíamos usar épsilon al cuadrado en su lugar; pero las operaciones que estamos realizando son las que producen el error de redondeo con números de coma flotante, por lo que usar el mismo épsilon en realidad tiene más sentido.
En un programa del mundo real, usted precalcularía $\hat{n}$ y $d$ cada vez que se mueva o añada un triángulo, por razones de eficiencia, pero si la velocidad no es un problema, podría implementar la prueba is-point-in-3D-triangle en Python simplificado como
def point_in_3d_triangle(p, v0, v1, v2, epsilon=0.0001):
# Plane unit normal and signed distance from origin
nx, ny, nz, d = Plane(v0[0], v0[1], v0[2],
v1[0], v1[1], v1[2],
v2[0], v2[1], v2[2])
# Point on the plane?
if not point_in_plane(p[0], p[1], p[2], nx, ny, nz, d, epsilon):
return False
# Point within the triangle?
return point_in_convex_polygon(p, [v0, v1, v2], n, epsilon)Para los triángulos (y paralelogramos), existe una prueba aún mejor, mediante coordenadas baricéntricas $(u, v)$ que corresponden al punto $$\vec{p} = \vec{v}_0 + u \bigr( \vec{v}_1 - \vec{v}_0 \bigr) + v \bigr( \vec{v}_2 - \vec{v}_0 \bigr) \tag{1a}\label{None1a}$$ es decir $$\left\lbrace ~ \begin{aligned} x &= x_0 + u ( x_1 - x_0 ) + v ( x_2 - x_0 ) \\ y &= y_0 + u ( y_1 - y_0 ) + v ( y_2 - y_0 ) \\ z &= z_0 + u ( z_1 - z_0 ) + v ( z_2 - z_0 ) \\ \end{aligned} \right . \tag{1b}\label{None1b}$$ Ahora, $(u, v)$ están dentro del triángulo si y sólo si $$\begin{array}{c} 0 \le u \le 1 \\ 0 \le v \le 1 \\ 0 \le u + v \le 1 \\ \end{array} \tag{2}\label{None2}$$ El caso es que tenemos tres ecuaciones, pero sólo dos incógnitas. Esencialmente, tenemos que excluir el eje con el ángulo más pequeño al plano como irrelevante.
Para cada triángulo, o después de mover un triángulo, calcula $$\begin{aligned} c_{xy} &= x_0 ( y_2 - y_1 ) + x_1 ( y_0 - y_2 ) + x_2 ( y_1 - y_0 ) \\ c_{xz} &= x_0 ( z_2 - z_1 ) + x_1 ( z_0 - z_2 ) + x_2 ( z_1 - z_0 ) \\ c_{yz} &= y_0 ( z_2 - z_1 ) + y_1 ( z_0 - z_2 ) + y_2 ( z_1 - z_0 ) \\ \end{aligned} \tag{3}\label{None3}$$ y elige la de mayor magnitud. Entonces,
$$\begin{aligned} u &= \displaystyle \frac{(y_0 - y_2) x + (x_2 - x_0) y + (x_0 y_2 - x_2 y_0)}{c_{xy}} \\ ~ &= \displaystyle \frac{(z_0 - z_2) x + (x_2 - x_0) z + (x_0 z_2 - x_2 z_0)}{c_{xz}} \\ ~ &= \displaystyle \frac{(z_0 - z_2) y + (y_2 - y_0) z + (y_0 z_2 - z_2 y_0)}{c_{yz}} \\ \end{aligned} \tag{4a}\label{None4a}$$ y $$\begin{aligned} v &= \displaystyle \frac{(y_1 - y_0) x + (x_0 - x_1) y + (x_1 y_0 - x_0 y_1)}{c_{xy}} \\ ~ &= \displaystyle \frac{(z_1 - z_0) x + (x_0 - x_1) z + (x_1 z_0 - x_0 z_1)}{c_{xz}} \\ ~ &= \displaystyle \frac{(z_1 - z_0) y + (y_0 - y_1) z + (y_1 z_0 - y_0 z_1)}{c_{yz}} \\ \end{aligned} \tag{4b}\label{None4b}$$
Sin embargo, como sólo necesitamos utilizar una de las tres anteriores, y son funciones lineales de las coordenadas del punto $\vec{p} = (x, y, z)$ podemos utilizar $$\left\lbrace ~ \begin{aligned} u &= \vec{p} \cdot \vec{u} + u_0 \\ v &= \vec{p} \cdot \vec{v} + v_0 \\ \end{aligned} \right . \quad \iff \quad \left\lbrace ~ \begin{aligned} u &= x u_x + y u_y + z u_z + u_0 \\ v &= x v_x + y v_y + z v_z + v_0 \\ \end{aligned} \right . \tag{5}\label{None5}$$ Ah-ha! Esto será seriamente eficiente. (Pero espero que a estas alturas puedas ver por qué he tenido que escribir semejante muro de texto para explicarlo: sólo con ver esto último me habría tenido al menos a mí mismo rascándome la cabeza con un ¿Qué?, a menos que entendiera cómo se llegó a esto, y que este es el mejor método).
Lo que tenemos que hacer en un programa de ordenador, es cada vez que se cambia el triángulo, calcular $c_{xy}$ , $c_{xz}$ y $c_{yz}$ .
Si $\lvert c_{xy} \rvert \ge \lvert c_{xz} \rvert$ y $\lvert c_{xy} \rvert \ge \lvert c_{yz} \rvert$ definir $$\begin{aligned} u_x &= \displaystyle \frac{y_0 - y_2}{c_{xy}} \\ u_y &= \displaystyle \frac{x_2 - x_0}{c_{xy}} \\ u_z &= 0 \\ u_0 &= \displaystyle \frac{x_0 y_2 - x_2 y_0}{c_{xy}} \\ \end{aligned} \quad \text{and} \quad \begin{aligned} v_x &= \displaystyle \frac{y_1 - y_0}{c_{xy}} \\ v_y &= \displaystyle \frac{x_0 - x_1}{c_{xy}} \\ v_z &= 0 \\ v_0 &= \displaystyle \frac{x_1 y_0 - x_0 y_1}{c_{xy}} \\ \end{aligned} \tag{6xy}\label{None6xy}$$ En caso contrario, si $\lvert c_{xz} \rvert \ge \lvert c_{xy} \rvert$ y $\lvert c_{xz} \rvert \ge \lvert c_{yz} \rvert$ definir $$\begin{aligned} u_x &= \displaystyle \frac{z_0 - z_2}{c_{xz}} \\ u_y &= 0 \\ u_z &= \displaystyle \frac{x_2 - x_0}{c_{xz}} \\ u_0 &= \displaystyle \frac{x_0 z_2 - x_2 z_0}{c_{xz}} \\ \end{aligned} \quad \text{and} \quad \begin{aligned} v_x &= \displaystyle \frac{z_1 - z_0}{c_{xz}} \\ v_y &= 0 \\ v_z &= \displaystyle \frac{x_0 - x_1}{c_{xz}} \\ v_0 &= \displaystyle \frac{x_1 z_0 - x_0 z_1}{c_{xz}} \\ \end{aligned} \tag{6xz}\label{None6xz}$$ Por lo demás, $\lvert c_{yz} \rvert \ge \lvert c_{xy} \rvert$ y $\lvert c_{yz} \rvert \ge \lvert c_{xz} \rvert$ y puede definir $$\begin{aligned} u_x &= 0 \\ u_y &= \displaystyle \frac{z_0 - z_2}{c_{yz}} \\ u_z &= \displaystyle \frac{y_2 - y_0}{c_{yz}} \\ u_0 &= \displaystyle \frac{y_0 z_2 - y_2 z_0}{c_{yz}} \\ \end{aligned} \quad \text{and} \quad \begin{aligned} v_x &= 0 \\ v_y &= \displaystyle \frac{z_1 - z_0}{c_{xz}} \\ v_z &= \displaystyle \frac{y_0 - y_1}{c_{xz}} \\ v_0 &= \displaystyle \frac{y_1 z_0 - y_0 z_1}{c_{xz}} \\ \end{aligned} \tag{6yz}\label{None6yz}$$
Aquí hay un dominio público (Creative Commons Zero) - con licencia de ejemplo adecuado en Python3 que implementa una clase Triángulo (y un ayudante de la clase Vector de aplicación básica de álgebra vectorial 3D) optimizado para la prueba de containmen:
# SPDX-License-Identifier: CC0-1.0
from math import sqrt
class Triangle(tuple):
"""3D Triangle class optimized for point-in-triangle testing"""
def __new__(cls, a, b, c):
a = Vector(a[0], a[1], a[2])
b = Vector(b[0], b[1], b[2])
c = Vector(c[0], c[1], c[2])
# First, calculate the unit normal vector (cross product).
n = ((b - a) ^ (c - a)).unit
if not n:
raise ValueError("Degenerate triangle")
# Calculate the signed distance from origin (dot product).
d = n | a
# Calculate the three possible divisors.
div_xy = a.x*(c.y-b.y) + b.x*(a.y-c.y) + c.x*(b.y-a.y)
div_xz = a.x*(c.z-b.z) + b.x*(a.z-c.z) + c.x*(b.z-a.z)
div_yz = a.y*(c.z-b.z) + b.y*(a.z-c.z) + c.y*(b.z-a.z)
abs_xy = abs(div_xy)
abs_xz = abs(div_xz)
abs_yz = abs(div_yz)
if abs_xy >= abs_xz and abs_xy >= abs_yz:
# d_xy has the largest absolute value; using xy plane
u_axis = Vector((a.y-c.y)/div_xy, (c.x-a.x)/div_xy, 0)
v_axis = Vector((b.y-a.y)/div_xy, (a.x-b.x)/div_xy, 0)
u_offset = (a.x*c.y - a.y*c.x)/div_xy
v_offset = (a.y*b.x - a.x*b.y)/div_xy
elif abs_xz >= abs_xy and abs_xz >= abs_yz:
# d_xz has the largest absolute value; using xz plane
u_axis = Vector((a.z-c.z)/div_xz, 0, (c.x-a.x)/div_xz)
v_axis = Vector((b.z-a.z)/div_xz, 0, (a.x-b.x)/div_xz)
u_offset = (a.x*c.z - a.z*c.x)/div_xz
v_offset = (a.z*b.x - a.x*b.z)/div_xz
else:
# d_yz has the largest absolute value; using yz plane
u_axis = Vector(0, (a.z-c.z)/div_yz, (c.y-a.y)/div_yz)
v_axis = Vector(0, (b.z-a.z)/div_yz, (a.y-b.y)/div_yz)
u_offset = (a.y*c.z - a.z*c.y)/div_yz
v_offset = (a.z*b.y - a.y*b.z)/div_yz
return tuple.__new__(cls, (a, b, c, n, d, u_axis, u_offset, v_axis, v_offset))
@property
def vertex(self):
"""Tuple of all three vertices"""
return (self[0], self[1], self[2])
@property
def normal(self):
"""Triangle plane unit normal vector"""
return self[3]
@property
def d(self):
"""Triangle plane signed distance from origin"""
return self[4]
@property
def plane(self):
"""Triangle plane (nx, ny, nz, d)"""
return self[3][0], self[3][1], self[3][2], self[4]
@property
def u_axis(self):
"""U axis vector"""
return self[5]
@property
def u_offset(self):
"""U axis offset"""
return self[6]
@property
def v_axis(self):
"""V axis vector"""
return self[7]
@property
def v_offset(self):
"""V axis offset"""
return self[8]
def distance_to_point(self, p):
"""Distance between point p and the triangle plane"""
return abs((Vector(p[0], p[1], p[2]) | self.normal) - self.d)
def contains(self, p, epsilon=0.0001):
"""Returns (u,v) coordinates of point p if it is inside this triangle, None otherwise."""
# Make sure p is a Vector.
p = Vector(p[0], p[1], p[2])
# Outside the plane of the triangle?
if abs((p | self.normal) - self.d) > epsilon:
return None
# U coordinate. We can use exact zero for u,v coordinates.
u = (p | self.u_axis) + self.u_offset
if u < 0 or u > 1:
return None
# V coordinate.
v = (p | self.v_axis) + self.v_offset
if v < 0 or v > 1:
return None
# This is a triangle, not a parallelogram; u+v <= 1. (We know u+v >= 0.)
if u + v > 1:
return None
# Yes, p is within this triangle.
return (u, v)
class Vector(tuple):
"""3D Vector class with basic vector algebra support"""
def __new__(cls, x, y, z):
return tuple.__new__(cls, (float(x), float(y), float(z)))
def __str__(self):
return "(%s,%s,%s)" % (str(self[0]), str(self[1]), str(self[2]))
def __bool__(self):
"""Nonzero vectors are True, zero vectors False."""
return (self[0]*self[0] + self[1]*self[1] + self[2]*self[2] > 0.0)
def __abs__(self):
"""abs(v): Absolute value is Euclidean length"""
return sqrt(self[0]*self[0] + self[1]*self[1] + self[2]*self[2])
def __neg__(self):
"""-v: Vector negation"""
return tuple.__new__(self.__class__, (-self[0], -self[1], -self[2]))
def __pos__(self):
"""+v: Vector itself"""
return self
def __add__(self, other):
"""v+w: Vector addition (left side term)"""
return tuple.__new__(self.__class__, (self[0]+other[0], self[1]+other[1], self[2]+other[2]))
def __radd__(self, other):
"""w+v: Vector addition (right side term)"""
return tuple.__new__(self.__class__, (other[0]+self[0], other[1]+self[1], other[2]+self[2]))
def __sub__(self, other):
"""v-w: Vector subtraction (left side term)"""
return tuple.__new__(self.__class__, (self[0]-other[0], self[1]-other[1], self[2]-other[2]))
def __rsub__(self, other):
"""w-v: Vector subtraction (right side term)"""
return tuple.__new__(self.__class__, (other[0]-self[0], other[1]-self[1], other[2]-self[2]))
def __mul__(self, scalar):
"""v*n: Vector-scalar multiplication (left side term)"""
if isinstance(scalar, (int, float)):
return tuple.__new__(self.__class__, (self[0]*scalar, self[1]*scalar, self[2]*scalar))
else:
return NotImplemented
def __rmul__(self, scalar):
"""n*v: Vector-scalar multiplication (right side term)"""
if isinstance(scalar, (int, float)):
return tuple.__new__(self.__class__, (scalar*self[0], scalar*self[1], scalar*self[2]))
else:
return NotImplemented
def __truediv__(self, scalar):
"""v/n: Vector-scalar division (left side term)"""
if isinstance(scalar, (int, float)):
return tuple.__new__(self.__class__, (self[0]/scalar, self[1]/scalar, self[2]/scalar))
else:
return NotImplemented
def __or__(self, other):
"""v|w: Vector dot product (left side term)"""
if isinstance(other, (tuple, list)) and len(other) >= 3:
return self[0]*other[0] + self[1]*other[1] + self[2]*other[2]
else:
return NotImplemented
def __ror__(self, other):
"""w|v: Vector dot product (right side term)"""
if isinstance(other, (tuple, list)) and len(other) >= 3:
return other[0]*self[0] + other[1]*self[1] + other[2]*self[2]
else:
return NotImplemented
def __xor__(self, other):
"""v^w: Vector cross product (left side term)"""
if isinstance(other, (tuple, list)) and len(other) >= 3:
return tuple.__new__(self.__class__, (self[1]*other[2] - self[2]*other[1],
self[2]*other[0] - self[0]*other[2],
self[0]*other[1] - self[1]*other[0]))
else:
return NotImplemented
def __rxor__(self, other):
"""w^v: Vector dot product (right side term)"""
if isinstance(other, (tuple, list)) and len(other) >= 3:
return tuple.__new__(self.__class__, (other[1]*self[2] - other[2]*self[1],
other[2]*self[0] - other[0]*self[2],
other[0]*self[1] - other[1]*self[0]))
else:
return NotImplemented
def __and__(self, notimplemented): return NotImplemented
def __rand__(self, notimplemented): return NotImplemented
def __rtruediv__(self, notimplemented): return NotImplemented
@property
def x(self):
return self[0]
@property
def y(self):
return self[1]
@property
def z(self):
return self[2]
@property
def norm(self):
"""Absolute value or Euclidean length"""
return sqrt(self[0]*self[0] + self[1]*self[1] + self[2]*self[2])
@property
def normsqr(self):
"""v|v: Euclidean length squared"""
return self[0]*self[0] + self[1]*self[1] + self[2]*self[2]
@property
def unit(self):
"""Vector scaled to unit length"""
n = sqrt(self[0]*self[0] + self[1]*self[1] + self[2]*self[2])
return tuple.__new__(self.__class__, (self[0]/n, self[1]/n, self[2]/n))
if __name__ == '__main__':
from random import Random
from sys import exit
rng = Random()
# Coordinate range for the random points.
Cmin = -9
Cmax = +9
N = 100
P = 500
epsilon = 0.0001
for i in range(0, N):
# Generate three random vertices, at least 0.01 distance away from each other.
while True:
v0 = Vector(rng.uniform(Cmin,Cmax), rng.uniform(Cmin,Cmax), rng.uniform(Cmin,Cmax))
v1 = Vector(rng.uniform(Cmin,Cmax), rng.uniform(Cmin,Cmax), rng.uniform(Cmin,Cmax))
v2 = Vector(rng.uniform(Cmin,Cmax), rng.uniform(Cmin,Cmax), rng.uniform(Cmin,Cmax))
if (v1-v0).norm >= 0.01 and (v2-v0).norm >= 0.01 and (v2-v1).norm >= 0.01:
break
tri = Triangle(v0, v1, v2)
# Unit normal vector to the triangle plane
n = ((v1 - v0) ^ (v2 - v0)).unit
for j in range(0, P):
# Generate a random point (u, v) within a triangle
while True:
u = rng.uniform(0, 1)
v = rng.uniform(0, 1)
if u + v <= 1:
break
# Calculate its position.
p = v0 + u*(v1 - v0) + v*(v2 - v0)
# Ask the Triangle class for the (u,v) coordinates.
c = tri.contains(p)
if c is None:
print("Error: %s is at (%.3f,%.3f) inside the triangle %s, %s, %S, but .contains() failed." % (str(p), u, v, str(v0), str(v1), str(v2)))
exit(1)
# Find q definitely outside the triangle
if Cmin < 0.1 or Cmax > 0.1:
while True:
r = rng.uniform(Cmin, Cmax)
if r < -0.1 or r > 0.1:
break
q = p + r*n
c = tri.contains(q)
if c is not None:
print("Error: %s is outside the triangle %s, %s, %s, but .contains returned %s." % (str(q), str(v0), str(v1), str(v2), str(c)))
exit(1)
print("Tested %d random triangles and %d points inside and outside successfully." % (N, P))Si lo guarda como triangle.py puede consultar la documentación de su interfaz utilizando pydoc3 triangle en el mismo directorio. Puedes ejecutarlo para hacer autocomprobaciones con triángulos y puntos aleatorios (tanto interiores como exteriores). Incluso puedes incluirlo en tus propios proyectos poniéndolo en el mismo directorio, y añadiendo from triangle import Vector, Triangle a su guión.
Después de crear una instancia de Triangle puede utilizar su clase .contains(point) para comprobar si point está dentro del triángulo o no. En caso contrario, devolverá None de lo contrario $(u,v)$ Coordenadas baricéntricas a point dentro del triángulo.
