Determine whether a vertex is in a specified triangle (1)

Problem: Determine whether point P is in the Triangle ABC

 

To determine whether a point is in a triangle, the two most common methods are Area Method and vector homophone method. Although the algorithm is simple, it is not easy to achieve efficiency. We need to consider the differences between two and three dimensions, but also whether the coordinates are represented by floating points or integers.

 

On a two-dimensional plane, the problem is relatively simple. Generally, only 6 multiplication is required. However, the problem of 3D planes is much more complicated. The algorithms seen on the Internet generally require 30 multiplication calculations (if the known point P is on the plane ABC, 21 times is required ). In fact, in a three-dimensional coordinate system, we can increase the number of comparisons and reduce the multiplication calculation to 13 (if the point P is on the plane ABC, we can only perform 8 multiplication calculation at most ).

 

The two most common methods are Area Method and vector same direction method, which are essentially equivalent.

Vector homography: If point P is in a triangle, three vectors are used: AB × ap, ap × ac, pb × pc parallel to the same direction (they are also in the same direction as the vector AB × ac), because these three vectors may be 0, it is quite troublesome to directly judge the parallel direction of them, but considering that AB × ac cannot be 0, it is difficult to directly judge the "vector: AB × ap, ap × ac, and pb × pc are all parallel to AB × ac.

 

Area Method: When point p is within the Triangle abc, the area of the four triangles is equal to abc = abp + apc + pbc.

We can use the vector's vector product to calculate the area abc = | AB × ac |/2.

On the surface, the area of the four triangles is calculated, but according to the following formula:

Ap × ap = 0, pb × pc = (AB-ap) × (ac-ap) = AB × ac-AB × ap-ap × ac

You can calculate less than one vector product.

Formula: | AB × ac | = | AB × ap | + | ap × ac | + | (AB × ac-AB × ap-ap × ac) |

 

For arbitrary vectors a, B, c: | a + B + c | = | a | + | B | + | c | <=> vectors a, B, and c are in the parallel direction.

Therefore, the area method and the vector same direction method are essentially equivalent.

 

 

Next we will first discuss the two-dimensional coordinate system (each point X is considered as the two-dimensional vector OX from the origin O to the point X ).

First, define a two-dimensional Vector Template:

 

Template <typename T> class Vec2 {

T x, y;

Public:

Typedef T value_type;

Vec2 (T xx = 0, T yy = 0): x (xx), y (yy ){};

T cross (const Vec2 & v) const {return x * v. y-y * v. x;} // Vector Product

Vec2 operator-(const Vec2 & v) const {return Vec2 (x-v. x, y-v. y );}

};

 

 

If the coordinates use floating point numbers, it is convenient to take the absolute value of the Floating Point Numbers (with special floating point commands), but there is an error in the size of the comparison between them. The area method is more convenient:

 

Typedef Vec2 <double> Vd2;

 

Bool is_in_triangle (const Vd2 & a, const Vd2 & B, const Vd2 & c, const Vd2 & p)

{

Vd2 AB (B-a), ac (c-a), ap (p-);

 

// Calculate the area using the vector product. The absolute values of the following four values are twice the area of the corresponding triangle,

Double abc = AB. cross (ac );

Double abp = AB. cross (ap );

Double apc = ap. cross (ac );

Double pbc = abc-abp-apc; // equal to pb. cross (pc)

// Area method: The area difference between the four triangles is equal to 0.

Double delta = fabs (abc)-fabs (abp)-fabs (apc)-fabs (pbc );

Return fabs (delta) <DBL_EPSILON;

}

 

 

If the coordinates are represented by integers, the code is relatively troublesome:

 

Typedef Vec2 <int> Vi2;

 

Bool is_in_triangle (const Vi2 & a, const Vi2 & B, const Vi2 & c, const Vi2 & p)

{

Vi2 AB (B-a), ac (c-a), ap (p-);

 

// Calculate the area using the vector product. The absolute values of the following four values are twice the area of the corresponding triangle,

Int abc = AB. cross (ac );

Int abp = AB. cross (ap );

Int apc = ap. cross (ac );

Int pbc = abc-abp-apc; // equal to pb. cross (pc)

 

// Method 1: Area Method: The area difference between the four triangles is equal to 0.

Return abs (abc) = abs (abp) + abs (apc) + abs (pbc)

 

// Method 2: vector co-direction method: abcapc pbc is in the same direction as abc:

If (abc <0) {abp =-abp; apc =-apc; pbc =-pbc ;}

Return (abp> = 0) & (apc> = 0) & (pbc> = 0 );

}

 

Method 1: To calculate four absolute values, it seems that four conditional jumps are required. However, most compilers can use bitwise operations to directly calculate the absolute values (Note: GCC requires additional parameters ), no conditional jump is required.

Method 2: less than method 1, but one more conditional jump.

Which method is more efficient depends on the code generated by the compiler.

 

In the above Code, two optimization methods can be used:

① Take the absolute value of integer x and use bitwise operations:

Set y = 0 (when x> = 0)

=-1 (when x <0)

(The compiler can use commands such as cdq or sar to directly calculate the y value from x)

Abs (x) = (x xor y)-y

Or: = (x + y) xor y

Or: = x-(2 * x & y)

 

② For integers a, B, c, a >=0 & B >=0 & c >=0 is equivalent

(A> = 0) & (B> = 0) & (c> = 0) is equivalent:

(A | B | c)> = 0

 

To prevent the compiler from performing optimization, you can directly perform manual optimization:

 

Inline int chg_sign (int x, int sign) // sign can only be 0 or-1. The functions return x and-x respectively.

{

Return (x + sign) ^ sign;

// Return (x ^ sign)-sign;

}

 

Bool is_in_triangle (const Vi2 & a, const Vi2 & B, const Vi2 & c, const Vi2 & p)

{

Vi2 AB (B-a), ac (c-a), ap (p-);

 

// Calculate the area using the vector product. The absolute values of the following four values are twice the area of the corresponding triangle,

Int abc = AB. cross (ac );

Int abp = AB. cross (ap );

Int apc = ap. cross (ac );

Int pbc = abc-abp-apc; // equal to pb. cross (pc)

 

// Method 3: vector homography (optimized version)

Const int sign = (abc> = 0)-1;

// Const int sign = abc> (sizeof (abc) * CHAR_BIT-1 );

Return (chg_sign (abp, sign) | chg_sign (apc, sign) | chg_sign (pbc, sign)> = 0;

}