Initial import

.gitignore

@@ -0,0 +1,3 @@
+.idea*
+lib*
+node_modules*

Makefile

@@ -0,0 +1,9 @@
+
+.PHONY: all
+all:
+	rm -rf lib
+	tsc
+
+.PHONY: run
+run: all
+	node lib/test.js

package.json

@@ -0,0 +1,5 @@
+{
+	"dependencies": {
+		"nodeplotlib": "^0.7.5"
+	}
+}

pnpm-lock.yaml

@@ -0,0 +1,25 @@
+lockfileVersion: 5.3
+
+specifiers:
+  nodeplotlib: ^0.7.5
+
+dependencies:
+  nodeplotlib: 0.7.5
+
+packages:
+
+  /@types/d3/3.5.45:
+    resolution: {integrity: sha512-wLICfMtjDEoAJie1MF6OuksAzOapRXgJy+l5HQVpyC1yMAlvHz2QKrrasUHru8xD6cbgQNGeO+CeyjOlKtly2A==}
+    dev: false
+
+  /@types/plotly.js/1.54.16:
+    resolution: {integrity: sha512-Kli9UqZwO12cvJ4dT2kGaBtpgq96y8TPsvsB87XTCO+jEMMVy0OZ4WMW+8plS/lBKCtMNtwRNYNlD3+8Stgn2Q==}
+    dependencies:
+      '@types/d3': 3.5.45
+    dev: false
+
+  /nodeplotlib/0.7.5:
+    resolution: {integrity: sha512-vCebCtxvzOkN8VnFFyeMShOZmVW/pNPHtl3WuqvvcXZylBLKoZXpmgDuL/M+OcAsBCgMJMA5pXfBHydnjwRrJw==}
+    dependencies:
+      '@types/plotly.js': 1.54.16
+    dev: false

src/delaunayRefine.ts

@@ -0,0 +1,92 @@
+import {addBoundingSquareTo, Angle, Graph, GraphBoundary, Point, rad2deg, rad2degp, Segment, Triangle} from "./pslg";
+
+export function addTriangle(graph: Graph, boundary: GraphBoundary, segment: Segment) {
+    let p1 = segment.from
+    let p2 = segment.to
+
+    // Want to order p1 and p2 such that p2 is the right neighbor of p1
+    if ( !boundary.rightNeighbor(p1).is(p2) ) {
+        let temp = p1
+        p1 = p2
+        p2 = temp
+    }
+
+    // Get a segment with the ordered from/to
+    // So, looking from p1 to p2, points on the "left side" of the line
+    // are within the boundary
+    const orderedSegment = new Segment(p1, p2)
+
+    // Collect all points that could form valid triangles
+    const possiblePoints: {angle: number, point: Point}[] = []
+
+    // Find some p3 such that p1p2p3 is Delaunay (<p1p3p2 is the largest angle)
+    for ( const point of graph.points ) {
+        if ( point.is(p1) || point.is(p2) ) continue
+
+        // Check if the point is within boundary -- excluding boundary points themselves
+        if ( !boundary.containsPointWithin(point) ) continue
+
+        // TODO Check if the point is within the search region
+
+        // Check if the point is on the left-side of segment
+        if ( !orderedSegment.pointIsLeftOf(point) ) continue
+
+        possiblePoints.push({
+            angle: Triangle.getAngleFrom(p1, point, p2),
+            point,
+        })
+    }
+
+    // Sort the possible points to find the one with the largest angle
+    const sortedPossiblePoints = possiblePoints.sort((a, b) => {
+        const aDeg = rad2degp(a.angle)
+        const bDeg = rad2degp(b.angle)
+        return bDeg - aDeg
+    })
+
+    // Go through the sorted points and find one that forms a valid triangle
+    for ( const {point} of sortedPossiblePoints ) {
+        const seg13 = new Segment(p1, point)
+        const seg23 = new Segment(p2, point)
+
+        // Make sure the legs that would form from this triangle wouldn't intersect
+        // with the boundary segment except at p1 and p2
+        if ( segment.getIntersectionWithin(seg13) || segment.getIntersectionWithin(seg23) ) {
+            continue
+        }
+
+        // We found a valid 3rd point for this triangle!
+        const graphSeg13 = graph.findExistingSegmentOrAdd(seg13)
+        const graphSeg23 = graph.findExistingSegmentOrAdd(seg23)
+        const graphSeg12 = graph.findExistingSegmentOrAdd(segment)
+
+        graph.findExistingTriangleOrAdd(new Triangle([graphSeg13, graphSeg23, graphSeg12]))
+        break
+    }
+}
+
+export function delaunay(originalGraph: Graph): Graph {
+    const graph = originalGraph.clone()
+    const boundary = addBoundingSquareTo(graph)
+
+    // for ( const segment of boundary.segments ) {
+    //     addTriangle(graph, boundary, segment)
+    //     break
+    // }
+    // const segment = boundary.segments[0]
+    // addDelaunayTriangle(graph, boundary, segment)
+
+    for ( const point of graph.points ) {
+
+    }
+
+    return graph
+}
+
+export function delaunayRefine(originalGraph: Graph, minAngle: Angle): Graph {
+    const graph = originalGraph.clone()
+
+    addBoundingSquareTo(graph)
+
+    return graph
+}

src/linear.ts

