@@ -30,33 +30,32 @@ public static bool pointAboveOrOnLine(Point pt, Point left, Point right)
3030 var Cx = pt . x ;
3131 var Cy = pt . y ;
3232
33- return ( Bx - Ax ) * ( Cy - Ay ) - ( By - Ay ) * ( Cx - Ax ) >= - eps ;
33+ var ABx = Bx - Ax ;
34+ var ABy = By - Ay ;
35+ var AB = Math . Sqrt ( ABx * ABx + ABy * ABy ) ;
36+ // algebraic distance of 'pt' to ('left', 'right') line is:
37+ // [ABx * (Cy - Ay) - ABy * (Cx - Ax)] / AB
38+ return ABx * ( Cy - Ay ) - ABy * ( Cx - Ax ) >= - eps * AB ;
3439 }
3540
3641 public static bool pointBetween ( Point pt , Point left , Point right )
3742 {
3843 // p must be collinear with left->right
3944 // returns false if p == left, p == right, or left == right
45+ if ( pointsSame ( pt , left ) || pointsSame ( pt , right ) ) return false ;
4046 var d_py_ly = pt . y - left . y ;
4147 var d_rx_lx = right . x - left . x ;
4248 var d_px_lx = pt . x - left . x ;
4349 var d_ry_ly = right . y - left . y ;
4450
4551 var dot = d_px_lx * d_rx_lx + d_py_ly * d_ry_ly ;
4652
47- // if `dot` is 0, then `p` == `left` or `left` == `right` (reject)
48- // if `dot` is less than 0, then `p` is to the left of `left` (reject)
49- if ( dot < eps )
50- return false ;
51-
53+ // dot < 0 is p is to the left of 'left'
54+ if ( dot < 0 ) return false ;
5255 var sqlen = d_rx_lx * d_rx_lx + d_ry_ly * d_ry_ly ;
5356
54- // if `dot` > `sqlen`, then `p` is to the right of `right` (reject)
55- // therefore, if `dot - sqlen` is greater than 0, then `p` is to the right of `right` (reject)
56- if ( dot - sqlen > - eps )
57- return false ;
58-
59- return true ;
57+ // dot <= sqlen is p is to the left of 'right'
58+ return dot <= sqlen ;
6059 }
6160
6261 public static bool pointsSameX ( Point p1 , Point p2 )
@@ -95,7 +94,18 @@ public static bool pointsCollinear(Point p1, Point p2, Point p3)
9594 var dx2 = p2 . x - p3 . x ;
9695 var dy2 = p2 . y - p3 . y ;
9796
98- return Math . Abs ( dx1 * dy2 - dx2 * dy1 ) < eps ;
97+ var n1 = Math . Sqrt ( dx1 * dx1 + dy1 * dy1 ) ;
98+ var n2 = Math . Sqrt ( dx2 * dx2 + dy2 * dy2 ) ;
99+ // Assuming det(u, v) = 0, we have:
100+ // |det(u + u_err, v + v_err)| = |det(u + u_err, v + v_err) - det(u,v)|
101+ // =|det(u, v_err) + det(u_err. v) + det(u_err, v_err)|
102+ // <= |det(u, v_err)| + |det(u_err, v)| + |det(u_err, v_err)|
103+ // <= N(u)N(v_err) + N(u_err)N(v) + N(u_err)N(v_err)
104+ // <= eps * (N(u) + N(v) + eps)
105+ // We have N(u) ~ N(u + u_err) and N(v) ~ N(v + v_err).
106+ // Assuming eps << N(u) and eps << N(v), we end with:
107+ // |det(u + u_err, v + v_err)| <= eps * (N(u + u_err) + N(v + v_err))
108+ return Math . Abs ( dx1 * dy2 - dx2 * dy1 ) <= eps * ( n1 + n2 ) ;
99109 }
100110
101111 public static bool linesIntersect ( Point a0 , Point a1 , Point b0 , Point b1 , out Intersection intersection )
@@ -125,7 +135,9 @@ public static bool linesIntersect(Point a0, Point a1, Point b0, Point b1, out In
125135 var bdy = b1 . y - b0 . y ;
126136
127137 var axb = adx * bdy - ady * bdx ;
128- if ( Math . Abs ( axb ) < eps )
138+ var n1 = Math . Sqrt ( adx * adx + ady * ady ) ;
139+ var n2 = Math . Sqrt ( bdx * bdx + bdy * bdy ) ;
140+ if ( Math . Abs ( axb ) <= eps * ( n1 + n2 ) )
129141 {
130142 intersection = Intersection . Empty ;
131143 return false ; // lines are coincident
@@ -137,39 +149,37 @@ public static bool linesIntersect(Point a0, Point a1, Point b0, Point b1, out In
137149 var A = ( bdx * dy - bdy * dx ) / axb ;
138150 var B = ( adx * dy - ady * dx ) / axb ;
139151
152+ var pt = new Point ( )
153+ {
154+ x = a0 . x + A * adx ,
155+ y = a0 . y + A * ady
156+ } ;
157+
140158 intersection = new Intersection ( )
141159 {
142160 alongA = 0 ,
143161 alongB = 0 ,
144- pt = new Point ( )
145- {
146- x = a0 . x + A * adx ,
147- y = a0 . y + A * ady
148- }
162+ pt = pt
149163 } ;
150164
151165 // categorize where intersection point is along A and B
152166
153- if ( A <= - eps )
154- intersection . alongA = - 2 ;
155- else if ( A < eps )
167+ if ( pointsSame ( pt , a0 ) )
156168 intersection . alongA = - 1 ;
157- else if ( A - 1 <= - eps )
158- intersection . alongA = 0 ;
159- else if ( A - 1 < eps )
169+ else if ( pointsSame ( pt , a1 ) )
160170 intersection . alongA = 1 ;
161- else
171+ else if ( A < 0 )
172+ intersection . alongA = - 2 ;
173+ else if ( A > 1 )
162174 intersection . alongA = 2 ;
163175
164- if ( B <= - eps )
165- intersection . alongB = - 2 ;
166- else if ( B < eps )
176+ if ( pointsSame ( pt , b0 ) )
167177 intersection . alongB = - 1 ;
168- else if ( B - 1 <= - eps )
169- intersection . alongB = 0 ;
170- else if ( B - 1 < eps )
178+ else if ( pointsSame ( pt , b1 ) )
171179 intersection . alongB = 1 ;
172- else
180+ else if ( B < 0 )
181+ intersection . alongB = - 2 ;
182+ else if ( B > 1 )
173183 intersection . alongB = 2 ;
174184
175185 return true ;
0 commit comments