diff --git a/nijigp/operator_io_paste.py b/nijigp/operator_io_paste.py
index 62c54d8..3410d27 100644
--- a/nijigp/operator_io_paste.py
+++ b/nijigp/operator_io_paste.py
@@ -79,20 +79,6 @@ def execute(self, context):
             context.object.data.materials.append(holdout_material)
             holdout_material_index = len(context.object.data.materials) -1
 
-            # Emitting rays vertically to get winding numbers
-            def pos_H_intersect(point, seg):
-                if seg[0][0]>point[0] and seg[1][0]>point[0]:
-                    return False
-                if seg[0][0]<point[0] and seg[1][0]<point[0]:
-                    return False
-                if math.isclose(seg[0][0], point[0]):
-                    return False
-                if math.isclose(seg[1][0], point[0]):
-                    return seg[1][1] > point[1]
-                ratio = (point[0] - seg[0][0]) / (seg[1][0] - seg[0][0])
-                h_intersect = seg[0][1] + ratio * (seg[1][1] - seg[0][1])
-                return h_intersect > point[1]
-
             # Detect holes using even-odd rule (roughly)
             for i in range(len(context.object.data.layers) - num_layers):
                 strokes = context.object.data.layers[i].active_frame.strokes
@@ -104,28 +90,26 @@ def pos_H_intersect(point, seg):
                 is_hole = False
 
                 while len(to_process)>0 and len(outer_shapes)>0:
-                    winding_number_list = []
+                    crossing_number_list = []
                     is_hole = (not is_hole)
 
                     # Judge whether a stroke is inside another one
                     # Roughly computed by sampling only one point on each curve to avoid too much calculation
                     for stroke_src in to_process:
                         sample_point = vec3_to_vec2(stroke_src.points[len(stroke_src.points)//2].co)
-                        winding_number = 0
+                        crossing_number = 0
                         for stroke_dst in outer_shapes:
                             if stroke_dst != stroke_src:
                                 poly_list, _ = stroke_to_poly([stroke_dst])
                                 co_list = poly_list[0]
                                 for j, co in enumerate(co_list):
-                                    if j!= 0:
-                                        seg = [co_list[j], co_list[j-1]]
-                                        winding_number += pos_H_intersect(sample_point, seg)       
-                        winding_number_list.append(winding_number)  
+                                    crossing_number += raycast_2d_up(sample_point[0], sample_point[1], co_list[j][0], co_list[j][1], co_list[j-1][0], co_list[j-1][1])
+                        crossing_number_list.append(crossing_number)  
 
                     # Process the inner strokes
                     outer_shapes = []
-                    for j,winding_number in enumerate(winding_number_list):
-                        if winding_number % 2 == 1:
+                    for j,crossing_number in enumerate(crossing_number_list):
+                        if crossing_number % 2 == 1:
                             if is_hole:
                                 to_process[j].material_index = holdout_material_index
                             to_process[j].select = True
diff --git a/nijigp/utils.py b/nijigp/utils.py
index cc2192f..5791aee 100644
--- a/nijigp/utils.py
+++ b/nijigp/utils.py
@@ -138,7 +138,30 @@ def intersecting_segments(x1,y1,x2,y2,x3,y3,x4,y4):
         if is_2d_point_on_segment(x1,y1,x2,y2,x_inter,y_inter) and is_2d_point_on_segment(x3,y3,x4,y4,x_inter,y_inter):
             return True
         return False
+
+def raycast_2d_up(x0, y0, x1, y1, x2, y2):
+    """
+    Emit ray from (x0, y0) in +y direction and check if it intersects with the line segment (x1,y1)->(x2,y2)
+    mathutils.geometry.intersect_line_line_2d sometimes has inconsistent results, therefore is not adopted here
+    """
+    if x1>x0 and x2>x0:
+        return False
+    if x1<x0 and x2<x0:
+        return False
     
+    # Several special cases to consider
+    if math.isclose(x1,x2):
+        return False
+    if x1<x2 and math.isclose(x1, x0):
+        return y1 > y0
+    if x1>x2 and math.isclose(x2, x0):
+        return y2 > y0
+
+    # The normal case
+    ratio = (x0 - x1) / (x2 - x1)
+    h_intersect = y1 + ratio * (y2 - y1)
+    return h_intersect > y0
+
 def overlapping_strokes(s1, s2):
     """
     Check if two strokes overlap with each other. Ignore the cases involving holes