diff --git a/src/elevation.ts b/src/elevation.ts new file mode 100644 index 0000000..7322f8a --- /dev/null +++ b/src/elevation.ts @@ -0,0 +1,87 @@ +export type Point = { + x: number; + y: number; + z: number; +}; + +export type SphericalPoint = { + altitude: number; + azimuth: number; +}; + +export function fillMissingAltitudes(maxAngles: SphericalPoint[]): void { + // First copy the maxAngles to a newAngles list, so that changes + // in the list do not affect the algorithm + let newAngles = maxAngles.map((angle) => ({ ...angle })); + for (let i = 0; i < newAngles.length; i++) { + if (newAngles[i].altitude != -Infinity) { + continue; + } + let distance = 1; + while (true) { + let prevIndex = (i - distance + newAngles.length) % newAngles.length; + let nextIndex = (i + distance) % newAngles.length; + + if (maxAngles[nextIndex].altitude !== -Infinity) { + newAngles[i].altitude = maxAngles[nextIndex].altitude; + break; + } else if (maxAngles[prevIndex].altitude !== -Infinity) { + newAngles[i].altitude = maxAngles[prevIndex].altitude; + break; + } else distance++; + } + } + // Overwrite the maxAngles to make changes in this vector global + for (let i = 0; i < maxAngles.length; i++) { + maxAngles[i] = newAngles[i]; + } +} + +/** + * + * @param start + * @param end + * @returns azimuth from 0 to 2*PI and altitude from 0 to PI/2, where altitude = 0 is facing directly upwards + */ +export function calculateSphericalCoordinates(start: Point, end: Point): { azimuth: number; altitude: number } { + const dx = end.x - start.x; + const dy = end.y - start.y; + const dz = end.z - start.z; + + const r = Math.sqrt(dx * dx + dy * dy + dz * dz); + const altitude = Math.acos(dz / r); + let azimuth = Math.atan2(dy, dx); + + if (azimuth < 0) { + azimuth += 2 * Math.PI; // Adjust azimuth to be from 0 to 2PI + } + + return { azimuth, altitude }; +} + +/** + * Calculates the maximum heights visible from an observer in a set of directions. + * Returns a list of spherical points of length numDirections. + * @param elevation list of points with x,y,z component + * @param observer Point of interest for which the elevation angles are calculated. + * @param numDirections Number of steps for the azimuth angle. + * @returns + */ +export function getMaxElevationAngles(elevation: Point[], observer: Point, numDirections: number = 360): SphericalPoint[] { + let maxAngles: SphericalPoint[] = Array.from({ length: numDirections }, (_, index) => ({ + azimuth: index * ((2 * Math.PI) / numDirections), + altitude: -Infinity, + })); + + for (let point of elevation) { + const { azimuth, altitude } = calculateSphericalCoordinates(observer, point); + console.log(azimuth, altitude); + const closestIndex = Math.round(azimuth / ((2 * Math.PI) / numDirections)) % numDirections; + + if (altitude > maxAngles[closestIndex].altitude) { + maxAngles[closestIndex].altitude = altitude; + } + } + fillMissingAltitudes(maxAngles); + return maxAngles; +} diff --git a/src/index.ts b/src/index.ts index 8247b0a..c65bea1 100644 --- a/src/index.ts +++ b/src/index.ts @@ -2,9 +2,9 @@ import * as THREE from 'three'; import { BufferAttribute, BufferGeometry, TypedArray } from 'three'; import * as BufferGeometryUtils from 'three/examples/jsm/utils/BufferGeometryUtils.js'; import { viridis } from './colormaps'; +import * as elevation from './elevation'; import { getRandomSunVectors } from './sun'; import * as triangleUtils from './triangleUtils.js'; -import { ArrayType, Triangle } from './triangleUtils.js'; // @ts-ignore import { rayTracingWebGL } from './rayTracingWebGL.js'; @@ -20,6 +20,8 @@ import { rayTracingWebGL } from './rayTracingWebGL.js'; export default class Scene { simulationGeometries: Array; shadingGeometries: Array; + elevationRaster: Array; + elevationRasterMidpoint: elevation.Point; latitude: number; longitude: number; @@ -31,6 +33,8 @@ export default class Scene { constructor(latitude: number, longitude: number) { this.simulationGeometries = []; this.shadingGeometries = []; + this.elevationRaster = []; + this.elevationRasterMidpoint = { x: 0, y: 0, z: 0 }; this.latitude = latitude; this.longitude = longitude; } @@ -57,6 +61,11 @@ export default class Scene { this.shadingGeometries.push(geometry); } + addElevationRaster(raster: elevation.Point[], midpoint: elevation.Point) { + this.elevationRaster = raster; + this.elevationRasterMidpoint = midpoint; + } + /** @ignore */ refineMesh(mesh: BufferGeometry, maxLength: number): BufferGeometry { const positions = mesh.attributes.position.array.slice(); @@ -102,6 +111,7 @@ export default class Scene { console.log('Simulation package was called to calculate'); let simulationGeometry = BufferGeometryUtils.mergeGeometries(this.simulationGeometries); let shadingGeometry = BufferGeometryUtils.mergeGeometries(this.shadingGeometries); + // TODO: This breaks everything, why? simulationGeometry = this.refineMesh(simulationGeometry, 0.5); // TODO: make configurable @@ -187,6 +197,7 @@ export default class Scene { * @param midpoints midpoints of triangles for which to calculate intensities * @param normals normals for each midpoint * @param meshArray array of vertices for the shading mesh + * @param numberSimulations * @return * @memberof Scene */ @@ -198,6 +209,15 @@ export default class Scene { progressCallback: (progress: number, total: number) => void, ) { let sunDirections = getRandomSunVectors(numberSimulations, this.latitude, this.longitude); + if (this.elevationRaster.length > 0) { + let shadingElevationAngles = elevation.getMaxElevationAngles( + this.elevationRaster, + this.elevationRasterMidpoint, + sunDirections.spherical.length / 2, + ); + } + //TODO: add shading of elevation here return rayTracingWebGL(midpoints, normals, meshArray, sunDirections, progressCallback); + } } diff --git a/src/rayTracingWebGL.ts b/src/rayTracingWebGL.ts index fdbbea7..56c7485 100644 --- a/src/rayTracingWebGL.ts +++ b/src/rayTracingWebGL.ts @@ -10,7 +10,10 @@ export function rayTracingWebGL( pointsArray: TypedArray, normals: TypedArray, trianglesArray: TypedArray, - sunDirections: Float32Array, + sunDirections: { + cartesian: Float32Array; + spherical: Float32Array; + }, progressCallback: (progress: number, total: number) => void, ): Float32Array | null { const N_TRIANGLES = trianglesArray.length / 9; @@ -179,11 +182,17 @@ export function rayTracingWebGL( var colorCodedArray = null; var isShadowedArray = null; - for (var i = 0; i < sunDirections.length; i += 3) { - progressCallback(i/3, sunDirections.length/3); + + for (var i = 0; i < sunDirections.cartesian.length; i += 3) { + progressCallback(i/3, sunDirections.cartesian.length/3); + // TODO: Iterate over sunDirection let sunDirectionUniformLocation = gl.getUniformLocation(program, 'u_sun_direction'); - gl.uniform3fv(sunDirectionUniformLocation, [sunDirections[i], sunDirections[i + 1], sunDirections[i + 2]]); + gl.uniform3fv(sunDirectionUniformLocation, [ + sunDirections.cartesian[i], + sunDirections.cartesian[i + 1], + sunDirections.cartesian[i + 2], + ]); drawArraysWithTransformFeedback(gl, tf, gl.POINTS, N_POINTS); diff --git a/src/sun.ts b/src/sun.ts index f6b41cf..71ed7cd 100644 --- a/src/sun.ts +++ b/src/sun.ts @@ -1,21 +1,39 @@ import { getPosition } from 'suncalc'; -export function getRandomSunVectors(Ndates: number, lat: number, lon: number): Float32Array { +/** + * Creates arrays of sun vectors. "cartesian" is a vector of length 3*Ndates where every three entries make up one vector. + * "spherical" is a vector of length 2*Ndates, where pairs of entries are altitude, azimuth. + * @param Ndates + * @param lat + * @param lon + * @returns + */ +export function getRandomSunVectors( + Ndates: number, + lat: number, + lon: number, +): { + cartesian: Float32Array; + spherical: Float32Array; +} { const sunVectors = new Float32Array(Ndates * 3); + const sunVectorsSpherical = new Float32Array(Ndates * 2); var i = 0; while (i < Ndates) { let date = getRandomDate(new Date(2023, 1, 1), new Date(2023, 12, 31)); - const pos = getPosition(date, lat, lon); - if (pos.altitude < 0.1 || pos.altitude == Number.NaN) { + const posSperical = getPosition(date, lat, lon); + if (posSperical.altitude < 0.1 || posSperical.altitude == Number.NaN) { continue; } - sunVectors[3 * i] = -Math.cos(pos.altitude) * Math.sin(pos.azimuth); - sunVectors[3 * i + 1] = -Math.cos(pos.altitude) * Math.cos(pos.azimuth); - sunVectors[3 * i + 2] = Math.sin(pos.altitude); + sunVectors[3 * i] = -Math.cos(posSperical.altitude) * Math.sin(posSperical.azimuth); + sunVectors[3 * i + 1] = -Math.cos(posSperical.altitude) * Math.cos(posSperical.azimuth); + sunVectors[3 * i + 2] = Math.sin(posSperical.altitude); + sunVectorsSpherical[2 * i] = posSperical.altitude; + sunVectorsSpherical[2 * i + 1] = posSperical.azimuth; i += 1; } - return sunVectors; + return { cartesian: sunVectors, spherical: sunVectorsSpherical }; } function getRandomDate(start: Date, end: Date): Date { diff --git a/tests/elevation.test.ts b/tests/elevation.test.ts new file mode 100644 index 0000000..6bb0715 --- /dev/null +++ b/tests/elevation.test.ts @@ -0,0 +1,61 @@ +import { describe, expect, test } from 'vitest'; +import * as elevation from '../src/elevation'; + +describe('calculateSphericalCoordinates', () => { + test('should calculate the correct spherical coordinates', () => { + const start = { x: 0, y: 0, z: 0 }; + const ends = [ + { x: 0, y: 0, z: 1 }, + { x: 1, y: 0, z: 1 }, + { x: 0, y: -1, z: 1 }, + ]; + const expectedResults = [ + { altitude: 0, azimuth: 0 }, + { altitude: Math.PI / 4, azimuth: 0 }, + { altitude: Math.PI / 4, azimuth: (3 / 2) * Math.PI }, + ]; + + ends.forEach((end, index) => { + const result = elevation.calculateSphericalCoordinates(start, end); + const expected = expectedResults[index]; + expect(result.azimuth).toBeCloseTo(expected.azimuth); + expect(result.altitude).toBeCloseTo(expected.altitude); + }); + }); +}); + +describe('fillMissingAltitudes', () => { + test('should fill negative infinity altitude with the nearest non-negative infinity altitude', () => { + const points: elevation.SphericalPoint[] = [ + { altitude: -Infinity, azimuth: 0 }, + { altitude: 10, azimuth: 90 }, + { altitude: -Infinity, azimuth: 180 }, + { altitude: -Infinity, azimuth: 230 }, + { altitude: -Infinity, azimuth: 240 }, + { altitude: 20, azimuth: 270 }, + ]; + + elevation.fillMissingAltitudes(points); + + expect(points[0].altitude).toBe(10); + expect(points[2].altitude).toBe(10); + expect(points[3].altitude).toBe(20); + expect(points[4].altitude).toBe(20); + }); +}); + +describe('getMaxElevationAngles', () => { + test('should correctly calculate the maximum elevation angles for given elevation points and observer', () => { + const elevations: elevation.Point[] = [ + { x: 1, y: 1, z: 2 }, + { x: 1, y: -1, z: 4 }, + { x: -1, y: -1, z: 6 }, + { x: -1, y: 1, z: 8 }, + ]; + const observer: elevation.Point = { x: 0, y: 0, z: 0 }; + const numDirections = 20; + const result: elevation.SphericalPoint[] = elevation.getMaxElevationAngles(elevations, observer, numDirections); + console.log(result); + expect(result).to.be.an('array').that.has.lengthOf(numDirections); + }); +}); diff --git a/tests/sun.test.ts b/tests/sun.test.ts index bcfa57b..03bf7d7 100644 --- a/tests/sun.test.ts +++ b/tests/sun.test.ts @@ -4,19 +4,24 @@ import * as sun from '../src/sun'; describe('Test functionalities from sun.ts: ', () => { const N = 50; let vectors = sun.getRandomSunVectors(N, 0, 0); - test('Get Correct number of positions.', () => { - expect(vectors.length).toStrictEqual(3 * N); + test('Get Correct number of positions for cartesian coordiantes.', () => { + expect(vectors.cartesian.length).toStrictEqual(3 * N); + }); + test('Get Correct number of positions for spherical coordiantes.', () => { + expect(vectors.spherical.length).toStrictEqual(2 * N); }); test('Get normalized sun vectors.', () => { for (let i = 0; i < N / 3; i++) { - let length = vectors[3 * i] ** 2 + vectors[3 * i + 1] ** 2 + vectors[3 * i + 2] ** 2; + let length = vectors.cartesian[3 * i] ** 2 + vectors.cartesian[3 * i + 1] ** 2 + vectors.cartesian[3 * i + 2] ** 2; expect(length).to.closeTo(1, 0.001); } }); test('Sun is always above the horizon.', () => { for (let i = 0; i < N / 3; i++) { - let z = vectors[3 * i + 2]; + let z = vectors.cartesian[3 * i + 2]; + let altitude = vectors.spherical[2 * i]; expect(z).toBeGreaterThan(0); + expect(altitude).toBeGreaterThan(0); } }); });