diff --git a/Makefile b/Makefile
index c786681..e725bfa 100644
--- a/Makefile
+++ b/Makefile
@@ -16,14 +16,14 @@
 #	gcc -Wall -Wextra -Werror -pedantic -g -o main.o main.c -c
 
 # Link command:
-voronoi2: common.o geometry.o dcel.o voronoi.o input.o main.o
-	gcc -Wall -Wextra -Werror -pedantic -g -o voronoi2 main.o input.o voronoi.o dcel.o geometry.o common.o
+voronoi2: common.o towers.o dcel.o voronoi.o input.o main.o
+	gcc -Wall -Wextra -Werror -pedantic -g -o voronoi2 main.o input.o voronoi.o dcel.o towers.o common.o -lm
 
 common.o: common.c
 	gcc -Wall -Wextra -Werror -pedantic -g -o common.o common.c -c
 
-geometry.o: geometry.c
-	gcc -Wall -Wextra -Werror -pedantic -g -o geometry.o geometry.c -c
+towers.o: towers.c
+	gcc -Wall -Wextra -Werror -pedantic -g -o towers.o towers.c -c
 
 dcel.o: dcel.c
 	gcc -Wall -Wextra -Werror -pedantic -g -o dcel.o dcel.c -c
diff --git a/common.h b/common.h
index 7709843..8527956 100644
--- a/common.h
+++ b/common.h
@@ -19,6 +19,13 @@
 
 #endif
 
+#ifndef MATH_HEADER
+#define MATH_HEADER
+
+#include <math.h>
+
+#endif
+
 #ifndef COMMON_HEADER
 #define COMMON_HEADER
 
diff --git a/common.o b/common.o
index 734d90f..d255b40 100644
Binary files a/common.o and b/common.o differ
diff --git a/dcel.c b/dcel.c
index 52af33b..9ba0348 100644
--- a/dcel.c
+++ b/dcel.c
@@ -2,22 +2,21 @@
  *
  * Created by Rory Healy (healyr@student.unimelb.edu.au)
  * Created on 25th August 2021
- * Last modified 9th September 2021
+ * Last modified 11th September 2021
  *
- * Contains functions for the DCEL data structure, including creating, 
+ * Contains functions for the DCEL data structure, including initialisation, 
  * splitting an edge, and identifying towers in faces.
  *
+ * Contains several functions and structures from Grady Fitzpatrick's Base
+ * Code on the Ed Discussion Forum.
+ * 
  */
 
 #ifndef DCEL_HEADER
 #include "dcel.h"
 #endif
 
-#ifndef COMMON_HEADER
-#include "common.h"
-#endif
-
-/* Reads the polygon file and stores the information in the vertices array. */
+/* Deprecated by Grady's sample solution */
 vertex_t **readPolygon(vertex_t **vertices, FILE *polygonFile, \
     int *numVertices) {
     double currentX, currentY;
@@ -47,3 +46,1073 @@ void freeVertices(vertex_t **vertices, int numVertices) {
     }
     free(vertices);
 }
