From 35b7f33097b3e53009f3988fa1d75b6fa29f8f37 Mon Sep 17 00:00:00 2001 From: garrettmills Date: Fri, 24 Sep 2021 03:20:39 -0500 Subject: [PATCH] Initial import --- .gitignore | 3 + Makefile | 9 + package.json | 5 + pnpm-lock.yaml | 25 ++ src/bounding.ts | 0 src/delaunayRefine.ts | 92 ++++ src/linear.ts | 58 +++ src/pslg.ts | 766 ++++++++++++++++++++++++++++++++++ src/test.ts | 36 ++ src/trapezoidTriangulation.ts | 249 +++++++++++ src/viz.ts | 25 ++ tsconfig.json | 11 + 12 files changed, 1279 insertions(+) create mode 100644 .gitignore create mode 100644 Makefile create mode 100644 package.json create mode 100644 pnpm-lock.yaml create mode 100644 src/bounding.ts create mode 100644 src/delaunayRefine.ts create mode 100644 src/linear.ts create mode 100644 src/pslg.ts create mode 100644 src/test.ts create mode 100644 src/trapezoidTriangulation.ts create mode 100644 src/viz.ts create mode 100644 tsconfig.json diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..49c3a86 --- /dev/null +++ b/.gitignore @@ -0,0 +1,3 @@ +.idea* +lib* +node_modules* diff --git a/Makefile b/Makefile new file mode 100644 index 0000000..fa7220f --- /dev/null +++ b/Makefile @@ -0,0 +1,9 @@ + +.PHONY: all +all: + rm -rf lib + tsc + +.PHONY: run +run: all + node lib/test.js diff --git a/package.json b/package.json new file mode 100644 index 0000000..d884c2b --- /dev/null +++ b/package.json @@ -0,0 +1,5 @@ +{ + "dependencies": { + "nodeplotlib": "^0.7.5" + } +} diff --git a/pnpm-lock.yaml b/pnpm-lock.yaml new file mode 100644 index 0000000..22a6066 --- /dev/null +++ b/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 diff --git a/src/bounding.ts b/src/bounding.ts new file mode 100644 index 0000000..e69de29 diff --git a/src/delaunayRefine.ts b/src/delaunayRefine.ts new file mode 100644 index 0000000..81d907d --- /dev/null +++ b/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 ( { + 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 +} diff --git a/src/linear.ts b/src/linear.ts new file mode 100644 index 0000000..67e9fa4 --- /dev/null +++ b/src/linear.ts @@ -0,0 +1,58 @@ +type TupleOf = N extends N ? number extends N ? T[] : _TupleOf : never +type _TupleOf = R['length'] extends N ? R : _TupleOf + +export class Vector { + public readonly values: TupleOf + + constructor( + public readonly length: N + ) { + this.values = Array(length).fill(0) as TupleOf + } +} + +export class Matrix { + public readonly values: TupleOf, R> + + constructor( + public readonly rows: R, + public readonly cols: C, + ) { + this.values = Array(rows).fill(0).map(() => Array(cols).fill(0) as TupleOf) as TupleOf, 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') + } + } +} diff --git a/src/pslg.ts b/src/pslg.ts new file mode 100644 index 0000000..8ea8234 --- /dev/null +++ b/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] + ) +} diff --git a/src/test.ts b/src/test.ts new file mode 100644 index 0000000..b8f4b33 --- /dev/null +++ b/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))*/ diff --git a/src/trapezoidTriangulation.ts b/src/trapezoidTriangulation.ts new file mode 100644 index 0000000..5dd1538 --- /dev/null +++ b/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 +} diff --git a/src/viz.ts b/src/viz.ts new file mode 100644 index 0000000..d672703 --- /dev/null +++ b/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) +} diff --git a/tsconfig.json b/tsconfig.json new file mode 100644 index 0000000..922f4ef --- /dev/null +++ b/tsconfig.json @@ -0,0 +1,11 @@ +{ + "compilerOptions": { + "target": "es6", + "module": "CommonJS", + "declaration": true, + "outDir": "./lib", + "strict": true, + }, + "include": ["src"], + "exclude": ["node_modules"] +}