/* dcel.c
 *
 * Created by Rory Healy (healyr@student.unimelb.edu.au)
 * Created on 25th August 2021
 * Last modified 13th September 2021
 *
 * 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

void freePoints(vertex_t **points, int numPoints) {
    for (int i = 0; i < numPoints; i++) {
        free(points[i]);
    }
    free(points);
}

void freeBisectors(bisector_t **bisectors, int numBisectors) {
    for (int i = 0; i < numBisectors; i++) {
        free(bisectors[i]->mid);
        free(bisectors[i]);
    }
    free(bisectors);
}

void freeIntersections(intersection_t **intersections, int numIntersections) {
    for (int i = 0; i < numIntersections; i++) {
        free(intersections[i]->fromPoint);
        free(intersections[i]->toPoint);
        free(intersections[i]);
    }
    free(intersections);
}

bisector_t **getBisectors(bisector_t **bisectors, vertex_t **points, \
    int numBisectors) {
    for (int i = 0; i < numBisectors; i++) {
        bisectors[i] = malloc(sizeof(*bisectors[i]));
        bisectors[i]->mid = malloc(sizeof(*bisectors[i]->mid));

        /* Calculate midpoint of the two points */
        bisectors[i]->mid->x = (points[2 * i]->x + points[2 * i + 1]->x) / 2;
        bisectors[i]->mid->y = (points[2 * i]->y + points[2 * i + 1]->y) / 2;

        /* Calculating bisector slope according to slope of AB */
        if (points[2 * i]->x == points[2 * i + 1]->x) {
            /* The line segment AB has an infinite gradient, 
             * so the orthogonal line will have zero slope.
             */
            bisectors[i]->slope = 0;
            bisectors[i]->isSlopeInfinite = 0;
        } else if (points[2 * i]->y == points[2 * i + 1]->y) {
            /* The line segment AB has gradient of zero, so 
             * the orthogonal line will have an infinite slope. 
             */
            bisectors[i]->isSlopeInfinite = 1;

            /* Not actually zero, just a placeholder to prevent 
             * accidental errors 
             */
            bisectors[i]->slope = 0;
        } else {
            /* The line segment AB has a non-zero and finite gradient, 
             * so the gradient of the bisector can be calculated.
             */
            bisectors[i]->isSlopeInfinite = 0;
            bisectors[i]->slope = -1 / \
                ((points[2 * i + 1]->y - points[2 * i]->y) / \
                (points[2 * i + 1]->x - points[2 * i]->x));
        }
    }
    
    return bisectors;
}

vertex_t *getBisectorPoint(double distance, bisector_t *b, int direction) {
    vertex_t *returnPoint = malloc(sizeof(*returnPoint));

    if (b->isSlopeInfinite) {
        /* Vertical line - just add vertical distance */
        returnPoint->x = b->mid->x;
        returnPoint->y = b->mid->y + (distance * direction);
    } else if (b->slope == 0) {
        /* Horizontal line - just add horizontal distance */
        returnPoint->x = b->mid->x + (distance * direction);
        returnPoint->y = b->mid->y;
    } else {
        /* Not horizontal or vertical - add distance to x, then find 
         * y-intercept of the bisector, and plug it in to y = mx + c to find
         * the corresponding y-value.
         */
        double c = b->mid->y - b->slope * b->mid->x;
        returnPoint->x = b->mid->x + (distance * direction);
        returnPoint->y = b->slope * returnPoint->x + c;
    }

    return returnPoint;
}

intersection_t *newIntersection() {
    intersection_t *intersection = malloc(sizeof(*intersection));
    checkNullPointer(intersection);

    intersection->fromPoint = malloc(sizeof(*intersection->fromPoint));
    checkNullPointer(intersection->fromPoint);

    intersection->toPoint = malloc(sizeof(*intersection->toPoint));
    checkNullPointer(intersection->toPoint);

    return intersection;
}

