diff --git a/examples/website/hierarchy/app.tsx b/examples/website/hierarchy/app.tsx index 4651cf39..4f89cda3 100644 --- a/examples/website/hierarchy/app.tsx +++ b/examples/website/hierarchy/app.tsx @@ -34,7 +34,7 @@ const App: React.FC<{showCellId?: boolean}> = ({showCellId = true}) => { }, []); // Calculate resolution based on zoom level - let resolution = Math.min(Math.floor(2 * viewState.zoom - 4), Math.floor(viewState.zoom + 3)); + let resolution = Math.min(Math.floor(2 * viewState.zoom - 4), Math.floor(viewState.zoom + 1)); resolution = Math.max(1, Math.min(MAX_RESOLUTION, resolution)); // Memoize the entire cells calculation diff --git a/examples/website/spherical-polygon/app.tsx b/examples/website/spherical-polygon/app.tsx new file mode 100644 index 00000000..3b0d093a --- /dev/null +++ b/examples/website/spherical-polygon/app.tsx @@ -0,0 +1,108 @@ +import React, { Suspense, useState } from 'react'; +import { createRoot } from 'react-dom/client'; +import { Canvas } from '@react-three/fiber'; +import { OrbitControls } from '@react-three/drei'; +import { hexToBigInt, cellToChildren, cellToParent } from 'a5/index'; +import { A5Pentagon, Marker, Sphere, sphericalPentagonFromCell } from './components'; +import { toCartesian } from 'a5/core/coordinate-transforms'; +import { Cartesian, Spherical } from 'a5/core/coordinate-systems'; +import { vec3 } from 'gl-matrix'; + +const initialSpherical = [2 * Math.PI * Math.random(), Math.PI * Math.random()] as Spherical; +const initialPoint = toCartesian(initialSpherical); +const camera = toCartesian(initialSpherical); +vec3.scale(camera, camera, 2); + +const a5CellHex = '1200000000000000'; +const a5cell = hexToBigInt(a5CellHex); + +function Scene({ resolution }: { resolution: number }) { + const [point, setPoint] = useState(initialPoint); + + // Get cells at current resolution + const a5cells = cellToChildren(cellToParent(a5cell), resolution); + + // Filter cells based on current point + const filteredCells = a5cells.filter(cell => { + const pentagon = sphericalPentagonFromCell(cell); + return pentagon.containsPoint(point) > 0; + }); + + const handleSphereClick = (intersection: Spherical) => { + setPoint(toCartesian(intersection)); + }; + + return ( + <> + + + + + {filteredCells.map(cell => )} + {a5cells.map(cell => )} + + + + ); +} + +const App: React.FC = () => { + const [resolution, setResolution] = useState(3); + + return ( +
+
+
+ Resolution: {resolution} +
+ setResolution(Number(e.target.value))} + style={{ + width: '100%', + height: '20px', + WebkitAppearance: 'none', + background: 'rgba(135, 206, 235, 0.2)', + borderRadius: '10px', + outline: 'none' + }} + /> +
+ + + + + + +
+ ); +}; + +export default App; + +export async function renderToDOM(container: HTMLDivElement) { + const root = createRoot(container); + root.render(); +} \ No newline at end of file diff --git a/examples/website/spherical-polygon/components.tsx b/examples/website/spherical-polygon/components.tsx new file mode 100644 index 00000000..32e17ece --- /dev/null +++ b/examples/website/spherical-polygon/components.tsx @@ -0,0 +1,73 @@ +import React, { useMemo } from 'react'; +import { BufferGeometry, Float32BufferAttribute, DoubleSide, Vector3, Raycaster, Sphere as ThreeSphere } from 'three'; +import { SphericalPentagonShape } from 'a5/core/spherical-pentagon'; +import { toCartesian, fromLonLat, toSpherical } from 'a5/core/coordinate-transforms'; +import { cellToBoundary } from 'a5/index'; +import type { Spherical, Radians, Cartesian } from 'a5/core/coordinate-systems'; +import { useThree } from '@react-three/fiber'; + +export function Sphere({ onSphereClick }: { onSphereClick: (point: Spherical) => void }) { + const { camera} = useThree(); + + const handleClick = (event: { clientX: number; clientY: number }) => { + const spherical = toSpherical(camera.position.toArray() as Cartesian); + onSphereClick(spherical); + }; + + return ( + + + + + ); +} + +export function PentagonLines(props: { sphericalPentagon: SphericalPentagonShape, disabled: boolean }) { + const { sphericalPentagon, disabled } = props; + const geometry = useMemo(() => { + // Use 20 segments per edge for smooth curves + const vertices = sphericalPentagon.getBoundary(20); + const geometry = new BufferGeometry(); + geometry.setAttribute('position', new Float32BufferAttribute(vertices.flatMap(p => [...p]), 3)); + return geometry; + }, [sphericalPentagon]); + + return ( + + + + ); +} + +export function sphericalPentagonFromCell(cell: bigint): SphericalPentagonShape { + const boundary = cellToBoundary(cell); + const cartesianBoundary = boundary.map(lonlat => toCartesian(fromLonLat(lonlat))); + return new SphericalPentagonShape(cartesianBoundary); +} + +export function A5Pentagon(props: { cell: bigint, disabled: boolean }) { + const {cell, disabled} = props; + const a5Pentagon = sphericalPentagonFromCell(cell); + return ; +} + +export function Marker(props: { cartesian: Cartesian }) { + const cartesian = props.cartesian; + return ( + + + + + ); +} \ No newline at end of file diff --git a/modules/core/cell.ts b/modules/core/cell.ts index 497d5cc4..78fe2866 100644 --- a/modules/core/cell.ts +++ b/modules/core/cell.ts @@ -6,15 +6,16 @@ import { mat2, vec2, glMatrix } from "gl-matrix"; glMatrix.setMatrixArrayType(Float64Array as any); import type { Face, LonLat } from "./coordinate-systems"; -import { FaceToIJ, fromLonLat, toFace } from "./coordinate-transforms"; +import { FaceToIJ, fromLonLat, toCartesian, toFace } from "./coordinate-transforms"; import { findNearestOrigin, quintantToSegment, segmentToQuintant } from "./origin"; import { unprojectDodecahedron } from "./dodecahedron"; -import { A5Cell, Pentagon, PentagonShape } from "./utils"; -import { getFaceVertices, getPentagonVertices, getQuintant, getQuintantPolar, getQuintantVertices } from "./tiling"; +import { A5Cell, PentagonShape } from "./utils"; +import { getFaceVertices, getPentagonVertices, getQuintantPolar, getQuintantVertices } from "./tiling"; import { PI_OVER_5 } from "./constants"; import { IJToS, sToAnchor } from "./hilbert"; -import { projectPentagon, projectPoint, reprojectPentagon } from "./project"; +import { projectPentagon, projectPoint } from "./project"; import { deserialize, serialize, FIRST_HILBERT_RESOLUTION } from "./serialization"; +import { SphericalPentagonShape } from "./spherical-pentagon"; // Reuse these objects to avoid allocation const rotation = mat2.create(); @@ -27,7 +28,7 @@ export function lonLatToCell(lonLat: LonLat, resolution: number): bigint { const hilbertResolution = 1 + resolution - FIRST_HILBERT_RESOLUTION; const samples: LonLat[] = [lonLat]; - const N = 100; + const N = 25; const scale = 50 / Math.pow(2, hilbertResolution); for (let i = 0; i < N; i++) { const R = (i / N) * scale; @@ -36,35 +37,31 @@ export function lonLatToCell(lonLat: LonLat, resolution: number): bigint { samples.push(coordinate as LonLat); } - const cells: {cell: A5Cell, distance: number}[] = []; - const estimates: {estimate: A5Cell, sample: LonLat}[] = []; - for (const sample of samples) { - const estimate = _lonLatToEstimate(sample, resolution); - estimates.push({estimate, sample}); - } - // Deduplicate estimates const estimateSet = new Set(); const uniqueEstimates: A5Cell[] = []; - for (const {estimate, sample} of estimates) { + + const cells: {cell: A5Cell, distance: number}[] = []; + for (const sample of samples) { + const estimate = _lonLatToEstimate(sample, resolution); const estimateKey = serialize(estimate); if (!estimateSet.has(estimateKey)) { + // Have new estimate, add to set and list estimateSet.add(estimateKey); uniqueEstimates.push(estimate); - } - } - for (const estimate of uniqueEstimates) { - const distance = a5cellContainsPoint(estimate, lonLat); - if (distance < 0) { - return serialize(estimate); - } else { - cells.push({cell: estimate, distance}); + // Check if we have a hit, storing distance if not + const distance = a5cellContainsPoint(estimate, lonLat); + if (distance > 0) { + return serialize(estimate); + } else { + cells.push({cell: estimate, distance}); + } } } - // Sort cells by distance and use the closest one - cells.sort((a, b) => a.distance - b.distance); + // As fallback, sort cells by distance and use the closest one + cells.sort((a, b) => b.distance - a.distance); return serialize(cells[0].cell); } @@ -129,17 +126,11 @@ export function cellToBoundary(cellId: bigint): LonLat[] { } export function a5cellContainsPoint(cell: A5Cell, point: LonLat): number { - const spherical = fromLonLat(point); - - // Important to use the same origin as the cell, so we unproject onto correct face - const {origin} = cell; - const polar = unprojectDodecahedron(spherical, origin.quat, origin.angle); - const face = toFace(polar); + const boundary = cellToBoundary(serialize(cell)); - // Required for points on pentagon that cross the origin boundary - const pentagon = _getPentagon(cell); - const reprojectedPentagon = reprojectPentagon(pentagon, origin); + const cartesian = toCartesian(fromLonLat(point)); + const sphericalBoundary = boundary.map(vertex => toCartesian(fromLonLat(vertex))); - // Perform containment test in Face coordinates, where cell edges are straight lines - return reprojectedPentagon.containsPoint(face); + const sphericalPentagon = new SphericalPentagonShape(sphericalBoundary); + return sphericalPentagon.containsPoint(cartesian); } \ No newline at end of file diff --git a/modules/core/coordinate-transforms.ts b/modules/core/coordinate-transforms.ts index aae0e659..85e017d4 100644 --- a/modules/core/coordinate-transforms.ts +++ b/modules/core/coordinate-transforms.ts @@ -42,13 +42,13 @@ export function toSpherical(xyz: Cartesian): Spherical { } export function toCartesian([theta, phi]: Spherical): Cartesian { - const x = Math.sin(phi) * Math.cos(theta); - const y = Math.sin(phi) * Math.sin(theta); + const sinPhi = Math.sin(phi); + const x = sinPhi * Math.cos(theta); + const y = sinPhi * Math.sin(theta); const z = Math.cos(phi); return [x, y, z] as Cartesian; } - /** * Determine the offset longitude for the spherical coordinate system * This is the angle between the Greenwich meridian and vector between the centers @@ -83,6 +83,12 @@ export function toLonLat([theta, phi]: Spherical): LonLat { return [longitude, latitude] as LonLat; } +/** + * Creates a quaternion representing a rotation + * from the north pole to a given axis. + * @param axis Spherical coordinate of axis to rotate to + * @returns quaternion + */ export function quatFromSpherical(axis: Spherical): quat { const cartesian = toCartesian(axis); const Q = quat.create(); diff --git a/modules/core/project.ts b/modules/core/project.ts index d1e701bf..7806b079 100644 --- a/modules/core/project.ts +++ b/modules/core/project.ts @@ -4,13 +4,13 @@ import { vec2, mat2, glMatrix } from 'gl-matrix'; glMatrix.setMatrixArrayType(Float64Array as any); -import { Pentagon, PentagonShape } from './utils'; +import { PentagonShape } from './utils'; import { Origin } from './utils'; import { movePointToFace, findNearestOrigin, isNearestOrigin } from './origin'; -import { projectDodecahedron, unprojectDodecahedron } from './dodecahedron'; -import type { Face, LonLat, Polar, Radians } from './coordinate-systems'; -import { fromLonLat, toFace, toLonLat, toPolar } from './coordinate-transforms'; -import { distanceToEdge, PI_OVER_5 } from './constants'; +import { projectDodecahedron } from './dodecahedron'; +import type { Face, LonLat, Radians } from './coordinate-systems'; +import { toLonLat, toPolar } from './coordinate-transforms'; +import { PI_OVER_5 } from './constants'; // Reusable matrices to avoid recreation const rotation = mat2.create(); @@ -51,38 +51,4 @@ export function projectPentagon(pentagon: PentagonShape, origin: Origin): LonLat // Normalize longitudes to handle antimeridian crossing const normalizedVertices = PentagonShape.normalizeLongitudes(rotatedVertices); return normalizedVertices; -} - -/** - * Reproject a pentagon so that all vertices are defined in the same Face coordinate system - * of a single origin. This is required for containment testing as the vertices of a cell - * can span multiple origins, which are warped differently. Applying this function allows - * us to perform containment testing in Face coordinates, where cell edges are straight lines. - * @param pentagon - * @param origin - * @returns - */ -export function reprojectPentagon(pentagon: PentagonShape, origin: Origin): PentagonShape { - // If all vertices are within a single origin, we can just return the pentagon - let withinSingleOrigin = true; - for (const vertex of pentagon.getVertices()) { - const D = vec2.length(vertex); - if (D > distanceToEdge) { - withinSingleOrigin = false; - break; - } - } - - if (withinSingleOrigin) { - return pentagon; - } - - // Project to sphere, projectPoint will handle "folding" vertices across the dodecahedron edges - const projectedPentagon = projectPentagon(pentagon, origin); - const sphericalPentagon = projectedPentagon.map(fromLonLat); - - // Unproject to dodecahedron, with all vertices now defined in the same Face coordinate system - const polarPentagon = sphericalPentagon.map(spherical => unprojectDodecahedron(spherical, origin.quat, origin.angle)); - const facePentagon = polarPentagon.map(polar => toFace(polar)) as Pentagon; - return new PentagonShape(facePentagon); } \ No newline at end of file diff --git a/modules/core/spherical-pentagon.ts b/modules/core/spherical-pentagon.ts new file mode 100644 index 00000000..94bf7a0c --- /dev/null +++ b/modules/core/spherical-pentagon.ts @@ -0,0 +1,133 @@ +// A5 +// SPDX-License-Identifier: Apache-2.0 +// Copyright (c) A5 contributors + +import {vec3, glMatrix, quat} from 'gl-matrix'; +glMatrix.setMatrixArrayType(Float64Array as any); +import type { Cartesian } from './coordinate-systems'; + +// Use Cartesian system for all calculations for greater accuracy +// Using [x, y, z] gives equal precision in all directions, unlike spherical coordinates +export type SphericalPolygon = Cartesian[]; +const UP = [0, 0, 1] as Cartesian; + +export class SphericalPentagonShape { + private vertices: SphericalPolygon; + + constructor(vertices: SphericalPolygon) { + this.vertices = vertices; + if (!this.isWindingCorrect()) { + this.vertices.reverse(); + } + } + + /** + * + * @param nSegments Returns a close boundary of the polygon, with nSegments points per edge + * @returns SphericalPolygon + */ + getBoundary(nSegments: number = 1): SphericalPolygon { + const points: SphericalPolygon = []; + const N = this.vertices.length; + for (let s = 0; s < N * nSegments; s++) { + const t = s / nSegments; + points.push(this.slerp(t)); + } + points.push(points[0]); + + return points; + } + + /** + * Interpolates along boundary of polygon. Pass t = 1.5 to get the midpoint between 2nd and 3rd vertices + * @param t + * @returns Cartesian coordinate + */ + slerp(t: number): Cartesian { + const N = this.vertices.length; + const f = t % 1; + const i = Math.floor(t % N); + const j = (i + 1) % N; + + // Points A & B + const A = this.vertices[i]; + const B = this.vertices[j]; + + // Quaternions + const identity = quat.create(); + const qOA = quat.rotationTo(quat.create(), UP, A); + const qAB = quat.rotationTo(quat.create(), A, B); + const qPartial = quat.slerp(quat.create(), identity, qAB, f); + const qCombined = quat.multiply(quat.create(), qPartial, qOA); + + const out = vec3.fromValues(0, 0, 1); + vec3.transformQuat(out, out, qCombined); + return out as Cartesian; + } + + /** + * Returns the vertex given by index t, along with the vectors: + * - VA: Vector from vertex to point A + * - VB: Vector from vertex to point B + * @param t + * @returns + */ + getTransformedVertices(t: number): [Cartesian, Cartesian, Cartesian] { + const N = this.vertices.length; + const i = Math.floor(t % N); + const j = (i + 1) % N; + const k = (i + N - 1) % N; + + // Points A & B (vertex before and after) + const V = vec3.clone(this.vertices[i]) as Cartesian; + const VA = vec3.clone(this.vertices[j]) as Cartesian; + const VB = vec3.clone(this.vertices[k]) as Cartesian; + vec3.sub(VA, VA, V); + vec3.sub(VB, VB, V); + return [V, VA, VB]; + } + + containsPoint(point: Cartesian): number { + // Adaption of algorithm from: + // 'Locating a point on a spherical surface relative to a spherical polygon' + // Using only the condition of 'necessary strike' + const N = this.vertices.length; + let thetaDeltaMin = Infinity; + + for (let i = 0; i < N; i++) { + // Transform point and neighboring vertices into coordinate system centered on vertex + const [V, VA, VB] = this.getTransformedVertices(i); + const VP = vec3.sub(vec3.create(), point, V); + + // Normalize to obtain unit direction vectors + vec3.normalize(VP, VP); + vec3.normalize(VA, VA); + vec3.normalize(VB, VB); + + // Cross products will point away from the center of the sphere when + // point P is within arc formed by VA and VB + const crossPA = vec3.cross(vec3.create(), VP, VA); + const crossBP = vec3.cross(vec3.create(), VB, VP); + + // Dot product will be positive when point P is within arc formed by VA and VB + // The magnitude of the dot product is the sine of the angle between the two vectors + // which is the same as the angle for small angles. + const sinPA = vec3.dot(V, crossPA); + const sinBP = vec3.dot(V, crossBP); + + // By returning the minimum value we find the arc where the point is closest to being outside + thetaDeltaMin = Math.min(thetaDeltaMin, sinPA, sinBP); + } + + // If point is inside all arcs, will return a position value + // If point is on edge of arc, will return 0 + // If point is outside all arcs, will return -1, the further away from 0, the further away from the arc + return thetaDeltaMin; + } + + private isWindingCorrect(): boolean { + const [V, VA, VB] = this.getTransformedVertices(0); + const cross = vec3.cross(vec3.create(), VA, VB); + return vec3.dot(V, cross) <= 0; + } +} \ No newline at end of file diff --git a/modules/core/tiling.ts b/modules/core/tiling.ts index cf792e76..ba2112f8 100644 --- a/modules/core/tiling.ts +++ b/modules/core/tiling.ts @@ -89,14 +89,6 @@ export function getFaceVertices(): PentagonShape { return new PentagonShape(vertices as Pentagon); } - -export function getQuintant(point: vec2): number { - // TODO perhaps quicker way without trigonometry - const angle = Math.atan2(point[1], point[0]); - const normalizedAngle = (angle - V + TWO_PI) % TWO_PI; - return Math.ceil(normalizedAngle / TWO_PI_OVER_5) % 5; -} - export function getQuintantPolar([_, gamma]: Polar): number { return (Math.round(gamma / TWO_PI_OVER_5) + 5) % 5; } \ No newline at end of file diff --git a/tests/cell.test.ts b/tests/cell.test.ts index 29efb3d7..f1da7b5a 100644 --- a/tests/cell.test.ts +++ b/tests/cell.test.ts @@ -120,6 +120,9 @@ describe('Cell Boundary Tests', () => { // Store failures for this resolution if any occurred if (resolutionFailures.length > 0) { + if (!failures[pointKey]) { + failures[pointKey] = {}; + } failures[pointKey][resolution] = resolutionFailures; } } diff --git a/tests/coordinate-transforms.test.ts b/tests/coordinate-transforms.test.ts index 637feba9..fc321788 100644 --- a/tests/coordinate-transforms.test.ts +++ b/tests/coordinate-transforms.test.ts @@ -9,6 +9,20 @@ import { } from 'a5/core/coordinate-transforms' import type { Degrees, LonLat, Radians, Spherical } from 'a5/core/coordinate-systems' +const TEST_POINTS: Array = [ + [0, 0], // Equator + [90, 0], // Equator + [180, 0], // Equator + [0, 45], // Mid latitude + [0, -45], // Mid latitude + [-90, -45], // West hemisphere mid-latitude + [180, 45], // Date line mid-latitude + [90, 45], // East hemisphere mid-latitude + [0, 90], // North pole + [0, -90], // South pole + [123, 45], // Arbitrary point +] as LonLat[]; + describe('angle conversions', () => { it('converts degrees to radians', () => { expect(degToRad(180 as Degrees)).toBe(Math.PI) @@ -74,14 +88,6 @@ describe('LonLat to/from spherical', () => { it('converts spherical to longitude/latitude coordinates', () => { // Test round trip conversion - const TEST_POINTS: Array<[number, number]> = [ - [0, 0], // Greenwich equator - [0, 90], // North pole - [0, -90], // South pole - [180, 45], // Date line mid-latitude - [-90, -45], // West hemisphere mid-latitude - ]; - TEST_POINTS.forEach(([lon, lat]) => { const spherical = fromLonLat([lon, lat] as LonLat) const [newLon, newLat] = toLonLat(spherical) @@ -90,4 +96,4 @@ describe('LonLat to/from spherical', () => { expect(newLat).toBeCloseTo(lat) }) }) -}) \ No newline at end of file +}); \ No newline at end of file diff --git a/website/src/examples-sidebar.js b/website/src/examples-sidebar.js index edae6f3a..d2210ad4 100644 --- a/website/src/examples-sidebar.js +++ b/website/src/examples-sidebar.js @@ -13,7 +13,7 @@ const sidebars = { { type: 'category', label: 'Technical', - items: ['pentagon', 'teohedron-dodecahedron', 'lattice', 'hilbert', 'globe'] + items: ['pentagon', 'teohedron-dodecahedron', 'spherical-polygon', 'lattice', 'hilbert', 'globe'] }, { type: 'category', diff --git a/website/src/examples/spherical-polygon.js b/website/src/examples/spherical-polygon.js new file mode 100644 index 00000000..727f494e --- /dev/null +++ b/website/src/examples/spherical-polygon.js @@ -0,0 +1,32 @@ +import React, {Component} from 'react'; +import {GITHUB_TREE} from '../constants/defaults'; +import App from 'website-examples/spherical-polygon/app'; + +import {makeExample} from '../components'; + +class SphericalPolygonDemo extends Component { + static title = 'Spherical Polygon'; + + static code = `${GITHUB_TREE}/examples/website/spherical-polygon`; + + static parameters = {}; + + static renderInfo(meta) { + return ( +
+

Visualization of a hit-testing of spherical polygons on the globe.

+

Rotate the globe to move the hit-testing point.

+
+ ); + } + + render() { + return ( +
+ +
+ ); + } +} + +export default makeExample(SphericalPolygonDemo); \ No newline at end of file diff --git a/website/src/examples/spherical-polygon.mdx b/website/src/examples/spherical-polygon.mdx new file mode 100644 index 00000000..204d90f7 --- /dev/null +++ b/website/src/examples/spherical-polygon.mdx @@ -0,0 +1,5 @@ +# Spherical Polygon + +import Demo from './spherical-polygon'; + + \ No newline at end of file