@@ -0,0 +1,58 @@
+type TupleOf<T, N extends number> = N extends N ? number extends N ? T[] : _TupleOf<T, N, []> : never
+type _TupleOf<T, N extends number, R extends unknown[]> = R['length'] extends N ? R : _TupleOf<T, N, [T, ...R]>
+
+export class Vector<N extends number> {
+    public readonly values: TupleOf<number, N>
+
+    constructor(
+        public readonly length: N
+    ) {
+        this.values = Array(length).fill(0) as TupleOf<number, N>
+    }
+}
+
+export class Matrix<R extends number, C extends number> {
+    public readonly values: TupleOf<TupleOf<number, C>, R>
+
+    constructor(
+        public readonly rows: R,
+        public readonly cols: C,
+    ) {
+        this.values = Array(rows).fill(0).map(() => Array(cols).fill(0) as TupleOf<number, C>) as TupleOf<TupleOf<number, C>, R>
+    }
+
+    isSquare() {
+        return Number(this.rows) === Number(this.cols)
+    }
+
+    set(row: number, col: number, value: number) {
+        this.validateCoordinates(row, col)
+        this.values[row-1 as R][col-1 as C] = value
+    }
+
+    get(row: number, col: number): number {
+        this.validateCoordinates(row, col)
+        return this.values[row-1 as R][col-1 as C] as number
+    }
+
+    getDeterminant() {
+        if ( !this.isSquare() ) {
+            throw new RangeError('Cannot find determinant of non-square matrix.')
+        }
+
+        const determinant = (m: number[][]): number =>
+            m.length == 1 ? m[0][0] :
+                m.length == 2 ? m[0][0]*m[1][1]-m[0][1]*m[1][0] :
+                    m[0].reduce((r,e,i) =>
+                        r+(-1)**(i+2)*e*determinant(m.slice(1).map(c =>
+                            c.filter((_,j) => i != j))),0)
+
+        return determinant(this.values)
+    }
+
+    protected validateCoordinates(row: number, col: number) {
+        if ( row < 1 || row > Number(this.rows) || col < 1 || col > Number(this.cols) ) {
+            throw new RangeError('Invalid matrix coordinates')
+        }
+    }
+}

src/pslg.ts