+
+bisector_t *getBisector(vertex_t *pointA, vertex_t *pointB) {
+    bisector_t *newBisector = malloc(sizeof(*newBisector));
+    checkNullPointer(newBisector);
+
+    double midpointX = (pointA->x + pointB->x) / 2;
+    double midpointY = (pointA->y + pointB->y) / 2;
+
+    newBisector->midX = midpointX;
+    newBisector->midY = midpointY;
+
+    /* Calculating bisector slope according to slope of AB */
+    if (pointA->x == pointB->x) {
+        /* The line segment AB has an infinite gradient, 
+         * so the orthogonal line will have zero slope.
+         */
+        newBisector->slope = 0;
+        newBisector->isSlopeInfinite = 0;
+    } else if (pointA->y == pointB->y) {
+        /* The line segment AB has gradient of zero, so 
+         * the orthogonal line will have an infinite slope. 
+         */
+        newBisector->isSlopeInfinite = 1;
+
+        /* Not actually zero, just a placeholder to prevent 
+         * accidental errors 
+         */
+        newBisector->slope = 0;
+    } else {
+        /* The line segment AB has a non-zero and finite gradient, 
+         * so the gradient of the bisector can be calculated.
+         */
+        newBisector->isSlopeInfinite = 0;
+        newBisector->slope = -1 / \
+            ((pointB->y - pointA->y) / (pointB->x - pointA->x));
+    }
+
+    return newBisector;
+}
+
+/* Here on out is the base code from Grady, with some alterations from me */
+
+#define INITIALVERTICES 4
+#define INITIALEDGES 4
+#define INITIALFACES 1
+#define NOVERTEX (-1)
+#define NOEDGE (-1)
+
+#define DIR_UNDECIDED (0)
+#define INSIDE (1)
+#define OUTSIDE (-1)
+#define NODIAMETER (-1)
+
+double dist_2D(double x1, double y1, double x2, double y2) {
+    return sqrt(pow(x2 - x1, 2) + pow(y2 - y1, 2));
+}
+
+double inv_dist_2D(double x1, double y1, double m, double c1, double d) {
+    /* Solve sqrt((y - y1)^2 + (x - x1)^2) = d^2 for x by substituting 
+     * y = mx + c1, which gives an equation for x 
+     * Using this, you can find a point that lies on the line y = mx + c1
+     * that is of distance d away from (x1, y1).
+     */
+    double a = m * m + 1;
+    double b = 2 * c1 * m - 2 * x1 - 2 * m * y1;
+    double c = c1 * c1 + x1 * x1 + y1 * y1 - 2 * c1 * y1 - d * d;
+    return (-b + sqrt(b * b - 4 * a * c)) / (2 * a);
+}
+
+vertex_t *getBisectorPoint(double distance, bisector_t *b) {
+    vertex_t *returnPoint = malloc(sizeof(*returnPoint));
+    if (b->isSlopeInfinite) {
+        returnPoint->x = b->midX;
+        returnPoint->y = b->midY + distance;
+    } else if (b->slope == 0) {
+        returnPoint->x = b->midX + distance;
+        returnPoint->y = b->midY;
+    } else {
+        /* Find the y-intercept of the bisector, use it to find the general 
+         * form of the bisector in the form of y = mx + c
+         */
+        double c = b->midY - b->slope * b->midX;
+        returnPoint->x = inv_dist_2D(b->midX, b->midY, b->slope, c, distance);
+        returnPoint->y = b->slope * returnPoint->x + c;
+    }
+
+    /* Due to floating-point precision, there may be a situation where the 
+     * distance between the bisector midpoint and returnPoint is slightly less 
+     * than distance, but unless getBisectorPoint needs to guarentee that this 
+     * distance is greater than distance, it shouldn't be an issue.
+     */
+
+    return returnPoint;
+}
+
+/* 
+This intersection is based on code by Joseph O'Rourke and is provided for use in 
+COMP20003 Assignment 2.
+
+The approach for intersections is:
+- Use the bisector to construct a finite segment and test it against the half-edge.
+- Use O'Rourke's segseg intersection 
+    (https://hydra.smith.edu/~jorourke/books/ftp.html) 
+  to check if the values overlap/intersect
+
+Generates a segment with each end at least minLength away in each direction 
+from the bisector midpoint.
+Returns 1 if b intersects the given half-edge on this segment, 0 otherwise.
+Sets the intersection point to the given x, y positions.
+*/
+
+int areaSign(double sx, double sy, double ex, double ey, double x, double y) {
+    double areaSq;
+    /* |AB x AC|^2, squared area */
+    /* See https://mathworld.wolfram.com/CrossProduct.html */
+    areaSq = (ex - sx) * (y  - sy) -
+             (x  - sx) * (ey - sy);
+    
+    if (areaSq > 0.0) {
+        return 1;
+    } else if (areaSq == 0.0) {
+        return 0;
+    } else {
+        return -1;
+    }
+}
+
+int between(double sx, double sy, double ex, double ey, double x, double y) {
+    if (sx != ex) {
+        /* If not vertical, check whether between x. */
+        if ((sx <= x && x <= ex) || (sx >= x && x >= ex)) {
+            return 1;
+        } else {
+            return 0;
+        }
+    } else {
+        /* Vertical, so can't check _between_ x-values. Check y-axis. */
+        if ((sy <= y && y <= ey) || (sy >= y && y >= ey)) {
+            return 1;
+        } else {
+            return 0;
+        }
+    }
+}
+    
+int collinear(double sx, double sy, double ex, double ey, double x, double y) {
+    /* If area of the parallelogram is 0, then the points 
+     * are in the same line. 
+     */
+    if (areaSign(sx, sy, ex, ey, x, y) == 0) {
+        return 1;
+    } else {
+        return 0;
+    }
+}
+
+enum intersectType parallelIntersects(double heSx, double heSy, \
+    double heEx, double heEy, double bSx, double bSy, \
+    double bEx, double bEy, double *x, double *y) {
+    if (!collinear(heSx, heSy, heEx, heEy, bSx, bSy)) {
+        /* Parallel, no intersection so don't set (x, y) */
+        return DOESNT_INTERSECT;
+    }
+    /* bS between heS and heE */
+    if (between(heSx, heSy, heEx, heEy, bSx, bSy)) {
+        *x = bSx; 
+        *y = bSy;
+        return SAME_LINE_OVERLAP;
+    }
+    /* bE between heS and heE */
+    if (between(heSx, heSy, heEx, heEy, bEx, bEy)) {
+        *x = bEx;
+        *y = bEy;
+        return SAME_LINE_OVERLAP;
+    }
+    /* heS between bS and bE */
+    if (between(bSx, bSy, bEx, bEy, heSx, heSy)) {
+        *x = heSx;
+        *y = heSy;
+        return SAME_LINE_OVERLAP;
+    }
+    /* heE between bS and bE */
+    if (between(bSx, bSy, bEx, bEy, heEx, heEy)) {
+        *x = heEx; 
+        *y = heEy;
+        return SAME_LINE_OVERLAP;
+    }
+    
+    return DOESNT_INTERSECT;
+}
+
+enum intersectType intersects(halfEdge_t *he, bisector_t *b, \
+    DCEL_t *dcel, double minLength, double *x, double *y) {
+    /* Half-edge x, y twin */
+    double heSx = dcel->vertices[he->startVertex].x;
+    double heSy = dcel->vertices[he->startVertex].y;
+    double heEx = dcel->vertices[he->endVertex].x;
+    double heEy = dcel->vertices[he->endVertex].y;
+    
+    /* Bisector x, y twin */
+    double bSx = b->startX;
+    double bSy = b->startY;
+    double bEx = b->endX;
+    double bEy = b->endY;
+    
+    /* Fill in segment. */
+    // min-length not used here?
+    printf("In intersets(): minlength = %lf\n", minLength);
+    
+    /* Parametric equation parameters */
+    double t1, t2;
+
+    /* Numerators for X and Y coordinate of intersection. */
+    double numeratorX, numeratorY;
+
+    /* Denominators of intersection coordinates. */
+    double denominator;
+    
+    /*
+    See http://www.cs.jhu.edu/~misha/Spring20/15.pdf
+    for explanation and intuition of the algorithm here.
+    x_1 = heSx, y_1 = heSy    |    p_1 = heS
+    x_2 = heEx, y_2 = heEy    |    q_1 = heE
+    x_3 = bSx , y_3 = bSy     |    p_2 =  bS
+    x_4 = bEx , y_4 = bEy     |    q_2 =  bE
+    ----------------------------------------
+    So the parameters t1 and t2 are given by:
+    | t1 |   | heEx - heSx  bSx - bEx | -1  | bSx - heSx |
+    |    | = |                        |     |            |
+    | t2 |   | heEy - heSy  bSy - bEy |     | bSy - heSy |
+    
+    Hence:
+    | t1 |       1     | bSy - bEy        bEx - bSx |  | bSx - heSx |
+    |    | = --------- |                            |  |            |
+    | t2 |    ad - bc  | heSy - heEy    heEx - heSx |  | bSy - heSy |
+    
+        where 
+        a = heEx - heSx
+        b = bSx  -  bEx
+        c = heEy - heSy
+        d = bSy  -  bEy
+    */
+    
+    /* Here we calculate ad - bc */
+    denominator = heSx * (bEy  -  bSy) +
+                  heEx * (bSy  -  bEy) +
+                  bEx  * (heEy - heSy) +
+                  bSx  * (heSy - heEy);
+    
+    if (denominator == 0) {
+        /* In this case the two are parallel */
+        return parallelIntersects(heSx, heSy, heEx, heEy, \
+            bSx, bSy, bEx, bEy, x, y);
+    }
+    
+    /*
+    Here we calculate the top row.
+    | bSy - bEy        bEx - bSx |  | bSx - heSx |
+    |                            |  |            |
+    |                            |  | bSy - heSy |
+    */
+    numeratorX = heSx * (bEy  -  bSy) +
+                 bSx  * (heSy -  bEy) +
+                 bEx  * (bSy  - heSy);
+    
+    /*
+    Here we calculate the bottom row.
+    |                            |  | bSx - heSx |
+    |                            |  |            |
+    | heSy - heEy    heEx - heSx |  | bSy - heSy |
+    */
+    numeratorY = -(heSx * (bSy  -  heEy) +
+                   heEx * (heSy -  bSy) +
+                   bSx  * (heEy  - heSy));
+    
+    /* Use parameters to convert to the intersection point */
+    t1 = numeratorX/denominator;
+    t2 = numeratorY/denominator;
+    *x = heSx + t1 * (heEx - heSx);
+    *y = heSy + t1 * (heEy - heSy);
+    
+    /* Make final decision - if point is on segments, parameter values will be
+    between 0, the start of the line segment, and 1, the end of the line segment.
+    */
+    if (0.0 < t1 && t1 < 1.0 && 0.0 < t2 && t2 < 1.0) {
+        return INTERSECT;
+    } else if (t1 < 0.0 || 1.0 < t1 || t2 < 0.0 || 1.0 < t2) {
+        /* s or t outside of line segment. */
+        return DOESNT_INTERSECT;
+    } else {
+        /* 
+        ((numeratorX == 0) || (numeratorY == 0) || 
+         (numeratorX == denominator) || (numeratorY == denominator))
+        */
+        return ENDS_OVERLAP;
+    }
+}
+
+char *getIntersectionString(intersection_t *intersection) {
+    /*
+    
+    
+    
+    FILL IN
+
+
+
+    */
+    if (!intersection) {
+      return NULL;
+    }
+    char *returnString = NULL;
+    if (0 <= 0) {
+        /* Find out memory needed. */
+        int stringLength = snprintf(returnString, 0, 
+            "From Edge %d (%lf, %lf) to Edge %d (%lf, %lf)", 
+            0, 0.0, 0.0,
+            0, 0.0, 0.0);
+        returnString = malloc(sizeof(*returnString) * (stringLength + 1));
+        checkNullPointer(returnString);
+        sprintf(returnString, "From Edge %d (%lf, %lf) to Edge %d (%lf, %lf)", 
+            0, 0.0, 0.0,
+            0, 0.0, 0.0);
+    } else {
+        /* Find out memory needed. */
+        int stringLength = snprintf(returnString, 0, 
+            "From Edge %d (%lf, %lf) to Edge %d (%lf, %lf)", 
+            0, 0.0, 0.0,
+            0, 0.0, 0.0);
+        returnString = malloc(sizeof(*returnString) * (stringLength + 1));
+        checkNullPointer(returnString);
+        sprintf(returnString, "From Edge %d (%lf, %lf) to Edge %d (%lf, %lf)", 
+            0, 0.0, 0.0,
+            0, 0.0, 0.0);
+    }
+    return returnString;
+}
+
+void freeIntersection(intersection_t *intersection) {
+    if (!intersection) {
+        return;
+    }
+    free(intersection);
+}
+
+DCEL_t *newDCEL() {
+    /* Setup DCEL. */
+    DCEL_t *dcel = malloc(sizeof(*dcel));
+    checkNullPointer(dcel);
+
+    dcel->edges = NULL;
+    dcel->edgesUsed = 0;
+    dcel->edgesAllocated = 0;
+
+    dcel->faces = NULL;
+    dcel->facesUsed = 0;
+    dcel->facesAllocated = 0;
+
+    dcel->vertices = NULL;
+    dcel->verticesUsed = 0;
+    dcel->verticesAllocated = 0;
+    
+    return dcel;
+}   
+
+halfEdge_t *newHalfEdge() {
+    halfEdge_t *he = malloc(sizeof(*he));
+    checkNullPointer(he);
+    he->startVertex = NOVERTEX;
+    he->endVertex = NOVERTEX;
+    he->next = NULL;
+    he->previous = NULL;
+    he->twin = NULL;
+    he->face = NOFACE;
+    he->edge = NOEDGE;
+    return he;
+}
+
+void ensureSpaceForVertex(DCEL_t *dcel) {
+    if (!(dcel->vertices)) {
+        dcel->vertices = malloc(sizeof(*dcel->vertices) * INITIALVERTICES);
+        checkNullPointer(dcel->vertices);
+        dcel->verticesAllocated = INITIALVERTICES;
+    } else if ((dcel->verticesUsed + 1) > dcel->verticesAllocated) {
+        dcel->vertices = realloc(dcel->vertices, \
+            sizeof(*dcel->vertices) * dcel->verticesAllocated * 2);
+        checkNullPointer(dcel->vertices);
+        dcel->verticesAllocated = dcel->verticesAllocated * 2;
+    }
+}
+
+void ensureSpaceForEdge(DCEL_t *dcel) {
+    if (!(dcel->edges)) {
+        dcel->edges = malloc(sizeof(*dcel->edges) * INITIALEDGES);
+        checkNullPointer(dcel->edges);
+        dcel->edgesAllocated = INITIALEDGES;
+    } else if ((dcel->edgesUsed + 1) > dcel->edgesAllocated) {
+        dcel->edges = realloc(dcel->edges, \
+            sizeof(*dcel->edges) * dcel->edgesAllocated * 2);
+        checkNullPointer(dcel->edges);
+        dcel->edgesAllocated = dcel->edgesAllocated * 2;
+    }
+}
+
+void ensureSpaceForFace(DCEL_t *dcel) {
+    if (!(dcel->faces)) {
+        dcel->faces = malloc(sizeof(*dcel->faces) * INITIALFACES);
+        checkNullPointer(dcel->faces);
+        dcel->facesAllocated = INITIALFACES;
+    } else if ((dcel->facesUsed + 1) > dcel->facesAllocated) {
+        dcel->faces = realloc(dcel->faces, \
+            sizeof(*dcel->faces) * dcel->facesAllocated * 2);
+        checkNullPointer(dcel->faces);
+        dcel->facesAllocated = dcel->facesAllocated * 2;
+    }
+}
+
+void addEdge(DCEL_t *dcel, int startVertex, int endVertex) {
+    ensureSpaceForEdge(dcel);
+    
+    int newEdge = dcel->edgesUsed;
+    
+    halfEdge_t *newHE = newHalfEdge();
+    newHE->startVertex = startVertex;
+    newHE->endVertex = endVertex;
+    // newHE->next = NULL;
+    // newHE->previous = NULL;
+    // newHE->twin = NULL;
+    // newHE->face = NOFACE;
+    newHE->edge = newEdge;
+    
+    (dcel->edges)[newEdge].halfEdge = newHE;
+    
+    dcel->edgesUsed = dcel->edgesUsed + 1;
+}
+
+void addFace(DCEL_t *dcel, halfEdge_t *he) {
+    ensureSpaceForFace(dcel);
+    (dcel->faces)[dcel->facesUsed].start = he;
+    /* Set the face in the half-edges. */
+    he->face = dcel->facesUsed;
+    
+    (dcel->faces)[dcel->facesUsed].tower = NULL;
+    
+    halfEdge_t *current = he->next;
+    while(current != he) {
+        current->face = dcel->facesUsed;
+        current = current->next;
+    }
+    
+    dcel->facesUsed = dcel->facesUsed + 1;
+}
+
+DCEL_t *readPolygonFile(char *polygonfileName) {
+    DCEL_t *dcel = newDCEL();
+    
+    FILE *polygonFile = fopen(polygonfileName, "r");
+    checkNullPointer(polygonFile);
+    double x;
+    double y;
+    
+    int startVertex = NOVERTEX;
+    int endVertex = NOVERTEX;
+    
+    /* Used to finish off the polygon in the first face. */
+    int firstVertex = NOVERTEX;
+    int firstEdge = NOEDGE;
+    
+    while(fscanf(polygonFile, "%lf %lf", &x, &y) == 2) {
+        ensureSpaceForVertex(dcel);
+
+        (dcel->vertices)[dcel->verticesUsed].x = x;
+        (dcel->vertices)[dcel->verticesUsed].y = y;
+        dcel->verticesUsed = dcel->verticesUsed + 1;
+
+        if (startVertex == NOVERTEX) {
+            startVertex = dcel->verticesUsed - 1;
+            firstVertex = startVertex;
+        } else if (endVertex == NOVERTEX) {
+            /* First edge */
+            endVertex = dcel->verticesUsed - 1;
+            firstEdge = dcel->edgesUsed;
+            addEdge(dcel, startVertex, endVertex);
+        } else {
+            /* Start from last vertex. */
+            startVertex = endVertex;
+            endVertex = dcel->verticesUsed - 1;
+            addEdge(dcel, startVertex, endVertex);
+
+            /* Connect last edge added to newest edge */
+            ((dcel->edges)[dcel->edgesUsed - 2].halfEdge)->next = \
+                (dcel->edges)[dcel->edgesUsed - 1].halfEdge;
+            
+            /* Connect newest edge to last edge added */
+            ((dcel->edges)[dcel->edgesUsed - 1].halfEdge)->previous = \
+                (dcel->edges)[dcel->edgesUsed - 2].halfEdge;
+        }
+    }
+    
+    if (firstEdge == NOEDGE) {
+        fputs("Error: Unable to create DCEL structure for polygon.\n", stderr);
+        exit(EXIT_FAILURE);
+    }
+
+    /* Finalise polygon by adding edge back to first vertex. */
+    int finalEdge = dcel->edgesUsed;
+    addEdge(dcel, endVertex, firstVertex);
+
+    /* Connect previous edge to this edge. */
+    ((dcel->edges)[dcel->edgesUsed - 2].halfEdge)->next = \
+        (dcel->edges)[dcel->edgesUsed - 1].halfEdge;
+
+    /* Connect newest edge to last edge added */
+    ((dcel->edges)[dcel->edgesUsed - 1].halfEdge)->previous = \
+        (dcel->edges)[dcel->edgesUsed - 2].halfEdge;
+    
+    /* Connect final edge back to start edge. */
+    ((dcel->edges)[finalEdge].halfEdge)->next = \
+        (dcel->edges)[firstEdge].halfEdge;
+
+    /* Connect start edge back to final edge. */
+    ((dcel->edges)[firstEdge].halfEdge)->previous = \
+        (dcel->edges)[finalEdge].halfEdge;
+    
+    /* Add face to DCEL - could be any edge we constructed, 
+     * so may as well be the first. 
+     */
+    addFace(dcel, (dcel->edges)[firstEdge].halfEdge);
+    if (polygonFile) {
+        fclose(polygonFile);
+    }
+        
+    return dcel;
+}
+
+split_t *readNextSplit(FILE *splitfile) {
+    int firstEdge;
+    int secondEdge;
+    if (fscanf(splitfile, "%d %d", &firstEdge, &secondEdge) != 2) {
+        return NULL;
+    }
+    split_t *split = malloc(sizeof(*split));
+    split->startEdge = firstEdge;
+    split->endEdge = secondEdge;
+    split->verticesSpecified = 0;
+    return split;
+}
+
+void freeSplit(split_t *split) {
+    if (split) {
+        free(split);
+    }
+}
+
+int vertexMatch(vertex_t *v1, vertex_t *v2) {
+    if (v1->x != v2->x) {
+        return 0;
+    }
+    if (v1->y != v2->y) {
+        return 0;
+    }
+    return 1;
+}
+
+void applySplit(split_t *split, DCEL_t *dcel) {
+    int isAdjacent;
+    double midpointX;
+    double midpointY;
+    halfEdge_t *startHE;
+    halfEdge_t *endHE;
+    halfEdge_t *newJoinHE;
+    halfEdge_t *newJoinHEPair;
+    halfEdge_t *newStartHEToMid;
+    halfEdge_t *newStartHEToMidPair;
+    halfEdge_t *newMidHEToEnd;
+    halfEdge_t *newMidHEToEndPair;
+    /* Temporary holders for old twin edges */
+    halfEdge_t *oldStartPairPrev;
+    halfEdge_t *oldEndPairNext;
+    /* Temporary holder for old pairs */
+    halfEdge_t *oldStartPair;
+    halfEdge_t *oldEndPair;
+    
+    int newVertexMidStart;
+    int newVertexMidEnd;
+    /* The vertex representing the end of the original starting edge */
+    int oldVertexStart;
+    /* The vertex representing the start of the original ending edge */
+    int oldVertexEnd;
+    
+    /* Each split creates exactly 3 edges, so we can set up space for these now. */
+    int joinEdge;
+    int newStartEdge;
+    int newEndEdge;
+    
+    ensureSpaceForEdge(dcel);
+    joinEdge = dcel->edgesUsed;
+    dcel->edgesUsed = dcel->edgesUsed + 1;
+    
+    ensureSpaceForEdge(dcel);
+    newStartEdge = dcel->edgesUsed;
+    dcel->edgesUsed = dcel->edgesUsed + 1;
+    
+    ensureSpaceForEdge(dcel);
+    newEndEdge = dcel->edgesUsed;
+    dcel->edgesUsed = dcel->edgesUsed + 1;
+    
+    /* Get vertices for MidStart and MidEnd */
+    ensureSpaceForVertex(dcel);
+    newVertexMidStart = dcel->verticesUsed;
+    dcel->verticesUsed = dcel->verticesUsed + 1;
+    
+    ensureSpaceForVertex(dcel);
+    newVertexMidEnd = dcel->verticesUsed;
+    dcel->verticesUsed = dcel->verticesUsed + 1;
+
+    /* Work out what half-edges we need to use. */
+    startHE = (dcel->edges)[split->startEdge].halfEdge;
+    endHE = (dcel->edges)[split->endEdge].halfEdge;
+    
+    /* Set midpoint of start */
+    double startX = (dcel->vertices)[startHE->startVertex].x;
+    double startY = (dcel->vertices)[startHE->startVertex].y;
+    double endX = (dcel->vertices)[startHE->endVertex].x;
+    double endY = (dcel->vertices)[startHE->endVertex].y;
+    if (split->verticesSpecified) {
+        /* See if vertex needs to be reused */
+        if (vertexMatch(&(dcel->vertices)[startHE->endVertex], 
+                       &split->startPoint)) {
+            newVertexMidStart = startHE->endVertex;
+        } else if (vertexMatch(&(dcel->vertices)[startHE->startVertex], 
+                              &split->startPoint)) {
+            newVertexMidStart = startHE->startVertex;
+        } else {
+            (dcel->vertices)[newVertexMidStart].x = split->startPoint.x;
+            (dcel->vertices)[newVertexMidStart].y = split->startPoint.y;
+        }
+    } else {
+        (dcel->vertices)[newVertexMidStart].x = (startX + endX) / 2.0;
+        (dcel->vertices)[newVertexMidStart].y = (startY + endY) / 2.0;
+    }
+    
+    
+    /* Set midpoint of end */
+    startX = (dcel->vertices)[endHE->startVertex].x;
+    startY = (dcel->vertices)[endHE->startVertex].y;
+    endX = (dcel->vertices)[endHE->endVertex].x;
+    endY = (dcel->vertices)[endHE->endVertex].y;
+    if (split->verticesSpecified) {
+        /* See if vertex needs to be reused */
+        if (vertexMatch(&(dcel->vertices)[endHE->startVertex], 
+                       &split->endPoint)) {
+            newVertexMidEnd = endHE->startVertex;
+        } else if (vertexMatch(&(dcel->vertices)[endHE->endVertex], 
+                              &split->endPoint)) {
+            newVertexMidEnd = endHE->endVertex;
+        } else {
+            (dcel->vertices)[newVertexMidEnd].x = split->endPoint.x;
+            (dcel->vertices)[newVertexMidEnd].y = split->endPoint.y;
+        }
+    } else {
+        (dcel->vertices)[newVertexMidEnd].x = (startX + endX) / 2.0;
+        (dcel->vertices)[newVertexMidEnd].y = (startY + endY) / 2.0;
+    }
+    
+    
+    /* Get point halfway between both midpoints */
+    double x1 = (dcel->vertices)[newVertexMidStart].x;
+    double x2 = (dcel->vertices)[newVertexMidEnd].x;
+    double y1 = (dcel->vertices)[newVertexMidStart].y;
+    double y2 = (dcel->vertices)[newVertexMidEnd].y;
+    midpointX = (x1 + x2) / 2.0;
+    midpointY = (y1 + y2) / 2.0;
+    
+    /* Work out whether on correct side. */
+    vertex_t *v1 = &((dcel->vertices)[startHE->startVertex]);
+    vertex_t *v2 = &((dcel->vertices)[startHE->endVertex]);
+    if (getRelativeDir(midpointX, midpointY, v1, v2) == OUTSIDE) {
+        startHE = startHE->twin;
+    }
+    v1 = &((dcel->vertices)[endHE->startVertex]);
+    v2 = &((dcel->vertices)[endHE->endVertex]);
+    if (getRelativeDir(midpointX, midpointY, v1, v2) == OUTSIDE) {
+        endHE = endHE->twin;
+    }
+    
+    /* Work out whether edges are adjacent. */
+    if (startHE->next == endHE) {
+        isAdjacent = 1;
+    } else {
+        isAdjacent = 0;
+    }
+    
+    /* Store old previous and next from start and end edges for convenience */
+    halfEdge_t *oldEndPrev = endHE->previous;
+    halfEdge_t *oldStartNext = startHE->next;
+    oldVertexEnd = endHE->startVertex;
+    oldVertexStart = startHE->endVertex;
+    
+    /* Update vertices of endHE and startHE */
+    endHE->startVertex = newVertexMidEnd;
+    startHE->endVertex = newVertexMidStart;
+    
+    /* Add bridging edges */
+    newJoinHE = newHalfEdge();
+    
+    newJoinHE->startVertex = newVertexMidStart;
+    newJoinHE->endVertex = newVertexMidEnd;
+    newJoinHE->next = endHE;
+    endHE->previous = newJoinHE;
+    newJoinHE->previous = startHE;
+    startHE->next = newJoinHE;
+    newJoinHE->twin = NULL; // Will be set later
+    /* joinHE is same face as startHE and endHE */
+    newJoinHE->face = startHE->face;
+    newJoinHE->edge = joinEdge;
+    
+    /* Set joinEdge to relevant halfEdge */
+    (dcel->edges)[joinEdge].halfEdge = newJoinHE;
+    
+    newJoinHEPair = newHalfEdge();
+    /* twin is in opposite direction. */
+    newJoinHEPair->startVertex = newVertexMidEnd;
+    newJoinHEPair->endVertex = newVertexMidStart;
+    newJoinHEPair->next = NULL; // Will join to new HEs
+    newJoinHEPair->previous = NULL; // Will join to new HEs
+    newJoinHEPair->twin = newJoinHE;
+    newJoinHE->twin = newJoinHEPair;
+    newJoinHEPair->face = NOFACE; // Will be new face set later
+    newJoinHEPair->edge = joinEdge;
+    
+    /* Set up what we can of new edges */
+    newStartHEToMid = newHalfEdge();
+    newStartHEToMid->startVertex = newVertexMidStart;
+    newStartHEToMid->endVertex = oldVertexStart;
+    newStartHEToMid->next = NULL; // Different setting based on adjacency, set below.
+    newStartHEToMid->previous = newJoinHEPair;
+    newJoinHEPair->next = newStartHEToMid;
+    newStartHEToMid->twin = NULL; // Will be set up later if needed.
+    newStartHEToMid->face = NOFACE; // Will be new face set later
+    newStartHEToMid->edge = newStartEdge;
+    
+    /* Set newStartEdge to relevant halfEdge */
+    (dcel->edges)[newStartEdge].halfEdge = newStartHEToMid;
+    
+    newMidHEToEnd = newHalfEdge();
+    newMidHEToEnd->startVertex = oldVertexEnd;
+    newMidHEToEnd->endVertex = newVertexMidEnd;
+    newMidHEToEnd->next = newJoinHEPair;
+    newJoinHEPair->previous = newMidHEToEnd;
+    newMidHEToEnd->previous = NULL; // Different setting based on adjacency, set below.
+    newMidHEToEnd->twin = NULL; // Will be set up later if needed.
+    newMidHEToEnd->face = NOFACE;
+    newMidHEToEnd->edge = newEndEdge;
+    
+    /* Set newEndEdge to relevant halfEdge */
+    (dcel->edges)[newEndEdge].halfEdge = newMidHEToEnd;
+    
+    /* If either start or end HEs have paired Half-Edges, we also need to split those. */
+    if (startHE->twin) {
+        oldStartPairPrev = startHE->twin->previous;
+        oldStartPair = startHE->twin;
+        
+        newStartHEToMidPair = newHalfEdge();
+        /* Reverse of twin */
+        newStartHEToMidPair->startVertex = oldVertexStart;
+        newStartHEToMidPair->endVertex = newVertexMidStart;
+        newStartHEToMidPair->next = oldStartPair;
+        newStartHEToMidPair->previous = oldStartPairPrev;
+        startHE->twin->previous = newStartHEToMidPair;
+        oldStartPair->previous = newStartHEToMidPair;
+        oldStartPair->startVertex = newVertexMidStart;
+        oldStartPairPrev->next = newStartHEToMidPair;
+        newStartHEToMid->twin = newStartHEToMidPair;
+        newStartHEToMidPair->twin = newStartHEToMid;
+        newStartHEToMidPair->face = startHE->twin->face;
+        newStartHEToMidPair->edge = newStartEdge;
+    } else {
+        newStartHEToMidPair = NULL;
+    }
+    if (endHE->twin) {
+        oldEndPairNext = endHE->twin->next;
+        oldEndPair = endHE->twin;
+        
+        newMidHEToEndPair = newHalfEdge();
+        newMidHEToEndPair->startVertex = newVertexMidEnd;
+        newMidHEToEndPair->endVertex = oldVertexEnd;
+        newMidHEToEndPair->next = oldEndPairNext; // endHE->twin ?
+        oldEndPair->next = newMidHEToEndPair;
+        oldEndPairNext->previous = newMidHEToEndPair; // Next?
+        oldEndPair->endVertex = newVertexMidEnd;
+        newMidHEToEndPair->previous = oldEndPair;
+        newMidHEToEnd->twin = newMidHEToEndPair;
+        newMidHEToEndPair->twin = newMidHEToEnd;
+        newMidHEToEndPair->face = endHE->twin->face;
+        newMidHEToEndPair->edge = newEndEdge;
+    } else {
+        newMidHEToEndPair = NULL;
+    }
+    
+    /* Set up remaining edges. */
+    if (isAdjacent) {
+        newStartHEToMid->next = newMidHEToEnd;
+        newMidHEToEnd->previous = newStartHEToMid;
+    } else {
+        /* Edges are old start and end edges (maybe the same edge). */
+        newStartHEToMid->next = oldStartNext;
+        oldStartNext->previous = newStartHEToMid;
+        newMidHEToEnd->previous = oldEndPrev;
+        oldEndPrev->next = newMidHEToEnd;
+    }
+    
+    /* Setup new face. */
+    addFace(dcel, newJoinHEPair);
+    
+    /* Check if face has overwritten other face */
+    int joinFace = startHE->face;
+    if ((dcel->faces)[joinFace].start->face != joinFace) {    
+        (dcel->faces)[joinFace].start = startHE;
+    }
+}
+
+void freeDCEL(DCEL_t *dcel) {
+    if (!dcel) {
+        return;
+    }
+    int i;
+    if (dcel->edges) {
+        for(i = 0; i < dcel->edgesUsed; i++) {
+            if ((dcel->edges)[i].halfEdge) {
+                if (((dcel->edges)[i]).halfEdge->twin) {
+                    /* Free if edge has two halves. */
+                    free(((dcel->edges)[i]).halfEdge->twin);
+                }
+                free(((dcel->edges)[i]).halfEdge);
+            }
+        }
+        free(dcel->edges);
+    }
+    if (dcel->faces) {
+        /* All edges are freed above, so no need to free each edge here. */
+        free(dcel->faces);
+    }
+    if (dcel->vertices) {
+        free(dcel->vertices);
+    }
+    free(dcel);
+}
+
+int getFaceCount(DCEL_t *dcel) {
+    if (!dcel) {
+        return 0;
+    } else {
+        return dcel->facesUsed;
+    }
+}
+
+int getRelativeDir(double x, double y, vertex_t *v1, vertex_t *v2) {
+    /* Here we're doing a simple half-plane check against the vector v1->v2. */
+    double x1 = v1->x;
+    double x2 = v2->x;
+    double y1 = v1->y;
+    double y2 = v2->y;
+    if (x1 == x2 && y1 == y2) {
+        /* Same point. */
+        return DIR_UNDECIDED;
+    } else if (x1 == x2) {
+        /* y = c line */
+        /* Work out whether line is going up or down. */
+        if (y2 > y1) {
+            if (x > x1) {
+                return INSIDE;
+            } else if (x < x1) {
+                return OUTSIDE;
+            } else {
+                return DIR_UNDECIDED;
+            }
+        } else {
+            if (x < x1) {
+                return INSIDE;
+            } else if (x > x1) {
+                return OUTSIDE;
+            } else {
+                return DIR_UNDECIDED;
+            }
+        }
+    } else if (y1 == y2) {
+        /* x = c line */
+        /* Work out whether line is going left or right. */
+        if (x2 > x1) {
+            if (y < y1) {
+                return INSIDE;
+            } else if (y > y1) {
+                return OUTSIDE;
+            } else {
+                return DIR_UNDECIDED;
+            }
+        } else {
+            if (y > y1) {
+                return INSIDE;
+            } else if (y < y1) {
+                return OUTSIDE;
+            } else {
+                return DIR_UNDECIDED;
+            }
+        }
+    }
+    
+    /* 
+        x1, x2, y1, y2 distinct, so see whether point being tested is
+        above or below gradient line.
+    */
+    double m = (y2 - y1)/(x2 - x1);
+    double c = y1 - m*x1;
+    
+    double predictedY = x * m + c;
+    double residual = y - predictedY;
+    
+    /*
+        Being inside or outside the polygon depends on the direction
+        the half-edge is going.
+    */
+    if (x2 > x1) {
+        if (residual < 0) {
+            return INSIDE;
+        } else if (residual > 0) {
+            return OUTSIDE;
+        } else {
+            return DIR_UNDECIDED;
+        }
+    } else {
+        if (residual > 0) {
+            return INSIDE;
+        } else if (residual < 0) {
+            return OUTSIDE;
+        } else {
+            return DIR_UNDECIDED;
+        }
+    }
+}
+
+int directionOrUndecided(int decidedDirection, int direction) {
+    if (direction == decidedDirection || direction == DIR_UNDECIDED) {
+        return 1;
+    } else {
+        return 0;
+    }
+}
+
+int inFace(DCEL_t *dcel, double x, double y, int faceIndex) {
+    if (dcel->facesUsed < faceIndex || !(dcel->faces)[faceIndex].start) {
+        return OUTSIDE;
+    }
+    halfEdge_t *start = (dcel->faces)[faceIndex].start;
+    int first = 1;
+    int direction = DIR_UNDECIDED;
+    
+    halfEdge_t *current = start;
+    while(start != current || first) {
+        if (direction == DIR_UNDECIDED) {
+            /* Doesn't matter where the point is until we find it on one side or the 
+            other. */
+            direction = getRelativeDir(x, y, &(dcel->vertices)[current->startVertex],
+                &(dcel->vertices)[current->endVertex]);
+        } else {
+            if (!directionOrUndecided(direction, 
+                getRelativeDir(x, y, &(dcel->vertices)[current->startVertex],
+                    &(dcel->vertices)[current->endVertex]))) {
+                /* If the point is on the different side of any edge, it be inside 
+                    the face, because the face is convex. */
+                return 0;
+            }
+        }
+        current = current->next;
+        first = 0;
+    }
+    
+    return 1;
+}
+
+int getDCELPointCount(DCEL_t *dcel) {
+    if (!dcel) {
+        return 0;
+    }
+    return dcel->verticesUsed;
+}
+
+double getDCELVertexX(DCEL_t *dcel, int vertex) {
+    return (dcel->vertices)[vertex].x;
+}
+
+double getDCELVertexY(DCEL_t *dcel, int vertex) {
+    return (dcel->vertices)[vertex].y;
+}
+
+int getDCELEdgeCount(DCEL_t *dcel) {
+    if (!dcel) {
+        return 0;
+    }
+    return dcel->edgesUsed;
+}
+
+int getDCELEdgeVertexStart(DCEL_t *dcel, int edge) {
+    if (!dcel) {
+        return 0;
+    }
+    return (dcel->edges)[edge].halfEdge->startVertex;
+}
+
+int getDCELEdgeVertexEnd(DCEL_t *dcel, int edge) {
+    if (!dcel) {
+        return 0;
+    }
+    return (dcel->edges)[edge].halfEdge->endVertex;
+}
+
+int getDCELEdgeVertexPairStart(DCEL_t *dcel, int edge) {
+    if (!dcel) {
+        return 0;
+    }
+    return (dcel->edges)[edge].halfEdge->twin->startVertex;
+}
+
+int getDCELEdgeVertexPairEnd(DCEL_t *dcel, int edge) {
+    if (!dcel) {
+        return 0;
+    }
+    return (dcel->edges)[edge].halfEdge->twin->endVertex;
+}
+
+int DCELhasEdge(DCEL_t *dcel, int edge) {
+    if ((dcel->edges)[edge].halfEdge->face != NOFACE) {
+        return 1;
+    } else {
+        return 0;
+    }
+}
+
+int DCELhasEdgePair(DCEL_t *dcel, int edge) {
+    if ((dcel->edges)[edge].halfEdge->face == NOFACE) {
+        return 0;
+    }
+    if ((dcel->edges)[edge].halfEdge->twin) {
+        return 1;
+    } else {
+        return 0;
+    }
+}
+
+intersection_t *getIntersection(bisector_t *b, DCEL_t *dcel, int face,
+    double minLength) {
+    // FILL IN
+    printf("In getIntersection(): bisector start: %lf, dcel numfaces = %d" \
+        "face = %d, minLength = %lf\n", b->startX, dcel->facesAllocated, \
+        face, minLength);
+    return NULL;
+}
+
+double getDiameter(DCEL_t *dcel, int faceIndex) {
+    // FILL IN
+    printf("in getDiamter(): faceIndex = %d, dcel faces = %d\n", \
+        faceIndex, dcel->facesAllocated);
+    return NODIAMETER;
+}
+
+void incrementalVoronoi(DCEL_t *dcel, tower_t *tower) {
+    // FILL IN
+    printf("In incrementalVoronoi(): tower id = %s, dcel faces = %d\n", \
+        tower->id, dcel->facesAllocated);
+}
diff --git a/dcel.h b/dcel.h
index cabca01..d13827a 100644
--- a/dcel.h
+++ b/dcel.h
@@ -1,32 +1,261 @@
-#ifndef GEOMETRY_HEADER
-#include "geometry.h"
-#endif
-
-#ifndef STDIO_HEADER
-#include <stdio.h>
+#ifndef TOWERS_HEADER
+#include "towers.h"
 #endif
 
 #ifndef DCEL_HEADER
 #define DCEL_HEADER
 
