diff --git a/masque/file/gdsii.py b/masque/file/gdsii.py
index b9cfd05..cf1eab5 100644
--- a/masque/file/gdsii.py
+++ b/masque/file/gdsii.py
@@ -284,49 +284,54 @@ def read(filename: str,
                 if element.width == 0:
                     continue
 
-                #TODO add extension
                 v = remove_colinear_vertices(numpy.array(element.xy, dtype=float), closed_path=False)
                 dv = numpy.diff(v, axis=0)
                 dvdir = dv / numpy.sqrt((dv * dv).sum(axis=1))[:, None]
+                perp = dvdir[:, ::-1] * [[1, -1]] * element.width / 2
+
+                # add extension
                 v[0] -= dvdir[0] * extension
                 v[-1] += dvdir[-1] * extension
-                perp = dvdir[:, ::-1] * [1, -1] * element.width / 2
+
+                # Find intersections of expanded sides
+                As = numpy.stack((dv[:-1], -dv[1:]), axis=2)
+                bs = v[1:-1] - v[:-2] + perp[1:] - perp[:-1]
+                ds = v[1:-1] - v[:-2] - perp[1:] + perp[:-1]
+
+                rp = numpy.linalg.solve(As, bs)[:, 0, None]
+                rn = numpy.linalg.solve(As, ds)[:, 0, None]
+
+                print(rp)
+
+                intersection_p = v[:-2] + rp * dv[:-1] + perp[:-1]
+                intersection_n = v[:-2] + rn * dv[:-1] - perp[:-1]
+
+                towards_perp = (dv[1:] * perp[:-1]).sum(axis=1) > 0 # path bends towards previous perp?
+#               straight = (dv[1:] * perp[:-1]).sum(axis=1) == 0    # path is straight
+                acute = (dv[1:] * dv[:-1]).sum(axis=1) < 0          # angle is acute?
+
 
                 o0 = [v[0] + perp[0]]
                 o1 = [v[0] - perp[0]]
-
-                for i in range(1, dv.shape[0]):
-                    towards_perp = numpy.dot(perp[i - 1], dv[i]) > 0 # bends towards previous perp
-                    #straight = numpy.dot(perp[i - 1], dv[i]) == 0   # TODO maybe run without cleaning?
-                    acute = numpy.dot(dv[i - 1], dv[i]) < 0
-
-                    A = numpy.column_stack((dv[i - 1], -dv[i]))
-                    ab = numpy.linalg.solve(A, v[i] + perp[i] - v[i - 1] - perp[i - 1])
-                    cd = numpy.linalg.solve(A, v[i] - perp[i] - v[i - 1] + perp[i - 1])
-                    perpside_intersection = v[i - 1] + ab[0] * dv[i - 1] + perp[i - 1]
-                    otherside_intersection = v[i - 1] + cd[0] * dv[i - 1] - perp[i - 1]
-                    if towards_perp:
-                        o0.append(perpside_intersection)
-                        if acute:
-                            # Opposite is 180 > angle > 270
-                            o1.append(otherside_intersection)
+                for i in range(dv.shape[0] - 1):
+                    if towards_perp[i]:
+                        o0.append(intersection_p[i])
+                        if acute[i]:
+                            o1.append(intersection_n[i])
                         else:
                             # Opposite is >270
-                            pt0 = v[i] - perp[i - 1] + dvdir[i - 1] * element.width / 2
-                            pt1 = v[i] - perp[i] - dvdir[i] * element.width / 2
+                            pt0 = v[i + 1] - perp[i + 0] + dvdir[i + 0] * element.width / 2
+                            pt1 = v[i + 1] - perp[i + 1] - dvdir[i + 1] * element.width / 2
                             o1 += [pt0, pt1]
                     else:
-                        # > 180, opposite is <180
-                        o1.append(otherside_intersection)
-                        print('oi', otherside_intersection)
-                        if acute:
-                            # > 270, opposite is <90
-                            pt0 = v[i] + perp[i - 1] + dvdir[i - 1] * element.width / 2
-                            pt1 = v[i] + perp[i] - dvdir[i] * element.width / 2
+                        o1.append(intersection_n[i])
+                        if acute[i]:
+                            # > 270
+                            pt0 = v[i + 1] + perp[i + 0] + dvdir[i + 0] * element.width / 2
+                            pt1 = v[i + 1] + perp[i + 1] - dvdir[i + 1] * element.width / 2
                             o0 += [pt0, pt1]
                         else:
-                            # 180 > angle >270
-                            o0.append(perpside_intersection)
+                            o0.append(intersection_p[i])
                 o0.append(v[-1] + perp[-1])
                 o1.append(v[-1] - perp[-1])
                 verts = numpy.vstack((o0, o1[::-1]))