Merge branch 'main' into djairo-vtk-camera

vtk/src/CMakeLists.txt

@@ -2,6 +2,9 @@ cmake_minimum_required(VERSION 3.12 FATAL_ERROR)
 
 project(VtkBase)
 
+set(CMAKE_CXX_STANDARD 23)
+set(CMAKE_CXX_STANDARD_REQUIRED ON)
+
 set(CMAKE_EXPORT_COMPILE_COMMANDS ON)
 
 find_package(VTK COMPONENTS 
@@ -51,6 +54,20 @@ add_executable(VtkBase MACOSX_BUNDLE main.cpp
   layers/LGlyphLayer.h
   Program.cpp
   Program.h
+        CartographicTransformation.cpp
+        advection/AdvectionKernel.h
+        advection/EulerAdvectionKernel.cpp
+        advection/EulerAdvectionKernel.h
+        advection/interpolate.cpp
+        advection/interpolate.h
+        advection/readdata.cpp
+        advection/readdata.h
+        advection/RK4AdvectionKernel.cpp
+        advection/RK4AdvectionKernel.h
+        advection/UVGrid.cpp
+        advection/UVGrid.h
+        advection/Vel.cpp
+        advection/Vel.h
 )
 
 execute_process(

vtk/src/CartographicTransformation.cpp

@@ -15,11 +15,11 @@ vtkSmartPointer<vtkCamera> createNormalisedCamera() {
     return camera;
 }
 
-vtkSmartPointer<vtkMatrix4x4> getCartographicTransformMatrix() {
-    const double XMin = -15.875;
-    const double XMax = 12.875;
-    const double YMin = 46.125;
-    const double YMax = 62.625;
+vtkSmartPointer<vtkMatrix4x4> getCartographicTransformMatrix(const std::shared_ptr<UVGrid> uvGrid) {
+    const double XMin = uvGrid->lons.front();
+    const double XMax = uvGrid->lons.back();
+    const double YMin = uvGrid->lats.front();
+    const double YMax = uvGrid->lats.back();
 
     double eyeTransform[] = {
             2/(XMax-XMin), 0, 0, -(XMax+XMin)/(XMax-XMin),
@@ -34,10 +34,10 @@ vtkSmartPointer<vtkMatrix4x4> getCartographicTransformMatrix() {
 }
 
 // Assumes Normalised camera is used
-vtkSmartPointer<vtkTransformFilter> createCartographicTransformFilter() {
+vtkSmartPointer<vtkTransformFilter> createCartographicTransformFilter(const std::shared_ptr<UVGrid> uvGrid) {
     vtkNew<vtkTransform> transform;
 
-    transform->SetMatrix(getCartographicTransformMatrix());
+    transform->SetMatrix(getCartographicTransformMatrix(uvGrid));
 
     vtkSmartPointer<vtkTransformFilter> transformFilter = vtkSmartPointer<vtkTransformFilter>::New();
     transformFilter->SetTransform(transform);

vtk/src/CartographicTransformation.h

@@ -1,5 +1,7 @@
+#include <memory>
 #include <vtkCamera.h>
 #include <vtkTransformFilter.h>
+#include "advection/UVGrid.h"
 
 #ifndef VTKBASE_NORMALISEDCARTOGRAPHICCAMERA_H
 #define VTKBASE_NORMALISEDCARTOGRAPHICCAMERA_H
@@ -16,15 +18,14 @@ vtkSmartPointer<vtkCamera> createNormalisedCamera();
 /**
  * Constructs a 4x4 projection matrix that maps homogenious (longitude, latitude, 0, 1) points
  * to the normalised space.
- * TODO: This will soon require UVGrid as a parameter after the advection code is merged properly.
  * TODO: This transformation has room for improvement see:
  * https://github.com/MakeNEnjoy/interactive-track-and-trace/issues/12
  * @return pointer to 4x4 matrix
  */
-vtkSmartPointer<vtkMatrix4x4> getCartographicTransformMatrix();
+vtkSmartPointer<vtkMatrix4x4> getCartographicTransformMatrix(const std::shared_ptr<UVGrid> uvGrid);
 
 /**
  * Convenience function that converts the 4x4 projection matrix into a vtkTransformFilter
  * @return pointer to transform filter
  */
-vtkSmartPointer<vtkTransformFilter> createCartographicTransformFilter();
+vtkSmartPointer<vtkTransformFilter> createCartographicTransformFilter(const std::shared_ptr<UVGrid> uvGrid);

vtk/src/Program.cpp

@@ -34,9 +34,10 @@ void Program::setWinProperties() {
   interact->SetInteractorStyle(style);
 }
 
-void Program::setupTimer() {
+void Program::setupTimer(int dt) {
   auto callback = vtkSmartPointer<TimerCallbackCommand>::New(this);
   callback->SetClientData(this);
+  callback->setDt(dt);
   this->interact->AddObserver(vtkCommand::TimerEvent, callback);
   this->interact->AddObserver(vtkCommand::KeyPressEvent, callback);
   this->interact->CreateRepeatingTimer(17); // 60 fps == 1000 / 60 == 16.7 ms per frame
@@ -49,14 +50,15 @@ void Program::setupCameraCallback() {
   this->interact->AddObserver(vtkCommand::KeyPressEvent, callback);
 }
 
-Program::Program() {
+
+Program::Program(int timerDT) {
   this->win = vtkSmartPointer<vtkRenderWindow>::New();
   this->interact = vtkSmartPointer<vtkRenderWindowInteractor>::New();
   this->cam = createNormalisedCamera();
 
   this->win->SetNumberOfLayers(0);
   setWinProperties();
-  setupTimer();
+  setupTimer(timerDT);
   setupCameraCallback();
 }
 

vtk/src/Program.h

@@ -35,7 +35,7 @@ private:
 
   /** This function sets up and connects a TimerCallbackCommand with the program.
     */ 
-  void setupTimer();
+  void setupTimer(int dt);
 
   /** This function adds all interactors of each layer to the interactor/window
     */ 
@@ -48,7 +48,7 @@ private:
 public:
   /** Constructor.
     */ 
-  Program();
+  Program(int timerDT);
 
   /** This function adds a new layer (and thus vtkRenderer) to the program.
     * The layer is expected to set its own position in the vtkRenderWindow layer system.

vtk/src/advection/AdvectionKernel.h

@@ -0,0 +1,31 @@
+#ifndef ADVECTION_ADVECTIONKERNEL_H
+#define ADVECTION_ADVECTIONKERNEL_H
+
+#include <tuple>
+
+#include "Vel.h"
+
+/*
+ * Implement this class for every integration method.
+ */
+class AdvectionKernel {
+public:
+    /**
+     * This function must take a time, latitude and longitude of a particle and must output
+     * a new latitude and longitude after being advected once for AdvectionKernel::DT time as defined above.
+     * @param time Time since the beginning of the data
+     * @param latitude Latitude of particle
+     * @param longitude Longitude of particle
+     * @return A pair of latitude and longitude of particle.
+     */
+    virtual std::pair<double, double> advect(int time, double latitude, double longitude, int dt) const = 0;
+
+    // Taken from Parcels https://github.com/OceanParcels/parcels/blob/daa4b062ed8ae0b2be3d87367d6b45599d6f95db/parcels/tools/converters.py#L155
+    const static double metreToDegrees(double metre) {
+        return metre / 1000. / 1.852 / 60.;
+    }
+
+    virtual ~AdvectionKernel() = default; // Apparently I need this, idk why
+};
+
+#endif //ADVECTION_ADVECTIONKERNEL_H

vtk/src/advection/EulerAdvectionKernel.cpp

@@ -0,0 +1,13 @@
+
+#include "EulerAdvectionKernel.h"
+#include "interpolate.h"
+
+using namespace std;
+
+EulerAdvectionKernel::EulerAdvectionKernel(std::shared_ptr<UVGrid> grid) : grid(grid) {}
+
+std::pair<double, double> EulerAdvectionKernel::advect(int time, double latitude, double longitude, int dt) const {
+  auto [u, v] = bilinearinterpolate(*grid, time, latitude, longitude);
+
+  return {latitude + metreToDegrees(v * dt), longitude + metreToDegrees(u * dt)};
+}

vtk/src/advection/EulerAdvectionKernel.h

@@ -0,0 +1,26 @@
+#ifndef ADVECTION_EULERADVECTIONKERNEL_H
+#define ADVECTION_EULERADVECTIONKERNEL_H
+
+#include "AdvectionKernel.h"
+#include "UVGrid.h"
+#include <memory>
+
+/**
+ * Implementation of AdvectionKernel.
+ * The basic equation is:
+ * new_latitude = latitude + v*DT
+ * new_longitude = longitude + u*DT
+ *
+ * Uses bilinear interpolation as implemented in interpolate.h
+ */
+class EulerAdvectionKernel: public AdvectionKernel {
+private:
+    std::shared_ptr<UVGrid> grid;
+public:
+    explicit EulerAdvectionKernel(std::shared_ptr<UVGrid> grid);
+    std::pair<double, double> advect(int time, double latitude, double longitude, int dt) const override;
+
+};
+
+
+#endif //ADVECTION_EULERADVECTIONKERNEL_H

vtk/src/advection/RK4AdvectionKernel.cpp

@@ -0,0 +1,35 @@
+#include "RK4AdvectionKernel.h"
+#include "interpolate.h"
+
+using namespace std;
+
+RK4AdvectionKernel::RK4AdvectionKernel(std::shared_ptr<UVGrid> grid): grid(grid) { }
+
+std::pair<double, double> RK4AdvectionKernel::advect(int time, double latitude, double longitude, int dt) const {
+    auto [u1, v1] = bilinearinterpolate(*grid, time, latitude, longitude);
+//    lon1, lat1 = (particle.lon + u1*.5*particle.dt, particle.lat + v1*.5*particle.dt);
+    double lon1 = longitude + metreToDegrees(u1 * 0.5*dt);
+    double lat1 = latitude + metreToDegrees(v1 * 0.5*dt);
+
+//    (u2, v2) = fieldset.UV[time + .5 * particle.dt, particle.depth, lat1, lon1, particle]
+    auto [u2, v2] = bilinearinterpolate(*grid, time + 0.5 * dt, lat1, lon1);
+
+//    lon2, lat2 = (particle.lon + u2*.5*particle.dt, particle.lat + v2*.5*particle.dt)
+    double lon2 = longitude + metreToDegrees(u2 * 0.5 * dt);
+    double lat2 = latitude + metreToDegrees(v2 * 0.5 * dt);
+
+//    (u3, v3) = fieldset.UV[time + .5 * particle.dt, particle.depth, lat2, lon2, particle]
+    auto [u3, v3] = bilinearinterpolate(*grid, time + 0.5 * dt, lat2, lon2);
+
+//    lon3, lat3 = (particle.lon + u3*particle.dt, particle.lat + v3*particle.dt)
+    double lon3 = longitude + metreToDegrees(u3 * dt);
+    double lat3 = latitude + metreToDegrees(v3 * dt);
+
+//    (u4, v4) = fieldset.UV[time + particle.dt, particle.depth, lat3, lon3, particle]
+    auto [u4, v4] = bilinearinterpolate(*grid, time + dt, lat3, lon3);
+
+    double lonFinal = longitude + metreToDegrees((u1 + 2 * u2 + 2 * u3 + u4) / 6.0 * dt);
+    double latFinal = latitude + metreToDegrees((v1 + 2 * v2 + 2 * v3 + v4) / 6.0 * dt);
+
+    return {latFinal, lonFinal};
+}

vtk/src/advection/RK4AdvectionKernel.h

@@ -0,0 +1,23 @@
+#ifndef ADVECTION_RK4ADVECTIONKERNEL_H
+#define ADVECTION_RK4ADVECTIONKERNEL_H
+
+#include "AdvectionKernel.h"
+#include "UVGrid.h"
+#include <memory>
+
+/**
+ * Implementation of Advection kernel using RK4 integration
+ * See https://en.wikipedia.org/wiki/Runge%E2%80%93Kutta_methods for more details.
+ * Uses bilinear interpolation as implemented in interpolate.h
+ */
+class RK4AdvectionKernel: public AdvectionKernel {
+private:
+    std::shared_ptr<UVGrid> grid;
+public:
+    explicit RK4AdvectionKernel(std::shared_ptr<UVGrid> grid);
+    std::pair<double, double> advect(int time, double latitude, double longitude, int dt) const override;
+
+};
+
+
+#endif //ADVECTION_RK4ADVECTIONKERNEL_H

vtk/src/advection/UVGrid.cpp

@@ -0,0 +1,66 @@
+#include <ranges>
+
+#include "UVGrid.h"
+#include "readdata.h"
+
+#define sizeError2 "The sizes of the hydrodynamic data files are different"
+#define sizeError "The sizes of the hydrodynamicU or -V files does not correspond with the sizes of the grid file"
+
+using namespace std;
+
+UVGrid::UVGrid() {
+  auto us = readHydrodynamicU();
+  auto vs = readHydrodynamicV();
+  if (us.size() != vs.size()) {
+    throw domain_error(sizeError2);
+  }
+
+  tie(times, lats, lons) = readGrid();
+
+  timeSize = times.size();
+  latSize = lats.size();
+  lonSize = lons.size();
+
+  size_t gridSize = timeSize * latSize * lonSize;
+  if (gridSize != us.size()) {
+    throw domain_error(sizeError);
+  }
+
+  uvData.reserve(gridSize);
+
+  for (auto vel: views::zip(us, vs)) {
+    uvData.push_back(Vel(vel));
+  }
+}
+
+const Vel &UVGrid::operator[](size_t timeIndex, size_t latIndex, size_t lonIndex) const {
+  if (timeIndex < 0 or timeIndex >= timeSize
+      or latIndex < 0 or latIndex >= latSize
+      or lonIndex < 0 or lonIndex >= lonSize) {
+    throw std::out_of_range("Index out of bounds");
+  }
+  size_t index = timeIndex * (latSize * lonSize) + latIndex * lonSize + lonIndex;
+  return uvData[index];
+}
+
+double UVGrid::lonStep() const {
+  return lons[1] - lons[0];
+}
+
+double UVGrid::latStep() const {
+  return lats[1] - lats[0];
+}
+
+int UVGrid::timeStep() const {
+  return times[1] - times[0];
+}
+
+void UVGrid::streamSlice(ostream &os, size_t t) {
+  for (int x = 0; x < latSize; x++) {
+    for (int y = 0; y < lonSize; y++) {
+      auto vel = (*this)[t, x, y];
+      os << vel << " ";
+    }
+    os << endl;
+  }
+}

vtk/src/advection/UVGrid.h

@@ -0,0 +1,65 @@
+#ifndef ADVECTION_UVGRID_H
+#define ADVECTION_UVGRID_H
+
+#include <vector>
+#include "Vel.h"
+
+class UVGrid {
+private:
+    /**
+     * 1D data vector of all the us and vs
+     */
+    std::vector<Vel> uvData;
+public:
+    UVGrid();
+
+    /**
+     * The matrix has shape (timeSize, latSize, lonSize)
+     */
+    size_t timeSize;
+    size_t latSize;
+    size_t lonSize;
+
+    /**
+     * Assuming grid is a regular grid, gives the longitudinal spacing of grid.
+     * @return longitudinal spacing
+     */
+    double lonStep() const;
+
+    /**
+     * Assuming grid is a regular grid, gives the latitudinal spacing of grid.
+     * @return latitudinal spacing
+     */
+    double latStep() const;
+
+    /**
+     * Assuming grid is a regular grid, gives the time spacing of grid.
+     * @return time spacing
+     */
+    int timeStep() const;
+
+    /**
+     * times, lats, lons are vector of length timeSize, latSize, lonSize respectively.
+     * The maintain the following invariant:
+     * grid[timeIndex,latIndex,lonIndex] gives the u,v at the point with latitude at lats[latIndex],
+     * with longitude at lons[lonIndex], and with time at times[timeIndex].
+     */
+    std::vector<int> times;
+    std::vector<double> lats;
+    std::vector<double> lons;
+
+    /**
+     * The 3D index into the data. The array is sized by [8761][67][116]
+     * @return Velocity at that index
+     */
+    const Vel& operator[](size_t timeIndex, size_t latIndex, size_t lonIndex) const;
+
+    /**
+     * Streams a slice at timeIndex t of the matrix to the outstream given by os
+     * @param os outstream
+     * @param t index with which to slice matrix
+     */
+    void streamSlice(std::ostream &os, size_t t);
+};
+
+#endif //ADVECTION_UVGRID_H

vtk/src/advection/Vel.cpp

@@ -0,0 +1,40 @@
+#include "Vel.h"
+#include <stdexcept>
+#include <iomanip>
+
+using namespace std;
+
+Vel::Vel(double u, double v) : u(u), v(v) {}
+
+Vel::Vel(const std::pair<double, double>& p) : u(p.first), v(p.second) {}
+
+Vel& Vel::operator=(const std::pair<double, double>& p) {
+    u = p.first;
+    v = p.second;
+    return *this;
+}
+
+Vel Vel::operator+(const Vel& other) const {
+    return Vel(u + other.u, v + other.v);
+}
+
+Vel& Vel::operator+=(const Vel& other) {
+    u += other.u;
+    v += other.v;
+    return *this;
+}
+
+template<typename Scalar>
+Vel Vel::operator/(Scalar scalar) const {
+    if (scalar == 0) throw std::runtime_error("Division by zero");
+    return Vel(u / scalar, v / scalar);
+}
+
+std::ostream& operator<<(ostream& os, const Vel& vel) {
+    os << "(";
+    os << fixed << setprecision(2) << setw(5) << vel.u;
+    os << ", ";
+    os << fixed << setprecision(2) << setw(5) << vel.v;
+    os << ")";
+    return os;
+}

vtk/src/advection/Vel.h

@@ -0,0 +1,44 @@
+#ifndef ADVECTION_VEL_H
+#define ADVECTION_VEL_H
+
+#include <utility>
+#include <stdexcept>
+#include <iostream>
+#include <format>
+
+class Vel {
+public:
+    double u; // Eastward Current Velocity in the Water Column
+    double v; // Northward Current Velocity in the Water Column
+
+    Vel(double u, double v);
+    Vel(const std::pair<double, double>& p);  // Conversion constructor
+    Vel& operator=(const std::pair<double, double>& p);
+
+    // Operator + to add two Vel objects
+    Vel operator+(const Vel& other) const;
+
+    // Operator += to add another Vel object to this object
+    Vel& operator+=(const Vel& other);
+
+    // Operator * to multiply Vel by a scalar, defined as a member template
+    template<typename Scalar>
+    Vel operator*(Scalar scalar) const {
+        return Vel(u * scalar, v * scalar);
+    }
+
+    // Operator / to divide Vel by a scalar, defined as a member template
+    template<typename Scalar>
+    Vel operator/(Scalar scalar) const;
+
+    // Friend declaration for the stream insertion operator
+    friend std::ostream& operator<<(std::ostream& os, const Vel& vel);
+};
+
+// Non-member function for scalar multiplication on the left
+template<typename Scalar>
+Vel operator*(Scalar scalar, const Vel& p) {
+    return Vel(p.u * scalar, p.v * scalar);
+}
+
+#endif //ADVECTION_VEL_H

vtk/src/advection/interpolate.cpp

@@ -0,0 +1,47 @@
+#include "interpolate.h"
+
+using namespace std;
+
+Vel bilinearinterpolate(const UVGrid &uvGrid, int time, double lat, double lon) {
+    double latStep = uvGrid.latStep();
+    double lonStep = uvGrid.lonStep();
+    int timeStep = uvGrid.timeStep();
+
+    int latIndex = (lat - uvGrid.lats[0]) / latStep;
+    int lonIndex = (lon - uvGrid.lons[0]) / lonStep;
+    int timeIndex = (time - uvGrid.times[0]) / timeStep;
+
+    double timeRatio = (static_cast<double>(time) - uvGrid.times[timeIndex]) / timeStep;
+    double latRatio = (lat - uvGrid.lats[latIndex]) / latStep;
+    double lonRatio = (lon - uvGrid.lons[lonIndex]) / lonStep;
+
+    Vel point = {0, 0};
+    for (int timeOffset = 0; timeOffset <= 1; timeOffset++) {
+        for (int latOffset = 0; latOffset <= 1; latOffset++) {
+            for (int lonOffset = 0; lonOffset <= 1; lonOffset++) {
+                auto vertex = uvGrid[
+                        timeIndex + 1 < uvGrid.timeSize ? timeIndex + timeOffset : timeIndex,
+                                latIndex + 1 < uvGrid.latSize ? latIndex + latOffset : latIndex,
+                                lonIndex + 1 < uvGrid.lonSize ? lonIndex + lonOffset : lonIndex
+                ];
+
+                double timeRation = (1 - timeOffset) * (1 - timeRatio) + timeOffset * timeRatio;
+                double latRation = (1 - latOffset) * (1 - latRatio) + latOffset * latRatio;
+                double lonRation = (1 - lonOffset) * (1 - lonRatio) + lonOffset * lonRatio;
+                point += timeRation * latRation * lonRation * vertex;
+            }
+        }
+    }
+
+    return point;
+}
+
+vector<Vel> bilinearinterpolation(const UVGrid &uvGrid, vector<tuple<int, double, double>> points) {
+    vector<Vel> result;
+    result.reserve(points.size());
+    for (auto [time, lat, lon]: points) {
+        result.push_back(bilinearinterpolate(uvGrid, time, lat, lon));
+    }
+
+    return result;
+}

vtk/src/advection/interpolate.h

@@ -0,0 +1,28 @@
+#ifndef ADVECTION_INTERPOLATE_H
+#define ADVECTION_INTERPOLATE_H
+
+#include <vector>
+
+#include "UVGrid.h"
+
+/**
+ * Bilinearly interpolate the point (time, lat, lon) to produce the interpolated velocity.
+ * Since it is in 3D, this means that it interpolates against 8 points (excluding edges).
+ * As described in https://numerical.recipes/book.html Chapter 3.6
+ * @param uvGrid velocity grid
+ * @param time time of point
+ * @param lat latitude of point
+ * @param lon longitude of point
+ * @return interpolated velocity
+ */
+Vel bilinearinterpolate(const UVGrid &uvGrid, int time, double lat, double lon);
+
+/**
+ * Helper function for bilnearly interpolating a vector of points
+ * @param uvGrid velocity grid
+ * @param points vector of points
+ * @return interpolated velocities
+ */
+std::vector<Vel> bilinearinterpolation(const UVGrid &uvGrid, std::vector<std::tuple<int, double, double>> points);
+
+#endif //ADVECTION_INTERPOLATE_H

vtk/src/advection/readdata.cpp

@@ -0,0 +1,50 @@
+#include <stdexcept>
+
+#include <netcdf>
+
+#include "readdata.h"
+
+using namespace std;
+using namespace netCDF;
+
+template <typename T>
+vector<T> getVarVector(const NcVar &var) {
+    int length = 1;
+    for (NcDim dim : var.getDims()) {
+        length *= dim.getSize();
+    }
+
+    vector<T> vec(length);
+
+    var.getVar(vec.data());
+
+    return vec;
+}
+
+vector<double> readHydrodynamicU() {
+    // Vs and Us flipped cause the files are named incorrectly
+    netCDF::NcFile data("../../../../data/hydrodynamic_V.h5", netCDF::NcFile::read);
+
+    multimap< string, NcVar > vars = data.getVars();
+
+    return getVarVector<double>(vars.find("vo")->second);
+}
+
+vector<double> readHydrodynamicV() {
+    // Vs and Us flipped cause the files are named incorrectly
+    netCDF::NcFile data("../../../../data/hydrodynamic_U.h5", netCDF::NcFile::read);
+
+    multimap< string, NcVar > vars = data.getVars();
+
+    return getVarVector<double>(vars.find("uo")->second);
+}
+
+tuple<vector<int>, vector<double>, vector<double>> readGrid() {
+    netCDF::NcFile data("../../../../data/grid.h5", netCDF::NcFile::read);
+    multimap< string, NcVar > vars = data.getVars();
+    vector<int> time = getVarVector<int>(vars.find("times")->second);
+    vector<double> longitude = getVarVector<double>(vars.find("longitude")->second);
+    vector<double> latitude = getVarVector<double>(vars.find("latitude")->second);
+
+    return {time, latitude, longitude};
+}

vtk/src/advection/readdata.h

@@ -0,0 +1,22 @@
+#ifndef ADVECTION_READDATA_H
+#define ADVECTION_READDATA_H
+
+/**
+ * reads the file hydrodynamic_U.h5
+ * @return the data vector of us
+ */
+std::vector<double> readHydrodynamicU();
+
+/**
+ * reads the file hydrodynamic_V.h5
+ * @return the data vector of vs
+ */
+std::vector<double> readHydrodynamicV();
+
+/**
+ * Reads the file grid.h5
+ * @return a tuple of (times, latitude, longitude)
+ */
+std::tuple<std::vector<int>, std::vector<double>, std::vector<double>> readGrid();
+
+#endif //ADVECTION_READDATA_H

vtk/src/commands/SpawnPointCallback.cpp

@@ -9,11 +9,11 @@
 
 #include "../CartographicTransformation.h"
 
-void convertDisplayToWorld(vtkRenderer* renderer, int x, int y, double *worldPos) {
-    double displayPos[3] = {static_cast<double>(x), static_cast<double>(y), 0.0};
-    renderer->SetDisplayPoint(displayPos);
-    renderer->DisplayToWorld();
-    renderer->GetWorldPoint(worldPos);
+void convertDisplayToWorld(vtkRenderer *renderer, int x, int y, double *worldPos) {
+  double displayPos[3] = {static_cast<double>(x), static_cast<double>(y), 0.0};
+  renderer->SetDisplayPoint(displayPos);
+  renderer->DisplayToWorld();
+  renderer->GetWorldPoint(worldPos);
 }
 
 void SpawnPointCallback::Execute(vtkObject *caller, unsigned long evId, void *callData) {
@@ -54,23 +54,29 @@ void SpawnPointCallback::Execute(vtkObject *caller, unsigned long evId, void *ca
 }
 
 
-SpawnPointCallback::SpawnPointCallback() : data(nullptr), points(nullptr), inverseCartographicProjection(nullptr) {
-    inverseCartographicProjection = getCartographicTransformMatrix();
-    inverseCartographicProjection->Invert();
-}
+SpawnPointCallback::SpawnPointCallback() : data(nullptr),
+                                           points(nullptr),
+                                           inverseCartographicProjection(nullptr),
+                                           uvGrid(nullptr) { }
 
 SpawnPointCallback *SpawnPointCallback::New() {
     return new SpawnPointCallback;
 }
 
 void SpawnPointCallback::setData(const vtkSmartPointer<vtkPolyData> &data) {
-    this->data = data;
+  this->data = data;
 }
 
 void SpawnPointCallback::setPoints(const vtkSmartPointer<vtkPoints> &points) {
-    this->points = points;
+  this->points = points;
 }
 
 void SpawnPointCallback::setRen(const vtkSmartPointer<vtkRenderer> &ren) {
-    this->ren = ren;
+  this->ren = ren;
+}
+
+void SpawnPointCallback::setUVGrid(const std::shared_ptr<UVGrid> &uvGrid) {
+  this->uvGrid = uvGrid;
+  inverseCartographicProjection = getCartographicTransformMatrix(uvGrid);
+  inverseCartographicProjection->Invert();
 }

vtk/src/commands/SpawnPointCallback.h

@@ -2,31 +2,39 @@
 #define VTKBASE_SPAWNPOINTCALLBACK_H
 
 
+#include <memory>
 #include <vtkCallbackCommand.h>
 #include <vtkRenderWindowInteractor.h>
 #include <vtkPoints.h>
 #include <vtkPolyData.h>
 #include <vtkMatrix4x4.h>
+#include "../advection/UVGrid.h"
 
 class SpawnPointCallback : public vtkCallbackCommand {
 
 public:
-    static SpawnPointCallback *New();
-    SpawnPointCallback();
+  static SpawnPointCallback *New();
 
-    void setPoints(const vtkSmartPointer<vtkPoints> &points);
+  SpawnPointCallback();
 
-    void setData(const vtkSmartPointer<vtkPolyData> &data);
+  void setPoints(const vtkSmartPointer<vtkPoints> &points);
+
+  void setData(const vtkSmartPointer<vtkPolyData> &data);
+
+  void setRen(const vtkSmartPointer<vtkRenderer> &ren);
+
+  void setUVGrid(const std::shared_ptr<UVGrid> &uvGrid);
 
-    void setRen(const vtkSmartPointer<vtkRenderer> &ren);
 private:
-    vtkSmartPointer<vtkPolyData> data;
-    vtkSmartPointer<vtkPoints> points;
-    vtkSmartPointer<vtkRenderer> ren;
-    vtkSmartPointer<vtkMatrix4x4> inverseCartographicProjection;
+  vtkSmartPointer<vtkPolyData> data;
+  vtkSmartPointer<vtkPoints> points;
+  vtkSmartPointer<vtkRenderer> ren;
+  std::shared_ptr<UVGrid> uvGrid;
+  vtkSmartPointer<vtkMatrix4x4> inverseCartographicProjection;
 
-    void Execute(vtkObject *caller, unsigned long evId, void *callData) override;
-    bool dragging = false;
+  void Execute(vtkObject *caller, unsigned long evId, void *callData) override;
+
+  bool dragging = false;
 };
 
 

vtk/src/commands/TimerCallbackCommand.cpp

@@ -37,3 +37,7 @@ void TimerCallbackCommand::setProgram(Program *program) {
 void TimerCallbackCommand::setPaused(const bool val) {
   this->paused = val;
 }
+
+void TimerCallbackCommand::setDt(int dt) {
+  this->dt = dt;
+}

vtk/src/commands/TimerCallbackCommand.h

@@ -13,6 +13,8 @@ public:
   void setProgram(Program *program);
   void setPaused(const bool val);
 
+  void setDt(int dt);
+
 private:
   int time;
   int dt;

vtk/src/layers/EGlyphLayer.cpp

@@ -11,43 +11,20 @@
 #include <vtkProperty.h>
 #include <vtkProperty2D.h>
 #include <vtkVertexGlyphFilter.h>
-#include <netcdf>
 #include <vtkArrowSource.h>
 #include "../CartographicTransformation.h"
+#include "../advection/readdata.h"
+#include "../advection/interpolate.h"
 
-using namespace netCDF;
 using namespace std;
 
-template <typename T>
-vector<T> getVarVector(const NcVar &var) {
-    int length = 1;
-    for (NcDim dim : var.getDims()) {
-        length *= dim.getSize();
-    }
-
-    vector<T> vec(length);
-
-    var.getVar(vec.data());
-
-    return vec;
-}
-
-tuple<vector<int>, vector<double>, vector<double>> readGrid() {
-    netCDF::NcFile data("../../../../data/grid.h5", netCDF::NcFile::read);
-    multimap< string, NcVar > vars = data.getVars();
-    vector<int> time = getVarVector<int>(vars.find("times")->second);
-    vector<double> longitude = getVarVector<double>(vars.find("longitude")->second);
-    vector<double> latitude = getVarVector<double>(vars.find("latitude")->second);
-
-    return {time, latitude, longitude};
-}
-
-
-EGlyphLayer::EGlyphLayer() {
+EGlyphLayer::EGlyphLayer(std::shared_ptr<UVGrid> uvGrid) {
   this->ren = vtkSmartPointer<vtkRenderer>::New();
   this->ren->SetLayer(1);
   this->ren->InteractiveOff();
 
+  this->uvGrid = uvGrid;
+
   this->data = vtkSmartPointer<vtkPolyData>::New();
   this->direction = vtkSmartPointer<vtkDoubleArray>::New();
   this->direction->SetName("direction");
@@ -57,28 +34,30 @@ EGlyphLayer::EGlyphLayer() {
 
 void EGlyphLayer::readCoordinates() {
   vtkNew<vtkPoints> points;
-  auto [times, lats, lons] = readGrid(); // FIXME: import Robin's readData function and use it
-  this->numLats = lats.size();
-  this->numLons = lons.size();
+  this->numLats = uvGrid->lats.size();
+  this->numLons = uvGrid->lons.size();
 
   this->direction->SetNumberOfComponents(3);
-  this->direction->SetNumberOfTuples(numLats*numLons); 
-  points->Allocate(numLats*numLons);
+  this->direction->SetNumberOfTuples(numLats * numLons);
+  points->Allocate(numLats * numLons);
 
   int i = 0;
-  for (double lat : lats) {
-    for (double lon : lons) {
-        // cout << "lon: " << lon << " lat: " << lat << endl;
-      direction->SetTuple3(i, 0.45, 0.90, 0); //FIXME: read this info from file
+  int latIndex = 0;
+  for (double lat: uvGrid->lats) {
+    int lonIndex = 0;
+    for (double lon: uvGrid->lons) {
+      auto [u, v] = (*uvGrid)[0, latIndex, lonIndex];
+      direction->SetTuple3(i, 5*u, 5*v, 0);
       points->InsertPoint(i++, lon, lat, 0);
-      // see also https://vtk.org/doc/nightly/html/classvtkPolyDataMapper2D.html
+      lonIndex++;
     }
+    latIndex++;
   }
   this->data->SetPoints(points);
   this->data->GetPointData()->AddArray(this->direction);
   this->data->GetPointData()->SetActiveVectors("direction");
 
-  vtkSmartPointer<vtkTransformFilter> transformFilter = createCartographicTransformFilter();
+  vtkSmartPointer<vtkTransformFilter> transformFilter = createCartographicTransformFilter(uvGrid);
   transformFilter->SetInputData(data);
 
   vtkNew<vtkGlyphSource2D> arrowSource;
@@ -91,31 +70,33 @@ void EGlyphLayer::readCoordinates() {
   glyph2D->SetInputConnection(transformFilter->GetOutputPort());
   glyph2D->OrientOn();
   glyph2D->ClampingOn();
-  glyph2D->SetScaleModeToScaleByVector(); 
-  glyph2D->SetVectorModeToUseVector(); 
+  glyph2D->SetScaleModeToScaleByVector();
+  glyph2D->SetVectorModeToUseVector();
   glyph2D->Update();
 
-//  vtkNew<vtkCoordinate> coordinate;
-//  coordinate->SetCoordinateSystemToWorld();
-
   vtkNew<vtkPolyDataMapper>(mapper);
-  // mapper->SetTransformCoordinate(coordinate);
   mapper->SetInputConnection(glyph2D->GetOutputPort());
   mapper->Update();
 
   vtkNew<vtkActor> actor;
   actor->SetMapper(mapper);
 
-  actor->GetProperty()->SetColor(0,0,0);
+  actor->GetProperty()->SetColor(0, 0, 0);
   actor->GetProperty()->SetOpacity(0.2);
 
-  this->ren->AddActor(actor)  ;
+  this->ren->AddActor(actor);
 }
 
 
 void EGlyphLayer::updateData(int t) {
-  for (int i=0; i < numLats*numLons; i++) {
-    this->direction->SetTuple3(i, std::cos(t), std::sin(t), 0); // FIXME: fetch data from file.
+  int i = 0;
+  for (int lat = 0; lat < uvGrid->latSize; lat++) {
+    for (int lon = 0; lon < uvGrid->lonSize; lon++) {
+      auto [u, v] = (*uvGrid)[t/3600, lat, lon];
+      // TODO: 5*v scaling stuff should really be a filter transform
+      this->direction->SetTuple3(i, 5*u, 5*v, 0);
+      i++;
+    }
   }
   this->direction->Modified();
 }

vtk/src/layers/EGlyphLayer.h

@@ -2,8 +2,11 @@
 #define EGLYPHLAYER_H
 
 #include "Layer.h"
+#include <memory>
 #include <vtkPolyData.h>
 
+#include "../advection/UVGrid.h"
+
 /** Implements the Layer class for the case of a Eulerian visualization.
   * Specifically, this class models the eulerian flow-fields of the simulation using the 'glyph' mark and 'direction' and 'form' channels to denote direction and strength of velocities.
   */
@@ -11,6 +14,7 @@ class EGlyphLayer : public Layer {
 private:
   vtkSmartPointer<vtkPolyData> data;
   vtkSmartPointer<vtkDoubleArray> direction;
+  std::shared_ptr<UVGrid> uvGrid;
   int numLats;
   int numLons;
 
@@ -22,7 +26,7 @@ private:
 public:
   /** Constructor.
     */
-  EGlyphLayer();
+  EGlyphLayer(std::shared_ptr<UVGrid> uvGrid);
 
   /** updates the glyphs to reflect the given timestamp in the dataset.
     * @param t : the time at which to fetch the data.

vtk/src/layers/LGlyphLayer.cpp

@@ -18,13 +18,13 @@
 
 #include "../CartographicTransformation.h"
 
-
 vtkSmartPointer<SpawnPointCallback> LGlyphLayer::createSpawnPointCallback() {
-    auto newPointCallBack = vtkSmartPointer<SpawnPointCallback>::New();
-    newPointCallBack->setData(data);
-    newPointCallBack->setPoints(points);
-    newPointCallBack->setRen(ren);
-    return newPointCallBack;
+  auto newPointCallBack = vtkSmartPointer<SpawnPointCallback>::New();
+  newPointCallBack->setData(data);
+  newPointCallBack->setPoints(points);
+  newPointCallBack->setRen(ren);
+  newPointCallBack->setUVGrid(uvGrid);
+  return newPointCallBack;
 }
 
 // Further notes; current thinking is to allow tracking a particle's age by using a scalar array in the VtkPolyData. This would be incremented for every tick/updateData function call.
@@ -32,37 +32,39 @@ vtkSmartPointer<SpawnPointCallback> LGlyphLayer::createSpawnPointCallback() {
 //
 // TODO: modelling all this in vtkClasses is workable, but ideally i would want to work with a native C++ class. See if this is doable and feasible.
 
-LGlyphLayer::LGlyphLayer() {
-    this->ren = vtkSmartPointer<vtkRenderer>::New();
-    this->ren->SetLayer(2);
+LGlyphLayer::LGlyphLayer(std::shared_ptr<UVGrid> uvGrid, std::unique_ptr<AdvectionKernel> advectionKernel) {
+  this->ren = vtkSmartPointer<vtkRenderer>::New();
+  this->ren->SetLayer(2);
 
-    this->points = vtkSmartPointer<vtkPoints>::New();
-    this->data = vtkSmartPointer<vtkPolyData>::New();
-    this->data->SetPoints(this->points);
+  this->points = vtkSmartPointer<vtkPoints>::New();
+  this->data = vtkSmartPointer<vtkPolyData>::New();
+  this->data->SetPoints(this->points);
 
-    vtkSmartPointer<vtkTransformFilter> transformFilter = createCartographicTransformFilter();
-    transformFilter->SetInputData(data);
+  advector = std::move(advectionKernel);
+  this->uvGrid = uvGrid;
 
-    vtkNew<vtkGlyphSource2D> circleSource;
-    circleSource->SetGlyphTypeToCircle();
-    circleSource->SetScale(0.01);
-    circleSource->Update();
+  vtkSmartPointer<vtkTransformFilter> transformFilter = createCartographicTransformFilter(uvGrid);
+  transformFilter->SetInputData(data);
 
-    vtkNew<vtkGlyph2D> glyph2D;
-    glyph2D->SetSourceConnection(circleSource->GetOutputPort());
-    glyph2D->SetInputConnection(transformFilter->GetOutputPort());
-    glyph2D->SetColorModeToColorByScalar();
-    glyph2D->Update();
+  vtkNew<vtkGlyphSource2D> circleSource;
+  circleSource->SetGlyphTypeToCircle();
+  circleSource->SetScale(0.05);
+  circleSource->Update();
 
-    vtkNew<vtkPolyDataMapper> mapper;
-    mapper->SetInputConnection(glyph2D->GetOutputPort());
-    mapper->Update();
+  vtkNew<vtkGlyph2D> glyph2D;
+  glyph2D->SetSourceConnection(circleSource->GetOutputPort());
+  glyph2D->SetInputConnection(transformFilter->GetOutputPort());
+  glyph2D->SetColorModeToColorByScalar();
+  glyph2D->Update();
 
-    vtkNew<vtkActor> actor;
-    actor->SetMapper(mapper);
-  actor->GetProperty()->SetOpacity(0.8);
+  vtkNew<vtkPolyDataMapper> mapper;
+  mapper->SetInputConnection(glyph2D->GetOutputPort());
+  mapper->Update();
+  
+  vtkNew<vtkActor> actor;
+  actor->SetMapper(mapper);
 
-    this->ren->AddActor(actor);
+  this->ren->AddActor(actor);
 }
 
 // creates a few points so we can test the updateData function
@@ -82,27 +84,23 @@ void LGlyphLayer::spoofPoints() {
   this->points->Modified();
 }
 
-// returns new coords for a point; used to test the updateData function
-std::pair<double, double> advect(int time, double lat, double lon) {
-    return {lat + 0.01, lon + 0.01};
-}
-
 void LGlyphLayer::updateData(int t) {
-    double point[3];
-    for (vtkIdType n = 0; n < this->points->GetNumberOfPoints(); n++) {
-        this->points->GetPoint(n, point);
-        auto [xNew, yNew] = advect(n, point[0], point[1]);
-        this->points->SetPoint(n, xNew, yNew, 0);
+  const int SUPERSAMPLINGRATE = 4;
+  double point[3];
+  for (vtkIdType n = 0; n < this->points->GetNumberOfPoints(); n++) {
+    this->points->GetPoint(n, point);
+    for (int i = 0; i < SUPERSAMPLINGRATE; i++) {
+      std::tie(point[1], point[0]) = advector->advect(t, point[1], point[0], (t-lastT)/SUPERSAMPLINGRATE);
     }
-    this->points->Modified();
+    this->points->SetPoint(n, point[0], point[1], 0);
+  }
+  lastT = t;
+  this->points->Modified();
 }
 
 void LGlyphLayer::addObservers(vtkSmartPointer<vtkRenderWindowInteractor> interactor) {
-    auto newPointCallBack = vtkSmartPointer<SpawnPointCallback>::New();
-    newPointCallBack->setData(data);
-    newPointCallBack->setPoints(points);
-    newPointCallBack->setRen(ren);
-    interactor->AddObserver(vtkCommand::LeftButtonPressEvent, newPointCallBack);
-    interactor->AddObserver(vtkCommand::LeftButtonReleaseEvent, newPointCallBack);
-    interactor->AddObserver(vtkCommand::MouseMoveEvent, newPointCallBack);
+  auto newPointCallBack = createSpawnPointCallback();
+  interactor->AddObserver(vtkCommand::LeftButtonPressEvent, newPointCallBack);
+  interactor->AddObserver(vtkCommand::LeftButtonReleaseEvent, newPointCallBack);
+  interactor->AddObserver(vtkCommand::MouseMoveEvent, newPointCallBack);
 }

vtk/src/layers/LGlyphLayer.h

@@ -2,6 +2,7 @@
 #define LGLYPHLAYER_H
 
 #include "Layer.h"
+#include "../advection/AdvectionKernel.h"
 #include "../commands/SpawnPointCallback.h"
 #include <vtkPolyData.h>
 #include <vtkInteractorStyle.h>
@@ -13,13 +14,14 @@ class LGlyphLayer : public Layer {
 private:
   vtkSmartPointer<vtkPoints> points;
   vtkSmartPointer<vtkPolyData> data;
-
-
+  std::unique_ptr<AdvectionKernel> advector;
+  std::shared_ptr<UVGrid> uvGrid;
+  int lastT = 1000;
 
 public:
   /** Constructor.
     */
-  LGlyphLayer();
+  LGlyphLayer(std::shared_ptr<UVGrid> uvGrid, std::unique_ptr<AdvectionKernel> advectionKernel);
 
   /** This function spoofs a few points in the dataset. Mostly used for testing.
     */

vtk/src/main.cpp

@@ -7,25 +7,30 @@
 #include <vtkProperty2D.h>
 #include <vtkRenderer.h>
 #include <vtkVertexGlyphFilter.h>
+#include <memory>
 
 #include "layers/BackgroundImage.h"
 #include "layers/EGlyphLayer.h"
 #include "layers/LGlyphLayer.h"
 #include "Program.h"
-#include "CartographicTransformation.h"
-
+#include "advection/UVGrid.h"
+#include "advection/RK4AdvectionKernel.h"
 
 using namespace std;
 
+#define DT 60 * 60 // 60 sec/min * 60 mins
 
 int main() {
-  auto l = new LGlyphLayer();
-  l->spoofPoints();
+  cout << "reading data..." << endl;
+  shared_ptr<UVGrid> uvGrid = std::make_shared<UVGrid>();
+  auto kernelRK4 = make_unique<RK4AdvectionKernel>(uvGrid);
+  cout << "Starting vtk..." << endl;
 
+  auto l = new LGlyphLayer(uvGrid, std::move(kernelRK4));
 
-  Program *program = new Program();
+  Program *program = new Program(DT);
   program->addLayer(new BackgroundImage("../../../../data/map_661-661.png"));
-  program->addLayer(new EGlyphLayer());
+  program->addLayer(new EGlyphLayer(uvGrid));
   program->addLayer(l);
 
   program->render();