vertex_t **readPoints(vertex_t **points, FILE *pointsFile, int *numPoints) {
    /* 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));
    checkNullPointer(lineBuffer);

    /* Assuming that pointsFile is valid, there must be at least two points */
    int maxSizePoints = 2;

    while (getline(&lineBuffer, &lineBufferSize, pointsFile) > 0) {
        /* Ensure there is enough space in the points array */
        if (*numPoints == maxSizePoints) {
            maxSizePoints *= 2;
            vertex_t **temp = realloc(points, maxSizePoints * sizeof(*points));
            checkNullPointer(temp);
            points = temp;
        }

        double xCoordinateA, yCoordinateA, xCoordinateB, yCoordinateB;
        sscanf(lineBuffer, "%lf %lf %lf %lf", &xCoordinateA, &yCoordinateA, \
            &xCoordinateB, &yCoordinateB);
        
        points[*numPoints] = malloc(sizeof(*points[*numPoints]));
        points[*numPoints]->x = xCoordinateA;
        points[*numPoints]->y = yCoordinateA;
        *numPoints += 1;

        points[*numPoints] = malloc(sizeof(*points[*numPoints]));
        points[*numPoints]->x = xCoordinateB;
        points[*numPoints]->y = yCoordinateB;
        *numPoints += 1;
    }

    free(lineBuffer);
    return points;
}

/* 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)

/* 
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 */
    vertex_t *startPoint = getBisectorPoint(minLength, b, -1);
    vertex_t *endPoint = getBisectorPoint(minLength, b, 1);
    double bSx = startPoint->x;
    double bSy = startPoint->y;
    double bEx = endPoint->x;
    double bEy = endPoint->y;

    free(startPoint);
    free(endPoint);
    
    /* Fill in segment. */

    /* 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;
    }
}

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;
}

intersection_t **getIntersections(intersection_t **intersections, \
    int *numIntersections, bisector_t **bisectors, int numBisectors, \
    DCEL_t *dcel, int face, double minLength) {
    int maxSizeIntersections = 1;

    for (int i = 0; i < numBisectors; i++) {
        intersection_t *currentIntersection = NULL;

        /* Intersection coordinates */
        double x, y;

        /* Flag is raised when first intersection is found */
        int isIntersectionFound = 0;
        
        halfEdge_t *startHE = (dcel->faces)[face].start;   
        halfEdge_t *currentHE = startHE;
        int startVertex = startHE->startVertex;
        int first = 1;

        /* Loop through the face until the starting vertex is reached */
        while (first || currentHE->startVertex != startVertex) {
            enum intersectType typeOfIntersection = \
                intersects(currentHE, bisectors[i], dcel, minLength, &x, &y);
            
            switch (typeOfIntersection) {
                case INTERSECT:
                case SAME_LINE_OVERLAP:
                case ENDS_OVERLAP:
                    if (!isIntersectionFound) {
                        /* First point of intersection */
                        isIntersectionFound = 1;
                        currentIntersection = newIntersection();
                        currentIntersection->fromPoint->x = x;
                        currentIntersection->fromPoint->y = y;
                        currentIntersection->fromEdge = currentHE->edge;
                    } else {
                        /* Second point of intersection */
                        currentIntersection->toPoint->x = x;
                        currentIntersection->toPoint->y = y;
                        currentIntersection->toEdge = currentHE->edge;
                    }
                    break;
                case DOESNT_INTERSECT:
                default:
                    break;
            }
            currentHE = currentHE->next;
            first = 0;
        }

        /* If there is an intersection, add it to the array */
        if (currentIntersection != NULL) {
            /* Ensure there is enough space in the array */
            if (*numIntersections == maxSizeIntersections) {
                maxSizeIntersections *= 2;
                intersection_t **temp = realloc(intersections, maxSizeIntersections * sizeof(*intersections));
                checkNullPointer(temp);
                intersections = temp;
            }

            intersections[*numIntersections] = currentIntersection;
            *numIntersections += 1;
        }
    }

    return intersections;
}

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);
}