See whether a point is in a triangle.

I've been trying to solve a problem which states the following: Write a program, which, given the coordinates of a point, and the coordinates of the 2 remaining points of a triangle, find whether it is inside the triangle or not. (one point of the triangle is always 0, 0)

Example picture:
http://imgur.com/GOZoTXO

Let's say input is like this:
2 2 (coordinates of the point, let's say point A for easy reference)
1 5 3 1 (coordinates of the 2 remaining points of the triangle, B and C for easy reference)

I've thought of finding the slope between the point (0, 0) and point A. Also the equation of point (0, 0) with point A, and then the intercept point of point A with line point B with point C. If slope of point (0,0) and A is bigger than slope (B and (0,0)) and smaller than slope (C and (0,0)) and also the y value of the intercept point between lines (0,0, A) and (B, C) is bigger that y of A, the n it would be inside the triangle. So, I wrote the following program.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
  #include<iostream>
#include<algorithm>
#include<string>

using namespace std;

const double INF = 1000000.00;  //in case line is in the form of x = number

class Point{         //declare class for point
public:
double x, y;
void create(double x1, double y1){
x = x1;
y = y1;
}
void setOrigin(){
x = 0;
y = 0;
}
};

class Triangle{           //declare class for Triangle
public: 
Point a, b;
void create(Point a1, Point b1){
a = a1;
b = b1;
}
};

class Equation{   //declare class of equation 
public:
double slope;          
double intercept;      //intercept point, meaning b in y=ax + b
void create(double slope1, double intercept1){ 
slope = slope1;
intercept = intercept1;
}
double returnY(double x){
double y  = slope*x + intercept;
}
};

double findSlope (Point a, Point b){  //find slope in a line of a and b
double slope;
double Dy = a.y - b.y;
double Dx = a.x - b.x;
if (Dx==0){
    slope = INF;
}
else{
    slope = Dy / Dx;
}
return slope;
}

double findIntercept (Point a, Point b){  //find Intercept (b in y=ax+b)
double slope = findSlope (a, b);
double intercept = a.y - slope*a.x;
return intercept;
}

Equation findEquation(Point a, Point b){  //find equation of a line
Equation result;
double slope = findSlope(a, b);
double intercept = findIntercept (a, b);
result.create(slope, intercept);
return result;
}

Point findMeetPoint(Equation a, Equation b){     //find the intercept point between to lines
Point result;
double x = (a.intercept - b.intercept) / (b.slope - a.slope);
double y = a.returnY(x);
result.create(x, y);
}

bool isInTriangle(Triangle current, Point check){  //finally check if it's in the triangle or not
bool is = false;
Point origin;
origin.setOrigin();
double slopeEdge = findSlope(current.a, origin);
double slopeCheck = findSlope(check, origin);
double slopeBase = findSlope(current.a, origin);
Equation thirdLine = findEquation(current.a, current.b);
Equation checkLine = findEquation(check, origin);
Point meet = findMeetPoint(checkLine, thirdLine);
if (slopeEdge>slopeCheck && slopeBase<slopeCheck && findMeetPoint(findEquation(check, origin), findEquation(current.a,current.b).y>check.y) {
    is = true;
}
return is;
}

int main(){  //read input and print answer
int x, y;
Point check, a, b;
Triangle ab;
cin>>x>>y;
check.create(x, y);
cin>>x>>y;
a.create(x, y);
cin>>x>>y;
b.create(x, y);
ab.create(a, b);
if (isInTriangle(ab, check)){
    cout<<"Inside";
}
else{
    cout<<"Outside";
}
}
Last edited on
Courtesy of the barycentric coordinate system:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
#include <iostream>

struct Point {
	Point() : x(), y() {}
	Point(int _x, int _y) : x(_x), y(_y) {}
	int x, y;
};

struct Triangle {
	const static enum Index {First = 0, Second, Third, num_points};
	::Point points[Index::num_points] = {};

	bool isPointInside(::Point point);
};

bool Triangle::isPointInside(::Point point) {
	double denominator = ((points[Index::Second].y - points[Index::Third].y) * (points[Index::First].x - points[Index::Third].x) + (points[Index::Third].x - points[Index::Second].x) * (points[Index::First].y - points[Index::Third].y));
	double a = ((points[Index::Second].y - points[Index::Third].y) * (point.x - points[Index::Third].x) + (points[Index::Third].x - points[Index::Second].x) * (point.y - points[Index::Third].y)) / denominator;
	double b = ((points[Index::Third].y - points[Index::First].y) * (point.x - points[Index::Third].x) + (points[Index::First].x - points[Index::Third].x) * (point.y - points[Index::Third].y)) / denominator;
	double c = 1 - a - b;

	return ((0 <= a && a <= 1) && (0 <= b && b <= 1) && (0 <= c && c <= 1));
}

int main() {

	Triangle triangle;

	triangle.points[Triangle::Index::First] = { 1, 1 };
	triangle.points[Triangle::Index::Second] = { 4, 2 };
	triangle.points[Triangle::Index::Third] = { 2, 7 };

	::Point point_one = { 2, 3 };
	::Point point_two = { 1, 5 };

	std::cout << std::boolalpha;
	std::cout << "Point #1:\t" << triangle.isPointInside(point_one) << std::endl;
	std::cout << "Point #2:\t" << triangle.isPointInside(point_two) << std::endl;

	return 0;
}


Note, that this is just one way of doing it. It's also not 100% accurate, due to the inaccurate nature of floating-point numbers.
closed account (D80DSL3A)
Solved already? But, both versions have pesky potential division by zero problems (B,C both on x or y axis for eg).
Courtesy vector algebra:
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
// vector cross product UxV
double xProd( double Ux, double Uy, double Vx, double Vy ) { return Ux*Vy - Uy*Vx; }

bool isAinOBC( double Ax, double Ay, double Bx, double By, double Cx, double Cy )
{
    // require OAxOB, OAxOC opposite signs
    double OAxOB = xProd( Ax,Ay,Bx,By );
    double OAxOC = xProd( Ax,Ay,Cx,Cy );
    if( OAxOB/OAxOC > 0.0 ) return false;

    // require BAxBO, BAxBC opposite signs
    double BAxBO = xProd( Ax-Bx, Ay-By, -Bx, -By );
    double BAxBC = xProd( Ax-Bx, Ay-By, Cx-Bx, Cy-By );
    return BAxBO/BAxBC < 0.0;
}

int main()
{
    double bx = 15, by = 7;
    double cx = 5, cy = 13;
    char rpt = 'n';

    do// enter Ax and Ay to test if in triangle OBC
    {
        double ax, ay;
        cout << "A: ";
        cin >> ax >> ay;

        if( isAinOBC(ax,ay,bx,by,cx,cy) )
            cout << "A is in OBC\n";
        else
            cout << "A isn\'t in OBC\n";

        cout << "again?(y/n) "; cin >> rpt;
    }while( rpt == 'y' );

    return 0;
}

Actually, considering the thread title I feel justified in making sure this method gets presented.

Very brief justification... due to the high likelihood that I'm necro-posting again.
The cross product is similar to a slope comparison. Your "intermediate slope" requirement is being followed here.
Limit consideration to 1st quadrant for this simple proof (solution applies generally), so consider all Ax, Bx, Cx > 0.
You require:
slope of OB < slope of OA < slope of OC
or
By/Bx < Ay/Ax < Cy/Cx
splitting into 2 inequalities
Left: By/Bx < Ay/Ax
multiply both sides by Bx*Ax which is presumed > 0 to get
Ax*By < Ay*Bx --> Ax*By - Ay*Bx < 0
Define cross product of 2 vectors U,V as UxV = Ux*Vy - Uy*Vx
Then the fact that slope OB < slope OA implies AxB < 0
Repeat for the right inequality to get
AxC > 0
Generally (ie. no longer restricted to 1st quadrant) require that the signs are 1 positive, 1 negative, or:
AxB/AxC < 0
This accounts for 1st 3lines in isAinOAB function.
2nd set of 3 lines is repeat of test at vertex B

dit: Simplify code in main()
Last edited on
My method actually had a handling of division by zero.
1
2
3
4
5
6
7
8
9
10
11
12
double findSlope (Point a, Point b){  //find slope in a line of a and b
double slope;
double Dy = a.y - b.y;
double Dx = a.x - b.x;
if (Dx==0){
    slope = INF; //HERE
}
else{
    slope = Dy / Dx;
}
return slope;
}




I believe where my code was actually mistaken it was because I always took as a fact that the comparisson would always be like this
slopeEdge>slopeCheck && slopeBase<slopeCheck
and I didn't check whether slopeEdge was bigger than slopeBase, to make the correct comparisson, but again in a few test cases it wasn't correct.


I haven't really understood the mathematics that are used in neither of the methods.

@fun2code Could you explain a bit further how does this work, or link me to a webpage that explain this? Thanks.
closed account (D80DSL3A)
You did a nice job tackling this in terms of line equations and intercepts. I hadn't noticed your efforts to deal with division by 0.
The vector method removes some of the worry about division by 0. Was my brief justification useful? The vector cross product arises naturally in the slope comparison. The method makes it easy to repeat at a 2nd vertex, completing the test.

Vector algebra can be very useful. I don't have a link handy as I'm applying what I learned in school. My degree is in Physics, and there are lots of vectors in Physics!

Another quick example of using vectors to avoid exceptional cases.
Find the angle between 2 line segments (or vectors) A(Ax,Ay), B(Bx,By).
Algebraically we might find the angle each makes with the x-axis and subtract:
A-x angle = tan-1(Ay/Ax)
B-x angle = tan-1(By/Bx)
then angle A-B = a-x angle - B-x angle.
Division by 0 and other hazards apply.
Introducing the vector scalar (or dot) product:
A*B = |A|*|B|*cos( angle A-B ) = Ax*Bx + Ay*By
where |A| = length of A = sqrt(A*A)
Then:
angle A-B = cos-1(A*B/|A|*|B|)
No division by 0 cases to worry about.

edit: I'm actually underselling the vector based methods.
Vectors really get useful in 3D space.
Find the angle between 2 line segments in x,y,z space?
I think it would be somewhat harder to do algebraically than in the 2D case.
With vectors though it's still just
angle A-B = cos-1(A*B/|A|*|B|)
The dot product gets 1 more term is all:
A*B = Ax*Bx + Ay*By + Az*Bz
Last edited on
Topic archived. No new replies allowed.