@@ -0,0 +1,766 @@
+import {Matrix} from "./linear";
+
+export type Angle = number
+
+export function deg2rad(degrees: number): number {
+    return degrees * (Math.PI / 180)
+}
+
+export function rad2deg(radian: number): number {
+    return radian * (180 / Math.PI)
+}
+
+export function rad2degp(radian: number): number {
+    const deg = rad2deg(radian)
+    return deg < 0 ? deg + 180 : deg
+}
+
+export interface Coordinate {
+    x: number,
+    y: number,
+}
+
+export class Point {
+    public static from(x: number, y: number): Point {
+        return new Point({x, y})
+    }
+
+    constructor(
+        public readonly coordinate: Coordinate,
+        public readonly name?: string,
+    ) {}
+
+    get x() {
+        return this.coordinate.x
+    }
+
+    get y() {
+        return this.coordinate.y
+    }
+
+    clone(): Point {
+        return new Point({...this.coordinate}, this.name)
+    }
+
+    is(x: Point) {
+        return (
+            x.coordinate.x === this.coordinate.x
+            && x.coordinate.y === this.coordinate.y
+        )
+    }
+
+    public static midpoint(a: Point, b: Point): Point {
+        const x = (a.x + b.x) / 2
+        const y = (a.y + b.y) / 2
+        return Point.from(x, y)
+    }
+
+    public static distinct(points: Point[]): Point[] {
+        const distinct: Point[] = []
+        for ( const point of points ) {
+            if ( !distinct.some(other => other.is(point)) ) {
+                distinct.push(point)
+            }
+        }
+
+        return distinct
+    }
+}
+
+export class Segment {
+    constructor(
+        public readonly from: Point,
+        public readonly to: Point,
+    ) {
+        if ( from.is(to) ) {
+            console.log('new Segment', [from.coordinate, to.coordinate])
+            throw new RangeError('Cannot create segment with same from and to points')
+        }
+    }
+
+    get slope() {
+        return (this.to.y - this.from.y) / (this.to.x - this.from.x)
+    }
+
+    get yIntercept() {
+        return this.getXAtY(0)
+    }
+
+    getXAtY(y: number) {
+        const y1 = this.from.y
+        const x1 = this.from.x
+        const slope = this.slope
+        return ((y - y1) + (slope * x1)) / slope
+    }
+
+    getYAtX(x: number) {
+        const y1 = this.from.y
+        const x1 = this.from.x
+        const slope = this.slope
+        return (slope * (x - x1)) + y1
+    }
+
+    isHorizontal() {
+        return Math.abs(this.slope) === 0
+    }
+
+    isVertical() {
+        return Math.abs(this.slope) === Infinity
+    }
+
+    clone(): Segment {
+        return new Segment(this.from.clone(), this.to.clone())
+    }
+
+    is(x: Segment) {
+        return (
+            (
+                x.from.is(this.from)
+                && x.to.is(this.to)
+            )
+            || (
+                x.from.is(this.to)
+                && x.to.is(this.from)
+            )
+        )
+    }
+
+    containsPoint(x: Point) {
+        if ( this.isHorizontal() ) {
+            return (
+                x.y === this.ymax
+                && x.x >= this.xmin
+                && x.x <= this.xmax
+            )
+        }
+
+        if ( this.isVertical() ) {
+            return (
+                x.x === this.xmax
+                && x.y >= this.ymin
+                && x.y <= this.ymax
+            )
+        }
+
+        return (
+            this.getYAtX(x.x) === x.y
+            && this.getXAtY(x.y) === x.x
+        )
+    }
+
+    hasPoint(x: Point) {
+        return !(
+            this.from.is(x)
+            || this.to.is(x)
+        )
+    }
+
+    getIntersectionWith(x: Segment): Point | undefined {
+        const slope = this.slope
+        const xSlope = x.slope
+
+        if ( slope === xSlope ) return
+
+        const b1 = this.yIntercept
+        const b2 = x.yIntercept
+        const m1 = this.slope
+        const m2 = x.slope
+
+        if ( Math.abs(m1) === 0 && Math.abs(m2) === Infinity ) {
+            if ( x.xmin <= this.xmax && x.xmin >= this.xmin ) {
+                if ( this.ymin <= x.ymax && this.ymin >= x.ymin ) {
+                    return Point.from(x.xmin, this.ymin)
+                }
+            }
+        }
+
+        if ( Math.abs(m1) === Infinity && Math.abs(m2) === 0 ) {
+            if ( this.xmin <= x.xmax && this.xmin >= x.xmin ) {
+                if ( x.ymin <= this.ymax && x.ymin >= this.ymin ) {
+                    return Point.from(this.xmin, x.ymin)
+                }
+            }
+        }
+
+        const intersectX = (b1 - b2) / (m2 - m1)
+        const intersectY = this.getYAtX(intersectX)
+        if ( intersectX <= Math.max(this.from.x, this.to.x) && intersectX >= Math.min(this.from.x, this.to.x) ) {
+            return Point.from(intersectX, intersectY)
+        }
+    }
+
+    getIntersectionWithin(x: Segment): Point | undefined {
+        const slope = this.slope
+        const xSlope = x.slope
+
+        // FIXME account for overlapping parallel lines
+        if ( slope === xSlope ) return
+
+        const b1 = this.yIntercept
+        const b2 = x.yIntercept
+        const m1 = this.slope
+        const m2 = x.slope
+
+        const intersectX = (b1 - b2) / (m2 - m1)
+        const intersectY = this.getYAtX(intersectX)
+        if ( intersectX < Math.max(this.from.x, this.to.x) && intersectX > Math.min(this.from.x, this.to.x) ) {
+            return Point.from(intersectX, intersectY)
+        }
+    }
+
+    getLength(): number {
+        return Math.sqrt(
+            Math.pow(this.to.x - this.from.x, 2)
+            + Math.pow(this.to.y - this.from.y, 2)
+        )
+    }
+
+    /**
+     * Checks if a given point is to the left of this segment, if you were
+     * to look at the segment standing at "from" toward "to".
+     *
+     * This returns true if the point is on the segment itself.
+     *
+     * @param x
+     */
+    pointIsLeftOf(x: Point) {
+        // TODO replace with 3x3 determinant? See C. Du pg 24 (2)
+        const det = (this.to.x - this.from.x) * (x.y - this.from.y) - (this.to.y - this.from.y) * (x.x - this.from.x)
+        return det >= 0
+    }
+
+    split(): [Point, Segment, Segment] {
+        const midpoint = Point.midpoint(this.from, this.to)
+        return [midpoint, ...this.splitAt(midpoint)]
+    }
+
+    splitAt(x: Point): [Segment, Segment] {
+        if ( !this.containsPoint(x) ) {
+            console.log(
+                'splitAt',
+                [this.from.coordinate, this.to.coordinate],
+                x.coordinate,
+            )
+            throw new RangeError('Cannot split segment on point that does not occur on segment')
+        }
+
+        const segment1 = new Segment(this.from, x)
+        const segment2 = new Segment(x, this.to)
+        return [segment1, segment2]
+    }
+
+    getMidpoint(): Point {
+        return Point.midpoint(this.from, this.to)
+    }
+
+    get xmin() {
+        return Math.min(this.from.x, this.to.x)
+    }
+
+    get xmax() {
+        return Math.max(this.from.x, this.to.x)
+    }
+
+    get ymin() {
+        return Math.min(this.from.y, this.to.y)
+    }
+
+    get ymax() {
+        return Math.max(this.from.y, this.to.y)
+    }
+}
+
+export class Circle {
+    constructor(
+        public readonly center: Point,
+        public readonly radius: number,
+    ) { }
+
+    containsPoint(x: Point) {
+        const distance = (
+            Math.pow(x.x - this.center.x, 2)
+            + Math.pow(x.y - this.center.y, 2)
+        )
+
+        return distance <= Math.pow(this.radius, 2)
+    }
+
+    containsPointWithin(x: Point) {
+        const distance = (
+            Math.pow(x.x - this.center.x, 2)
+            + Math.pow(x.y - this.center.y, 2)
+        )
+
+        return distance < Math.pow(this.radius, 2)
+    }
+}
+
+export class Trapezoid {
+    constructor(
+        public readonly points: [Point, Point, Point, Point],
+        public readonly segments: [Segment, Segment, Segment, Segment],
+    ) {
+        if ( !this.validateSegments() ) {
+            throw new Error('Attempted to create Trapezoid with invalid points/segments.')
+        }
+    }
+
+    public validateSegments() {
+        const points: Point[] = []
+        this.segments.some(side => {
+            if ( !points.some(point => point.is(side.from)) ) points.push(side.from)
+            if ( !points.some(point => point.is(side.to)) ) points.push(side.to)
+        })
+
+        if ( !Trapezoid.distinctPoints(...points as [Point, Point, Point, Point]) ) {
+            return false
+        }
+
+        if ( points.length !== 4 ) {
+            return false
+        }
+
+        const sidesByPoint: Segment[][] = [[], [], [], []]
+        for ( const side of this.segments ) {
+            const toPointIndex = points.findIndex(point => point.is(side.to))
+            const fromPointIndex = points.findIndex(point => point.is(side.from))
+
+            if ( !sidesByPoint[toPointIndex].some(seg => seg.is(side)) )
+                sidesByPoint[toPointIndex].push(side)
+
+            if ( !sidesByPoint[fromPointIndex].some(seg => seg.is(side)) )
+                sidesByPoint[fromPointIndex].push(side)
+        }
+
+        return sidesByPoint.every(group => group.length === 2)
+    }
+
+    public static distinctPoints(p1: Point, p2: Point, p3: Point, p4: Point) {
+        return Point.distinct([p1, p2, p3, p4]).length === 4
+    }
+}
+
+/**
+ * Use the convention that:
+ *
+ * Vertex a: 0-from / 2-to
+ * Vertex b: 1-from / 0-to
+ * Vertex c: 2-from / 1-to
+ */
+export class Triangle {
+    get a(): Point {
+        return this.sides[0].from
+    }
+
+    get b(): Point {
+        return this.sides[1].from
+    }
+
+    get c(): Point {
+        return this.sides[2].from
+    }
+
+    get orderedSides(): [Segment, Segment, Segment] {
+        return this.sides.sort((a, b) => {
+            if ( a.from.x !== b.from.x ) {
+                return a.from.x - b.from.x
+            }
+
+            return a.from.y - b.from.y
+        })
+    }
+
+    constructor(
+        public readonly sides: [Segment, Segment, Segment]
+    ) {
+        if ( !this.validateSides() ) {
+            console.log(
+                [this.sides[0].from.coordinate, this.sides[0].to.coordinate],
+                [this.sides[1].from.coordinate, this.sides[1].to.coordinate],
+                [this.sides[2].from.coordinate, this.sides[2].to.coordinate],
+            )
+            throw new Error('Tried to create Triangle with invalid sides.')
+        }
+    }
+
+    public validateSides(): boolean {
+        const points: Point[] = []
+        this.sides.some(side => {
+            if ( !points.some(point => point.is(side.from)) ) points.push(side.from)
+            if ( !points.some(point => point.is(side.to)) ) points.push(side.to)
+        })
+
+        if ( !Triangle.distinctPoints(...points as [Point, Point, Point]) ) {
+            return false
+        }
+
+        if ( points.length !== 3 ) {
+            return false
+        }
+
+        const sidesByPoint: Segment[][] = [[], [], []]
+        for ( const side of this.sides ) {
+            const toPointIndex = points.findIndex(point => point.is(side.to))
+            const fromPointIndex = points.findIndex(point => point.is(side.from))
+
+            if ( !sidesByPoint[toPointIndex].some(seg => seg.is(side)) )
+                sidesByPoint[toPointIndex].push(side)
+
+            if ( !sidesByPoint[fromPointIndex].some(seg => seg.is(side)) )
+                sidesByPoint[fromPointIndex].push(side)
+        }
+
+        return sidesByPoint.every(group => group.length === 2)
+    }
+
+    clone(): Triangle {
+        const [s1, s2, s3] = this.sides.map(x => x.clone())
+        return new Triangle([s1, s2, s3])
+    }
+
+    is(x: Triangle) {
+        return this.orderedSides.every((side, i) => x.orderedSides[i].is(side))
+    }
+
+    hasSegment(x: Segment) {
+        return !this.sides.some(side => side.is(x))
+    }
+
+    hasPoint(x: Point) {
+        return !this.sides.some(side => side.hasPoint(x))
+    }
+
+    /** Get the angles of the vertices a, b, c, respectively. */
+    getAngles(): [Angle, Angle, Angle] {
+        const a = Triangle.getAngleFrom(this.b, this.a, this.c)
+        const b = Triangle.getAngleFrom(this.a, this.b, this.c)
+        const c = Triangle.getAngleFrom(this.b, this.c, this.a)
+        return [a, b, c]
+    }
+
+    getMinimumAngle(): Angle {
+        return deg2rad(
+            Math.min(...this.getAngles().map(x => rad2degp(x)))
+        )
+    }
+
+    /** Get the points of the triangle a, b, c, respectively. */
+    getPoints(): [Point, Point, Point] {
+        return [this.a, this.b, this.c]
+    }
+
+    getCircumcenter(): Point {
+        const [pointA, pointB, pointC] = this.getPoints()
+        const [angleA, angleB, angleC] = this.getAngles()
+
+        const [sin2A, sin2B, sin2C] = this.getAngles().map(x => Math.sin(2 * x))
+
+        const x = ((pointA.x * sin2A) + (pointB.x * sin2B) + (pointC.x * sin2C))
+            / (sin2A + sin2B + sin2C)
+
+        const y = ((pointA.y * sin2A) + (pointB.y * sin2B) + (pointC.y * sin2C))
+            / (sin2A + sin2B + sin2C)
+
+        return Point.from(x, y)
+    }
+
+    getCircumcenterCircle(): Circle {
+        const center = this.getCircumcenter()
+        const radius = (new Segment(center, this.a)).getLength()
+        return new Circle(center, radius)
+    }
+
+    public static distinctPoints(p1: Point, p2: Point, p3: Point) {
+        return Point.distinct([p1, p2, p3]).length === 3
+    }
+
+    public static getAngleFrom(p2: Point, p1: Point, p3: Point): Angle {
+        const numerator = p2.y * (p1.x - p3.x) + p1.y * (p3.x - p2.x) + p3.y * (p2.x - p1.x)
+        const denominator = (p2.x - p1.x) * (p1.x - p3.x) + (p2.y - p1.y) * (p1.y - p3.y)
+        const radio = numerator / denominator
+
+        return Math.atan(radio)
+    }
+}
+
+export enum GraphDirection {
+    UP,
+    DOWN,
+    LEFT,
+    RIGHT,
+}
+
+export interface SegmentWithIntersection {
+    segment: Segment,
+    intersect: Point,
+}
+
+export function getDirectionSorter(direction: GraphDirection) {
+    if ( direction === GraphDirection.UP ) {
+        return (a: SegmentWithIntersection, b: SegmentWithIntersection) => {
+            return a.intersect.y - b.intersect.y  // get the lowest y-value first
+        }
+    } else if ( direction === GraphDirection.DOWN ) {
+        return (a: SegmentWithIntersection, b: SegmentWithIntersection) => {
+            return b.intersect.y - a.intersect.y  // get the highest y-value first
+        }
+    } else if ( direction === GraphDirection.LEFT ) {
+        return (a: SegmentWithIntersection, b: SegmentWithIntersection) => {
+            return b.intersect.x - a.intersect.x  // get the right-most x value first
+        }
+    } else {
+        return (a: SegmentWithIntersection, b: SegmentWithIntersection) => {
+            return a.intersect.x - b.intersect.x  // get the left-most x-value first
+        }
+    }
+}
+
+export class GraphBoundary {
+    constructor(
+        public readonly points: Point[],
+        public readonly segments: Segment[],
+    ) { }
+
+    getLeftBoundary(): Segment {
+        const left = this.segments.find(segment => segment.from.x === this.xmin && segment.to.x === this.xmin)
+        if ( !left ) throw new RangeError('Unable to find left boundary')
+        return left
+    }
+
+    getRightBoundary(): Segment {
+        const right = this.segments.find(segment => segment.from.x === this.xmax && segment.to.x === this.xmax)
+        if ( !right ) throw new RangeError('Unable to find right boundary')
+        return right
+    }
+
+    getUpperBoundary(): Segment {
+        const top = this.segments.find(segment => segment.from.y === this.ymax && segment.to.y === this.ymax)
+        if ( !top ) throw new RangeError('Unable to find upper boundary')
+        return top
+    }
+
+    getLowerBoundary(): Segment {
+        const bottom = this.segments.find(segment => segment.from.y === this.ymin && segment.to.y === this.ymin)
+        if ( !bottom ) throw new RangeError('Unable to find lower boundary')
+        return bottom
+    }
+
+    getBoundary(direction: GraphDirection) {
+        if ( direction === GraphDirection.UP ) return this.getUpperBoundary()
+        else if ( direction === GraphDirection.DOWN ) return this.getLowerBoundary()
+        else if ( direction === GraphDirection.LEFT ) return this.getLeftBoundary()
+        else return this.getRightBoundary()
+    }
+
+    leftNeighbor(to: Point): Point {  // NP1
+        const i = this.points.findIndex(point => point.is(to))
+        if ( i < 0 ) throw new Error('Cannot find neighbor of point not on boundary.')
+
+        const neighbor = (i-1) < 0 ? this.points.length - 1 : i - 1
+        return this.points[neighbor]
+    }
+
+    rightNeighbor(to: Point) {  // NP2
+        const i = this.points.findIndex(point => point.is(to))
+        if ( i < 0 ) throw new Error('Cannot find neighbor of point not on boundary.')
+
+        const neighbor = (i+1) >= this.points.length ? 0 : i + 1
+        return this.points[neighbor]
+    }
+
+    isBoundaryPoint(x: Point) {
+        return this.points.some(point => point.is(x))
+    }
+
+    isBoundarySegment(x: Segment) {
+        return this.segments.some(segment => segment.is(x))
+    }
+
+    containsPoint(x: Point) {
+        return (
+            x.x >= this.xmin
+            && x.x <= this.xmax
+            && x.y >= this.ymin
+            && x.y <= this.ymax
+        )
+    }
+
+    containsPointWithin(x: Point) {
+        return (
+            x.x > this.xmin
+            && x.x < this.xmax
+            && x.y > this.ymin
+            && x.y < this.ymax
+        )
+    }
+
+    // TODO centralize into a ContainsPointsAndSegments base class
+    get xmin() {
+        return Math.min(
+            ...this.points.map(point => point.x)
+        )
+    }
+
+    get xmax() {
+        return Math.max(
+            ...this.points.map(point => point.x)
+        )
+    }
+
+    get ymin() {
+        return Math.min(
+            ...this.points.map(point => point.y)
+        )
+    }
+
+    get ymax() {
+        return Math.max(
+            ...this.points.map(point => point.y)
+        )
+    }
+}
+
+export class Graph {
+    public readonly vertices: Point[] = []
+    public readonly edges: Segment[] = []
+
+    constructor(
+        public points: Point[] = [],
+        public segments: Segment[] = [],
+        public triangles: Triangle[] = [],
+    ) {
+        this.vertices = [...points.map(x => x.clone())]
+        this.edges = [...segments.map(x => x.clone())]
+    }
+
+    get xmin() {
+        return Math.min(
+            ...this.points.map(point => point.x)
+        )
+    }
+
+    get xmax() {
+        return Math.max(
+            ...this.points.map(point => point.x)
+        )
+    }
+
+    get ymin() {
+        return Math.min(
+            ...this.points.map(point => point.y)
+        )
+    }
+
+    get ymax() {
+        return Math.max(
+            ...this.points.map(point => point.y)
+        )
+    }
+
+    get span() {
+        return Math.max(
+            this.xmax - this.xmin,
+            this.ymax - this.ymin
+        )
+    }
+
+    get center(): Point {
+        const xmax = Point.from(this.xmax, 0)
+        const xmin = Point.from(this.xmin, 0)
+        const xmid = Point.midpoint(xmin, xmax)
+
+        const ymax = Point.from(0, this.ymax)
+        const ymin = Point.from(0, this.ymin)
+        const ymid = Point.midpoint(ymin, ymax)
+
+        return Point.from(xmid.x, ymid.y)
+    }
+
+    findExistingPointOrAdd(x: Point): Point {
+        const existing = this.points.find(point => point.is(x))
+        if ( existing ) return existing
+
+        this.points.push(x)
+        return x
+    }
+
+    findExistingSegmentOrAdd(x: Segment): Segment {
+        const existing = this.segments.find(segment => segment.is(x))
+        if ( existing ) return existing
+
+        this.segments.push(x)
+        return x
+    }
+
+    findExistingTriangleOrAdd(x: Triangle): Triangle {
+        const existing = this.triangles.find(triangle => triangle.is(x))
+        if ( existing ) return existing
+
+        this.triangles.push(x)
+        return x
+    }
+
+    getFreePoints(): Point[] {
+        return this.points.filter(point => {
+            return !this.segments.some(segment => segment.hasPoint(point))
+        })
+    }
+
+    removeSegment(x: Segment) {
+        this.segments = this.segments.filter(segment => !segment.is(x))
+    }
+
+    clone() {
+        const newPoints: Point[] = this.points.map(point => point.clone())
+
+        const newSegments = this.segments.map(segment => {
+            const newFrom = newPoints.find(point => point.is(segment.from))
+            const newTo = newPoints.find(point => point.is(segment.to))
+            if ( !newFrom || !newTo ) {
+                throw new Error('Tried to clone segment, but could not match all points')
+            }
+
+            return new Segment(newFrom, newTo)
+        })
+
+        const newTriangles = this.triangles.map(triangle => {
+            const side0 = newSegments.find(segment => segment.is(triangle.sides[0]))
+            const side1 = newSegments.find(segment => segment.is(triangle.sides[1]))
+            const side2 = newSegments.find(segment => segment.is(triangle.sides[2]))
+            if ( !side0 || !side1 || !side2 ) {
+                throw new Error('Tried to clone triangle, but could not match all sides')
+            }
+
+            return new Triangle([side0, side1, side2])
+        })
+
+        return new Graph(newPoints, newSegments, newTriangles)
+    }
+}
+
+export function addBoundingSquareTo(graph: Graph): GraphBoundary {
+    const span = graph.span
+    const sideLength = 3 * span
+    const halfSideLength = sideLength / 2
+    const center = graph.center
+
+    const bottomLeftBound = Point.from(center.x - halfSideLength, center.y - halfSideLength)
+    const bottomRightBound = Point.from(center.x + halfSideLength, center.y - halfSideLength)
+    const topLeftBound = Point.from(center.x - halfSideLength, center.y + halfSideLength)
+    const topRightBound = Point.from(center.x + halfSideLength, center.y + halfSideLength)
+
+    const bottomLeft = graph.findExistingPointOrAdd(bottomLeftBound)
+    const bottomRight = graph.findExistingPointOrAdd(bottomRightBound)
+    const topLeft = graph.findExistingPointOrAdd(topLeftBound)
+    const topRight = graph.findExistingPointOrAdd(topRightBound)
+
+    const bottom = graph.findExistingSegmentOrAdd(new Segment(bottomLeft, bottomRight))
+    const right = graph.findExistingSegmentOrAdd(new Segment(bottomRight, topRight))
+    const left = graph.findExistingSegmentOrAdd(new Segment(bottomLeft, topLeft))
+    const top = graph.findExistingSegmentOrAdd(new Segment(topLeft, topRight))
+
+    return new GraphBoundary(
+        [bottomLeft, bottomRight, topRight, topLeft],
+        [bottom, right, top, left]
+    )
+}

