Add shading from rasterized elevation data #5

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

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<BufferGeometry>;
   shadingGeometries: Array<BufferGeometry>;
+  elevationRaster: Array<elevation.Point>;
+  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);
+
   }
 }

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

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 {

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

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