-typedef struct halfEdge {
-    struct halfEdge *previous;
-    struct halfEdge *next;
-    struct halfEdge *twin;
+/* Here is the Base Code from Grady, with some alterations by me */
+
+enum intersectType {
+    DOESNT_INTERSECT  = 0, // Doesn't intersect
+    INTERSECT         = 1, // Intersects
+    SAME_LINE_OVERLAP = 2, // Lines are the same
+    ENDS_OVERLAP      = 3  // Intersects at exactly one point (endpoint)
+};
+
+typedef struct vertex {
+    double x;
+    double y;
+} vertex_t;
+
+typedef struct bisector {
+    double startX;
+    double startY;
+    double endX;
+    double endY;
+
+    /* This refers to the midpoint of the line segment AB, not the midpoint
+     * of the bisector itself.
+     */
+    double midX;
+    double midY;
+
+    int isSlopeInfinite;
+    double slope;
+} bisector_t;
+
+typedef struct halfEdgeLabel {
+    struct halfEdgeLabel *previous;
+    struct halfEdgeLabel *next;
+    struct halfEdgeLabel *twin;
     int face;
     int edge;
+    int startVertex;
+    int endVertex;
 } halfEdge_t;
 
 typedef struct edge {
-    halfEdge_t halfEdge;
+    halfEdge_t *halfEdge;
 } edge_t;
 
 typedef struct face {
-    halfEdge_t start;
+    halfEdge_t *start;
+    tower_t *tower;
 } face_t;
 
+typedef struct split {
+    int startEdge;
+    int endEdge;
+    int verticesSpecified;
+    vertex_t startPoint;
+    vertex_t endPoint;
+} split_t;
+
+typedef struct intersection {
+    vertex_t intersectionPoint;
+} intersection_t;
+
+typedef struct DCEL {
+    edge_t *edges;
+    int edgesUsed;
+    int edgesAllocated;
+
+    face_t *faces;
+    int facesUsed;
+    int facesAllocated;
+
+    vertex_t *vertices;
+    int verticesUsed;
+    int verticesAllocated;
+} DCEL_t;
+
+/* Allocate a new DCEL and return it. */
+DCEL_t *newDCEL();
+
+/* Allocate a new halfEdge and return it. */
+halfEdge_t *newHalfEdge();
+
+/* Returns INSIDE if the points is on the INSIDE of the vector twin by the CW 
+ * winding order, OUTSIDE if it is OUTSIDE by the CW winding order, and 
+ * DIR_UNDECIDED if the point lies on the vector between the points v1 and v2. 
+ */
+int getRelativeDir(double x, double y, vertex_t *v1, vertex_t *v2);
+
+/* Takes an established direction and a new direction, and returns 1 if the 
+ * direction matches the decidedDirection or if the direction is undecided.
+ */
+int directionOrUndecided(int decidedDirection, int direction);
+
+/* Check there's space for another vertex in the DCEL, 
+ * or increase the allocated space. 
+ */
+void ensureSpaceForVertex(DCEL_t *dcel);
+
+/* Check there's space for another edge in the DCEL, 
+ * or increase the allocated space. 
+ */
+void ensureSpaceForEdge(DCEL_t *dcel);
+
+/* Check there's space for another face in the DCEL, 
+ * or increase the allocated space. 
+ */
+void ensureSpaceForFace(DCEL_t *dcel);
+
+/* Add an edge from the startVertex index vertex to the endVertex index. 
+ * Only fills one half-edge as other half-edges will always be added 
+ * through geometry construction.
+ */
+void addEdge(DCEL_t *dcel, int startVertex, int endVertex);
+
+/* Add a face to the DCEL given using the given halfEdge and sets the face. */
+void addFace(DCEL_t *dcel, halfEdge_t *he);
+
+/* Reads the polygon from the given file. */
+DCEL_t *readPolygonFile(char *polygonfileName);
+
+/* Reads the next split from the given file. */
+split_t *readNextSplit(FILE *splitfile);
+
+/* Frees a given split.  */
+void freeSplit(split_t *split);
+
+/* Returns 1 if vertices are sufficiently close, 0 otherwise. */
+int vertexMatch(vertex_t *v1, vertex_t *v2);
+
+/* Gets the string for the given bisector equation. */
+char *getBisectorEquation(bisector_t *b);
+
+/* Frees the given bisector. */
+void freeBisector(bisector_t *bisector);
+
+/* Representation of no face */
+#define NOFACE (-1)
+
+/* Default face for intersections. */
+#define DEFAULT_FACE 0
+
+/* Default minimum length for bisector in each direction */
+#define DEFAULTMINLENGTH (200)
+
+/* Gets the intersection between the given bisector and the given DCEL
+ * for the given face. 
+ */
+intersection_t *getIntersection(bisector_t *b, DCEL_t *dcel, \
+    int face, double minLength);
+
+/* Gets the string for the given intersection. */
+char *getIntersectionString(intersection_t *intersection);
+
+/* Frees a given intersection. */
+void freeIntersection(intersection_t *intersection);
+
+/* Applies a given split to the DCEL. */
+void applySplit(split_t *split, DCEL_t *dcel);
+
+/* Frees the given DCEL */
+void freeDCEL(DCEL_t *dcel);
+
+/* Gets the number of faces in the DCEL. */
+int getFaceCount(DCEL_t *dcel);
+
+/* Returns 1 if the given x,y point is inside the given face. */
+int inFace(DCEL_t *dcel, double x, double y, int faceIndex);
+
+/* Gets the diameter of the given face. */
+double getDiameter(DCEL_t *dcel, int faceIndex);
+
+/* Adds the watchtower to the Voronoi diagram represented by the given DCEL, 
+ * applying required splits and setting the watchtower as required.
+ */
+void incrementalVoronoi(DCEL_t *dcel, tower_t *watchTower);
+
+/* Returns the number of vertices in the DCEL. */
+int getDCELPointCount(DCEL_t *dcel);
+
+/* Get x value of given vertex in DCEL. */
+double getDCELVertexX(DCEL_t *dcel, int vertex);
+
+/* Get y value of given vertex in DCEL. */
+double getDCELVertexY(DCEL_t *dcel, int vertex);
+
+/* Returns the number of edges in the DCEL. */
+int getDCELEdgeCount(DCEL_t *dcel);
+
+/* Get start vertex of given edge in DCEL. */
+int getDCELEdgeVertexStart(DCEL_t *dcel, int edge);
+
+/* Get end vertex of given edge in DCEL. */
+int getDCELEdgeVertexEnd(DCEL_t *dcel, int edge);
+
+/* Get start vertex of paired given edge in DCEL. */
+int getDCELEdgeVertexPairStart(DCEL_t *dcel, int edge);
+
+/* Get end vertex of paired given edge in DCEL. */
+int getDCELEdgeVertexPairEnd(DCEL_t *dcel, int edge);
+
+/* Check if the DCEL has the given edge. */
+int DCELhasEdge(DCEL_t *dcel, int edge);
+
+/* Check if the DCEL has a pair for the given edge. */
+int DCELhasEdgePair(DCEL_t *dcel, int edge);
+
+/* My own functions */
+
+/* Reads the polygon file and stores the information in the vertices array */
 vertex_t **readPolygon(vertex_t **vertices, FILE *polygonFile, \
     int *numVertices);
+
+/* Frees an array of vertices */
 void freeVertices(vertex_t **vertices, int numVertices);
 
+/* Calculates and returns the equation of a bisector of two points */
+bisector_t *getBisector(vertex_t *pointA, vertex_t *pointB);
+
+/* Euclidian distance between two 2D points */
+double dist_2D(double x1, double y1, double x2, double y2);
+
+/* Returns a value for x that ensures the point (x, y) is at least d (distance) 
+ * away from (x1, y1) whilst ensuring (x, y) is on the line y = mx + c1.
+ */
+double inv_dist_2D(double x1, double y1, double m, double c, double d);
+
+/* Returns a point at least distance away from the midpoint of the bisector given */
+vertex_t *getBisectorPoint(double distance, bisector_t *b);
+
+/* O'Rourke's functions */
+/* Returns -1, 0 or 1, based on the area enclosed by the three points. 0 corresponds
+    to no area enclosed.
+*/
+int areaSign(double sx, double sy, double ex, double ey, double x, double y);
+
+/* Returns 1 if point (x, y) is between (sx, sy) and (ex, ey) */
+int between(double sx, double sy, double ex, double ey, double x, double y);
+
+/* Returns 1 if the point (x, y) is in the line from s(x, y) to e(x, y),
+ * 0 otherwise.
+ */
+int collinear(double sx, double sy, double ex, double ey, double x, double y);
+
+/* Tests if the half edge and bisector are parallel and overlapping, or not
+ * intersecting.
+ */
+enum intersectType parallelIntersects(double heSx, double heSy, \
+    double heEx, double heEy, double bSx, double bSy, \
+    double bEx, double bEy, double *x, double *y);
+
+/* Tests if  */
+enum intersectType intersects(halfEdge_t *he, bisector_t *b, \
+    DCEL_t *dcel, double minLength, double *x, double *y);
+
 #endif
diff --git a/dcel.o b/dcel.o
index 88feeaf..aca2613 100644
Binary files a/dcel.o and b/dcel.o differ
diff --git a/geometry.c b/geometry.c
deleted file mode 100644
index ea75e12..0000000
--- a/geometry.c
+++ /dev/null
@@ -1,56 +0,0 @@
-/* geometry.c
- *
- * Created by Rory Healy (healyr@student.unimelb.edu.au)
- * Created on 9th September 2021
- * Last modified 9th September 2021
- *
- * A small library of functions used in the geometric construction of a 
- * Voronoi diagram.
- *
- */
-
-#ifndef GEOMETRY_HEADER
-#include "geometry.h"
-#endif
-
-#ifndef COMMON_HEADER
-#include "common.h"
-#endif
-
-bisector_t getBisector(vertex_t *pointA, vertex_t *pointB) {
-    bisector_t newBisector;
-    double midpointX = (pointA->x + pointB->x) / 2;
-    double midpointY = (pointA->y + pointB->y) / 2;
-
-    newBisector.x = midpointX;
-    newBisector.y = midpointY;
-
-    /* Calculating bisector slope according to slope of AB */
-    if (pointA->x == pointB->x) {
-        /* The line segment AB has an infinite gradient, 
-         * so the orthogonal line will have zero slope.
-         */
-        newBisector.slope = 0;
-        newBisector.isSlopeInfinite = 0;
-    } else if (pointA->y == pointB->y) {
-        /* The line segment AB has gradient of zero, so 
-         * the orthogonal line will have an infinite slope. 
-         */
-        newBisector.isSlopeInfinite = 1;
-
-        /* Not actually zero, just a placeholder to prevent 
-         * accidental errors 
-         */
-        newBisector.slope = 0;
-    } else {
-        /* Slope of the line segment AB */
-        double segmentSlope = \
-            (pointB->y - pointA->y) / (pointB->x - pointA->x);
-
-        /* Calculate orthogonal slope */
-        newBisector.isSlopeInfinite = 0;
-        newBisector.slope = -1 / segmentSlope;
-    }
-
-    return newBisector;
-}
diff --git a/geometry.h b/geometry.h
deleted file mode 100644
index f199c17..0000000
--- a/geometry.h
+++ /dev/null
@@ -1,19 +0,0 @@
-#ifndef GEOMETRY_HEADER
-#define GEOMETRY_HEADER
-
-typedef struct vertex {
-    double x;
-    double y;
-} vertex_t;
-
-typedef struct bisector {
-    int isSlopeInfinite;    
-    double slope;
-    double x;
-    double y;
-} bisector_t;
-
-/* Calculates and returns the equation of a bisector of two points */
-bisector_t getBisector(vertex_t *vertexA, vertex_t *vertexB); 
-
-#endif
\ No newline at end of file
diff --git a/geometry.o b/geometry.o
deleted file mode 100644
index 5e75cfe..0000000
Binary files a/geometry.o and /dev/null differ
diff --git a/input.o b/input.o
index 325dbb4..a46796f 100644
Binary files a/input.o and b/input.o differ
diff --git a/main.o b/main.o
index af22448..d5ba5ac 100644
Binary files a/main.o and b/main.o differ
diff --git a/output.txt b/output.txt
index 7c243ae..e69de29 100644
--- a/output.txt
+++ b/output.txt
@@ -1,6 +0,0 @@
-y = 0.000000 * (x - 145.600000) + -34.700000
-y = 0.000000 * (x - 145.600000) + -35.200000
-x = 147.100000
-x = 147.100000
-x = 147.600000
-y = 1.000000 * (x - 147.600000) + -33.200000
diff --git a/towers.c b/towers.c
new file mode 100644
index 0000000..0d66bcb
--- /dev/null
+++ b/towers.c
@@ -0,0 +1,97 @@
+
+
+#ifndef TOWERS_HEADER
+#include "towers.h"
+#endif
+
+#define BASE_10 10
+#define MAX_CSV_ENTRY_LEN 512
+#define MAX_FIELD_LEN 128
+#define NUM_CSV_FIELDS 6
+
+tower_t **readTowers(tower_t **towers, FILE *datasetFile, int *numTowers) {
+    /* Maximum length of a single CSV line */
+    size_t lineBufferSize = MAX_CSV_ENTRY_LEN + 1;
+
+    /* Stores the current line from the CSV */
+    char *lineBuffer = malloc(lineBufferSize * sizeof(*lineBuffer));
+    checkNullPointer(lineBuffer);
+
+    int maxSizetowers = 1;
+
+    /* Discard the header line, then read the rest of the CSV */
+    getline(&lineBuffer, &lineBufferSize, datasetFile);
+    while (getline(&lineBuffer, &lineBufferSize, datasetFile) > 0) {
+        /* Check if there enough space in the towers array */
+        if (*numTowers == maxSizetowers) {
+            maxSizetowers *= 2;
+            tower_t **temp = realloc(towers, maxSizetowers * sizeof(*towers));
+            checkNullPointer(temp);
+            towers = temp;
+        }
+
+        /* The current tower being filled in with information */
+        towers[*numTowers] = malloc(sizeof(*towers[*numTowers]));
+
+        readCurrentTower(lineBuffer, towers[*numTowers]);
+        *numTowers += 1;    
+    }
+    free(lineBuffer);
+    return towers;
+}
+
+void readCurrentTower(char *lineBuffer, tower_t *tower) {
+    /* Stores the current CSV field for the current line */
+    char *token = strtok(lineBuffer, ",");
+    size_t tokenLength;
+
+    for (int i = 0; i < NUM_CSV_FIELDS; i++) {
+        tokenLength = strlen(token);
+        switch(i) {
+            /* Case 0, 1, and 3 deal with strings
+             * Case 2 deals with an integer
+             * Case 4 and 5 deal with doubles
+             * Each case is handled seperately to fill in the
+             * tower with no space wasted.
+             */
+            case 0:
+                tower->id = malloc(sizeof(char) * tokenLength + 1);
+                checkNullPointer(tower->id);
+                strcpy(tower->id, token);
+                tower->id[tokenLength] = '\0';
+                break;
+            case 1:
+                tower->postcode = malloc(sizeof(char) * tokenLength + 1);
+                checkNullPointer(tower->postcode);
+                strcpy(tower->postcode, token);
+                tower->postcode[tokenLength] = '\0';
+                break;
+            case 2:
+                tower->population = strtol(token, NULL, BASE_10);
+                break;
+            case 3:
+                tower->manager = malloc(sizeof(char) * tokenLength + 1);
+                checkNullPointer(tower->manager);
+                strcpy(tower->manager, token);
+                tower->manager[tokenLength] = '\0';
+                break;
+            case 4:
+                tower->x = strtod(token, NULL);
+                break;
+            case 5:
+                tower->y = strtod(token, NULL);
+                break;
+        }
+        token = strtok(NULL, ",");
+    }
+}
+
+void freeTowers(tower_t **towers, int numTowers) {
+    for (int i = 0; i < numTowers; i++) {
+        free(towers[i]->id);
+        free(towers[i]->manager);
+        free(towers[i]->postcode);
+        free(towers[i]);
+    }
+    free(towers);
+}
diff --git a/towers.h b/towers.h
new file mode 100644
index 0000000..2e9a696
--- /dev/null
+++ b/towers.h
@@ -0,0 +1,26 @@
+#ifndef COMMON_HEADER
+#include "common.h"
+#endif
+
+#ifndef TOWERS_HEADER
+#define TOWERS_HEADER
+
+typedef struct tower {
+    char *id;
+    char *postcode;
+    char *manager;
+    int population;
+    double x;
+    double y;
+} tower_t;
+
+/* Reads the CSV file and stores the information in the towers array */
+tower_t **readTowers(tower_t **towers, FILE *datasetFile, int *numTowers);
+
+/* Reads the current row from the CSV and converts it into a new tower */
+void readCurrentTower(char *lineBuffer, tower_t *tower);
+
+/* Frees all towers in a towers array */
+void freeTowers(tower_t **towers, int numTowers);
+
+#endif
diff --git a/towers.o b/towers.o
new file mode 100644
index 0000000..04f55ab
Binary files /dev/null and b/towers.o differ
diff --git a/voronoi.c b/voronoi.c
index c012bb0..78ef41f 100644
--- a/voronoi.c
+++ b/voronoi.c
@@ -2,7 +2,7 @@
  *
  * Created by Rory Healy (healyr@student.unimelb.edu.au)
  * Created on 12th August 2021
- * Last modified 9th September 2021
+ * Last modified 11th September 2021
  *
  * Contains functions involving the generation of the Voronoi diagram, and 
  * running each stage of the assignment.
@@ -17,101 +17,9 @@
 #include "voronoi.h"
 #endif
 
-#define BASE_10 10
-#define MAX_CSV_ENTRY_LEN 512
-#define MAX_FIELD_LEN 128
-#define NUM_CSV_FIELDS 6
-
-tower_t **readTowers(tower_t **towers, FILE *datasetFile, int *numTowers) {
-    /* Maximum length of a single CSV line */
-    size_t lineBufferSize = MAX_CSV_ENTRY_LEN + 1;
-
-    /* Stores the current line from the CSV */
-    char *lineBuffer = malloc(lineBufferSize * sizeof(*lineBuffer));
-    checkNullPointer(lineBuffer);
-
-    int maxSizetowers = 1;
-
-    /* Discard the header line, then read the rest of the CSV */
-    getline(&lineBuffer, &lineBufferSize, datasetFile);
-    while (getline(&lineBuffer, &lineBufferSize, datasetFile) > 0) {
-        /* Check if there enough space in the towers array */
-        if (*numTowers == maxSizetowers) {
-            maxSizetowers *= 2;
-            tower_t **temp = realloc(towers, maxSizetowers * sizeof(*towers));
-            checkNullPointer(temp);
-            towers = temp;
-        }
-
-        /* The current tower being filled in with information */
-        towers[*numTowers] = malloc(sizeof(*towers[*numTowers]));
-
-        readCurrentTower(lineBuffer, towers[*numTowers]);
-        *numTowers += 1;    
-    }
-    free(lineBuffer);
-    return towers;
-}
-
-void readCurrentTower(char *lineBuffer, tower_t *tower) {
-    /* Stores the current CSV field for the current line */
-    char *token = strtok(lineBuffer, ",");
-    size_t tokenLength;
-
-    for (int i = 0; i < NUM_CSV_FIELDS; i++) {
-        tokenLength = strlen(token);
-        switch(i) {
-            /* Case 0, 1, and 3 deal with strings
-             * Case 2 deals with an integer
-             * Case 4 and 5 deal with doubles
-             * Each case is handled seperately to fill in the
-             * tower with no space wasted.
-             */
-            case 0:
-                tower->id = malloc(sizeof(char) * tokenLength + 1);
-                checkNullPointer(tower->id);
-                strcpy(tower->id, token);
-                tower->id[tokenLength] = '\0';
-                break;
-            case 1:
-                tower->postcode = malloc(sizeof(char) * tokenLength + 1);
-                checkNullPointer(tower->postcode);
-                strcpy(tower->postcode, token);
-                tower->postcode[tokenLength] = '\0';
-                break;
-            case 2:
-                tower->population = strtol(token, NULL, BASE_10);
-                break;
-            case 3:
-                tower->manager = malloc(sizeof(char) * tokenLength + 1);
-                checkNullPointer(tower->manager);
-                strcpy(tower->manager, token);
-                tower->manager[tokenLength] = '\0';
-                break;
-            case 4:
-                tower->x = strtod(token, NULL);
-                break;
-            case 5:
-                tower->y = strtod(token, NULL);
-                break;
-        }
-        token = strtok(NULL, ",");
-    }
-}
-
-void freeTowers(tower_t **towers, int numTowers) {
-    for (int i = 0; i < numTowers; i++) {
-        free(towers[i]->id);
-        free(towers[i]->manager);
-        free(towers[i]->postcode);
-        free(towers[i]);
-    }
-    free(towers);
-}
-
 vertex_t **readPoints(vertex_t **points, FILE *pointsFile, int *numPoints) {
-    /* Current length of a single points line */
-    size_t lineBufferSize = 50;
+    /* Initial size of buffer is 1 as getline() reallocs as needed */
+    size_t lineBufferSize = 1;
 
     /* Stores the current line from the points file */
     char *lineBuffer = malloc(lineBufferSize * sizeof(*lineBuffer));
@@ -148,6 +56,16 @@ vertex_t **readPoints(vertex_t **points, FILE *pointsFile, int *numPoints) {
     return points;
 }
 
+void stage1PrintBisector(bisector_t *bisector, FILE *outputFile) {
+    if (bisector->isSlopeInfinite) {
+        /* Cannot print infinite slope, so print only x-intercept */
+        fprintf(outputFile, "x = %lf\n", bisector->midX);
+    } else {
+        fprintf(outputFile, "y = %lf * (x - %lf) + %lf\n", \
+            bisector->slope, bisector->midX,  bisector->midY);
+    }
+}
+
 void freePoints(vertex_t **points, int numPoints) {
     for (int i = 0; i < numPoints; i++) {
         free(points[i]);
@@ -168,17 +86,12 @@ void stage1(char *pointsFileName, char *outputFileName) {
 
     /* For each pair, calculate and print the bisector to outputFile */
     for (int i = 0; i < numPoints; i += 2) {
-        bisector_t currBisector = getBisector(points[i], points[i + 1]);
-
-        /* Cannot print infinite slope, so print only x-intercept */
-        if (currBisector.isSlopeInfinite) {
-            fprintf(outputFile, "x = %lf\n", currBisector.x);
-        } else {
-            fprintf(outputFile, "y = %lf * (x - %lf) + %lf\n", \
-                currBisector.slope, currBisector.x,  currBisector.y);
-        }
+        bisector_t *currBisector = getBisector(points[i], points[i + 1]);
+        stage1PrintBisector(currBisector, outputFile);
+        free(currBisector);
     }
 
+    /* Clean up */
     freePoints(points, numPoints);
     fclose(pointsFile);
     fclose(outputFile);
@@ -190,8 +103,21 @@ void stage2(char *pointsFileName, char *polygonFileName, char *outputFileName) {
     polygonFile = safeFileOpen(&polygonFile, polygonFileName, "r");
     outputFile = safeFileOpen(&outputFile, outputFileName, "w");
 
+    /* Points are given as pairs, so initialise 2 points */
+    vertex_t **points = malloc(sizeof(*points) * 2);
+    checkNullPointer(points);
+    int numPoints = 0;
+    points = readPoints(points, pointsFile, &numPoints);
+
+    /* Construct the DCEL from the polygon file */
+    DCEL_t *dcel = readPolygonFile(polygonFileName);
+
+    /* Calculate and print intersections to the output file */
     
 
+    /* Clean up */
+    freePoints(points, numPoints);
+    freeDCEL(dcel);
     fclose(pointsFile);
     fclose(polygonFile);
     fclose(outputFile);
@@ -205,6 +131,7 @@ void stage3(char *dataFileName, char *polygonFileName, char *outputFileName) {
 
 
 
+    /* Clean up */
     fclose(dataFile);
     fclose(polygonFile);
     fclose(outputFile);
@@ -217,7 +144,8 @@ void stage4(char *dataFileName, char *polygonFileName, char *outputFileName) {
     outputFile = safeFileOpen(&outputFile, outputFileName, "w");
 
 
-    
+
+    /* Clean up */
     fclose(dataFile);
     fclose(polygonFile);
     fclose(outputFile);
diff --git a/voronoi.h b/voronoi.h
index 905c1b8..8f059ac 100644
--- a/voronoi.h
+++ b/voronoi.h
@@ -5,27 +5,12 @@
 #ifndef VORONOI_HEADER
 #define VORONOI_HEADER
 
-typedef struct tower {
-    char *id;
-    char *postcode;
-    char *manager;
-    int population;
-    double x;
-    double y;
-} tower_t;
-
-/* Reads the CSV file and stores the information in the towers array */
-tower_t **readTowers(tower_t **towers, FILE *datasetFile, int *numTowers);
-
-/* Reads the current row from the CSV and converts it into a new tower */
-void readCurrentTower(char *lineBuffer, tower_t *tower);
-
-/* Frees all towers in a towers array */
-void freeTowers(tower_t **towers, int numTowers);
-
 /* Reads a points file and stores the information in the points array */
 vertex_t **readPoints(vertex_t **points, FILE *pointsFile, int *numPoints);
 
+/* Prints bisectors to an output file */
+void stage1PrintBisector(bisector_t *bisector, FILE *outputFile);
+
 /* Frees all points from a points array */
 void freePoints(vertex_t **points, int numPoints);
 
diff --git a/voronoi.o b/voronoi.o
index 707d2cc..4b32af2 100644
Binary files a/voronoi.o and b/voronoi.o differ
diff --git a/voronoi2 b/voronoi2
index 9096ef9..a1adcbe 100755
Binary files a/voronoi2 and b/voronoi2 differ