src/test.ts

@@ -0,0 +1,36 @@
+import {deg2rad, Graph, Point, Segment, Triangle} from "./pslg";
+import {delaunay, delaunayRefine} from "./delaunayRefine";
+import {plotGraph} from "./viz";
+import {triangulate} from "./trapezoidTriangulation";
+
+const pa = Point.from(1, 0)
+const pb = Point.from(3, 3)
+const pc = Point.from(5, 2)
+const pd = Point.from(4, 0)
+
+const sA = new Segment(pa, pb)
+const sB = new Segment(pb, pc)
+const sC = new Segment(pc, pd)
+const sD = new Segment(pd, pa)
+
+const g = new Graph([pa, pb, pc, pd], [sA, sB, sC, sD])
+
+console.log(g)
+// plotGraph(g)
+
+const refined = triangulate(g)
+
+console.log(refined)
+plotGraph(refined)
+
+/*
+const p1 = Point.from(1, 1)
+const p2 = Point.from(5, 1)
+
+const pLeft = Point.from(3, 3)
+const pRight = Point.from(3, -3)
+
+const seg = new Segment(p1, p2)
+
+console.log('det of pLeft', seg.pointIsLeftOf(pLeft))
+console.log('det of pRight', seg.pointIsLeftOf(pRight))*/

src/trapezoidTriangulation.ts

@@ -0,0 +1,249 @@
+import {
+    addBoundingSquareTo,
+    getDirectionSorter,
+    Graph,
+    GraphBoundary,
+    GraphDirection,
+    Point,
+    Segment,
+    SegmentWithIntersection, Triangle
+} from "./pslg";
+
+export function getFirstIntersectingSegmentInDirection(raySegment: Segment, boundary: GraphBoundary, graph: Graph, direction: GraphDirection): [Segment, Point] {
+    const intersectingSegment = boundary.getBoundary(direction)
+    const intersectingPoint = raySegment.getIntersectionWith(intersectingSegment)
+    if ( !intersectingPoint ) {
+        console.log(
+            'getFirstIntersectingSegmentInDirection',
+            [raySegment.from.coordinate, raySegment.to.coordinate],
+            [intersectingSegment.from.coordinate, intersectingSegment.to.coordinate],
+            intersectingPoint
+        )
+        throw new RangeError('Ray segment does not extend to boundary in the given direction!')
+    }
+
+    // Now, collect all non-boundary segments that intersect with the ray segment
+    const intersections = (graph.segments
+        .map(segment => {
+            if ( boundary.isBoundarySegment(segment) ) return undefined
+
+            return {
+                segment,
+                intersect: segment.getIntersectionWithin(raySegment)
+            }
+        })
+        .filter(x => x && x.intersect) as SegmentWithIntersection[])
+        .sort(getDirectionSorter(direction))
+
+    const intersection = intersections[0]
+    if ( intersection ) {
+        return [intersection.segment, intersection.intersect]
+    }
+
+    return [intersectingSegment, intersectingPoint]
+}
+
+export function triangulate(originalGraph: Graph): Graph {
+    const graph = originalGraph.clone()
+    const boundary = addBoundingSquareTo(graph)
+
+    const leftBound = boundary.getLeftBoundary()
+    const rightBound = boundary.getRightBoundary()
+
+    const trapezoidSegments: Segment[] = []
+
+    // For each vertex in the original graph, create a horizontal line that
+    // extends in both directions until it intersects with either (1) the boundary
+    // or (2) a segment in the graph.
+    for ( const point of graph.points ) {
+        if ( boundary.isBoundaryPoint(point) ) continue  // skip boundary points
+
+        // Create the segment extending out to the left boundary
+        const leftPoint = Point.from(leftBound.from.x, point.y)
+        let leftSegment = new Segment(point, leftPoint)
+
+        // Get segments that intersect with this
+        const leftIntersectingSegments = (graph.segments
+            .map(segment => {
+                if ( boundary.isBoundarySegment(segment) ) {
+                    // Exclude boundary segments
+                    return undefined
+                }
+
+                return {
+                    segment,
+                    intersect: segment.getIntersectionWithin(leftSegment),
+                }
+            })
+            .filter(group => group && group.intersect) as Array<{segment: Segment, intersect: Point}>)
+            .sort((a, b) => {
+                return b.intersect.x - a.intersect.x  // Sort by right-most x-value
+            })
+
+        // Check if there was a nearer intersecting segment
+        const firstLeftIntersectingSegment = leftIntersectingSegments?.[0]?.segment
+        if ( firstLeftIntersectingSegment ) {
+            // Modify the leftSegment to end at the intersection point
+            const leftIntersect = graph.findExistingPointOrAdd(leftIntersectingSegments[0].intersect)
+            leftSegment = new Segment(point, leftIntersect)
+        }
+
+        // Create the segment extending out to the right boundary
+        const rightPoint = Point.from(rightBound.from.x, point.y)
+        let rightSegment = new Segment(point, rightPoint)
+
+        // Get segments that intersect with this
+        const rightIntersectingSegments = (graph.segments
+            .map(segment => {
+                if ( boundary.isBoundarySegment(segment) ) {
+                    // Exclude boundary segments
+                    return undefined
+                }
+
+                return {
+                    segment,
+                    intersect: segment.getIntersectionWithin(rightSegment),
+                }
+            })
+            .filter(group => group && group.intersect) as Array<{segment: Segment, intersect: Point}>)
+            .sort((a, b) => {
+                return a.intersect.x - b.intersect.x  // Sort by left-most x-value
+            })
+
+        // Check if there was a nearer intersecting segment
+        const firstRightIntersectingSegment = rightIntersectingSegments?.[0]?.segment
+        if ( firstRightIntersectingSegment ) {
+            // Modify the leftSegment to end at the intersection point
+            const rightIntersect = graph.findExistingPointOrAdd(rightIntersectingSegments[0].intersect)
+            rightSegment = new Segment(point, rightIntersect)
+        }
+
+        const graphLeftSegment = graph.findExistingSegmentOrAdd(leftSegment)
+        const graphRightSegment = graph.findExistingSegmentOrAdd(rightSegment)
+        trapezoidSegments.push(graphLeftSegment, graphRightSegment)
+    }
+
+    // Now, go through and identify trapezoids for all the horizontal segments we just added
+    for ( const segment of trapezoidSegments ) {
+        // First, find the trapezoid formed with the segment as the bottom
+        // Create a vertical segment from the midpoint of the segment to the top boundary
+        const horizontalMidpoint = segment.getMidpoint()
+        let upperBoundaryPoint = Point.from(horizontalMidpoint.x, boundary.ymax)
+        let upperBoundaryVerticalSegment = new Segment(horizontalMidpoint, upperBoundaryPoint)
+
+        const [upperIntersectSegment, upperIntersectPoint] = getFirstIntersectingSegmentInDirection(
+            upperBoundaryVerticalSegment,
+            boundary,
+            graph,
+            GraphDirection.UP
+        )
+
+        upperBoundaryVerticalSegment = new Segment(horizontalMidpoint, upperIntersectPoint)
+
+        // Now we have the upper and lower boundaries of the trapezoid.
+        // So, we need to figure out the left and right boundaries next.
+        // Get the midpoint of the vertical segment
+        const verticalMidpoint = upperBoundaryVerticalSegment.getMidpoint()
+        let leftBoundaryPoint = Point.from(boundary.xmin, verticalMidpoint.y)
+        let leftBoundaryHorizontalSegment = new Segment(verticalMidpoint, leftBoundaryPoint)
+
+        const [leftIntersectSegment, leftIntersectPoint] = getFirstIntersectingSegmentInDirection(
+            leftBoundaryHorizontalSegment,
+            boundary,
+            graph,
+            GraphDirection.LEFT
+        )
+
+        leftBoundaryHorizontalSegment = new Segment(verticalMidpoint, leftIntersectPoint)
+
+        // Repeat to get the right boundary
+        let rightBoundaryPoint = Point.from(boundary.xmax, verticalMidpoint.y)
+        let rightBoundaryHorizontalSegment = new Segment(verticalMidpoint, rightBoundaryPoint)
+
+        const [rightIntersectSegment, rightIntersectPoint] = getFirstIntersectingSegmentInDirection(
+            rightBoundaryHorizontalSegment,
+            boundary,
+            graph,
+            GraphDirection.RIGHT,
+        )
+
+        rightBoundaryHorizontalSegment = new Segment(verticalMidpoint, rightIntersectPoint)
+
+        // Now, check if we actually have a 4-bound trapezoid, or if we have a triangle
+        const points = Point.distinct([
+            segment.from,
+            segment.to,
+            upperIntersectSegment.from,
+            upperIntersectSegment.to,
+        ])
+
+        if ( points.length === 3 ) {
+            // We found a triangle! Less work.
+            // Create the triangle and push it onto the graph
+            const [p1, p2, p3] = points.map(x => graph.findExistingPointOrAdd(x))
+            const s12 = graph.findExistingSegmentOrAdd(new Segment(p1, p2))
+            const s23 = graph.findExistingSegmentOrAdd(new Segment(p2, p3))
+            const s31 = graph.findExistingSegmentOrAdd(new Segment(p3, p1))
+            graph.findExistingTriangleOrAdd(new Triangle([s12, s23, s31]))
+            continue  // FIXME - remove to handle below-segment case
+        }
+
+        if ( points.length !== 4 ) {
+            throw new RangeError('Found shape with invalid number of distinct points!')
+        }
+
+        // Now, we have the 4 bounding segments of the trapezoid.
+        // Let's find the segments that make up the trapezoid
+        // We will do this by re-creating segments for the four sides of the trapezoid
+        // Split the left-side on the intersection point
+        let [leftSegment1, leftSegment2] = leftIntersectSegment.splitAt(leftIntersectPoint)  // This is not right. Needs to be `segment`'s intersect point, not the midpoint intersect point
+        graph.removeSegment(leftIntersectSegment)
+        leftSegment1 = graph.findExistingSegmentOrAdd(leftSegment1)
+        leftSegment2 = graph.findExistingSegmentOrAdd(leftSegment2)
+
+        // We care about the upper-segment from the split, as that is the bound of our trapezoid
+        const trapezoidLeftBoundSegment = leftSegment1.ymin === leftIntersectPoint.y ? leftSegment1 : leftSegment2
+
+        // Repeat this process for the right-side segment
+        let [rightSegment1, rightSegment2] = rightIntersectSegment.splitAt(rightIntersectPoint)
+        graph.removeSegment(rightIntersectSegment)
+        rightSegment1 = graph.findExistingSegmentOrAdd(rightSegment1)
+        rightSegment2 = graph.findExistingSegmentOrAdd(rightSegment2)
+
+        const trapezoidRightBoundSegment = rightSegment1.ymin === rightBoundaryPoint.y ? leftSegment1 : leftSegment2
+
+        // Now we have all 4 bounding segments. We find the bisector that creates
+        // triangles with the largest minimum angle.
+
+        // First, try making triangles from bottom-left to top-right
+        const lowerLeftPoint = graph.findExistingPointOrAdd(Point.from(segment.xmin, segment.ymin))
+        const upperRightPoint = graph.findExistingPointOrAdd(Point.from(upperIntersectSegment.xmax, upperIntersectSegment.ymax))
+
+        const bottomLeftBisectorSegment = new Segment(lowerLeftPoint, upperRightPoint)
+        // const bottomLeftBisectorUpperTriangle = new Triangle([bottomLeftBisectorSegment, upperIntersectSegment, trapezoidLeftBoundSegment])
+        // const bottomLeftBisectorLowerTriangle = new Triangle([bottomLeftBisectorSegment, segment, trapezoidRightBoundSegment])
+        // const bottomLeftBisectorMinAngle = Math.min(bottomLeftBisectorUpperTriangle.getMinimumAngle(), bottomLeftBisectorLowerTriangle.getMinimumAngle())
+
+        const upperLeftPoint = graph.findExistingPointOrAdd(Point.from(upperIntersectSegment.xmin, upperIntersectSegment.ymax))
+        const lowerRightPoint = graph.findExistingPointOrAdd(Point.from(segment.xmax, segment.ymin))
+
+        // const topRightBisectorSegment = new Segment(upperLeftPoint, lowerRightPoint)
+        // const upperRightBisectorUpperTriangle = new Triangle([topRightBisectorSegment, upperIntersectSegment, trapezoidRightBoundSegment])
+        // const upperRightBisectorLowerTriangle = new Triangle([topRightBisectorSegment, trapezoidLeftBoundSegment, segment])
+        // const upperRightBisectorMinAngle = Math.min(upperRightBisectorUpperTriangle.getMinimumAngle(), upperRightBisectorLowerTriangle.getMinimumAngle())
+
+        // const optimalBisectorUpperTriangle = upperRightBisectorMinAngle > bottomLeftBisectorMinAngle ? upperRightBisectorUpperTriangle : bottomLeftBisectorUpperTriangle
+        // const optimalBisectorLowerTriangle = upperRightBisectorMinAngle > bottomLeftBisectorMinAngle ? upperRightBisectorLowerTriangle : bottomLeftBisectorLowerTriangle
+
+        // Add the triangles to the graph
+        // const upperTriangleSegments = optimalBisectorUpperTriangle.sides.map(side => graph.findExistingSegmentOrAdd(side))
+        // graph.findExistingTriangleOrAdd(new Triangle(upperTriangleSegments as [Segment, Segment, Segment]))
+
+        // const lowerTriangleSegments = optimalBisectorLowerTriangle.sides.map(side => graph.findExistingSegmentOrAdd(side))
+        // graph.findExistingTriangleOrAdd(new Triangle(lowerTriangleSegments as [Segment, Segment, Segment]))
+    }
+
+    // FIXME handle the lower-trapezoid case
+
+    return graph
+}

src/viz.ts

@@ -0,0 +1,25 @@
+import {Graph} from "./pslg";
+import {Plot, plot} from "nodeplotlib";
+
+export function plotGraph(graph: Graph) {
+    const segmentTraces = graph.segments.map((segment): Plot => {
+        return {
+            x: [segment.from.x, segment.to.x],
+            y: [segment.from.y, segment.to.y],
+            type: 'scatter',
+        }
+    })
+
+    const traces = segmentTraces.concat(
+        graph.getFreePoints()
+            .map((point): Plot => {
+                return {
+                    x: [point.x],
+                    y: [point.y],
+                    type: 'scatter',
+                }
+            })
+    )
+
+    plot(traces)
+}

tsconfig.json

@@ -0,0 +1,11 @@
+{
+  "compilerOptions": {
+    "target": "es6",
+    "module": "CommonJS",
+    "declaration": true,
+    "outDir": "./lib",
+    "strict": true,
+  },
+  "include": ["src"],
+  "exclude": ["node_modules"]
+}