Merge pull request #15 from MakeNEnjoy/djairo-vtk-test

djairo vtk-test

vtk/.gitignore

@@ -0,0 +1,7 @@
+.DS_Store
+src/.DS_Store
+src/.cache
+src/build
+.idea
+src/cmake-build-debug
+compile_commands.json

vtk/README.md

@@ -0,0 +1,34 @@
+## Vtk
+This folder contains the Vtk program which actually displays the simulated data. The code is driven by the `Program` class, which contains the upper level of the vtk pipeline: the class has attributes for a vtkRenderWindow and vtkRenderWindowInteractor. vtkRenderers are managed through an abstract `Layer` class, which the program keeps track of trough a vector attribute.
+
+Each layer implementation contains and manages one vtkRenderer, this includes managing which layer of the vtkrenderwindow ths layer renders to. Currently implemented are three such layers:
+ * the `BackgroundImage` class reads in image data and displays this to the screen on the 0th layer - the background.
+ * the `EGlyphLayer` class renders a visualization of the Eulerian flow-velocities as a grid of arrow-glyphs (in which the direction and length of the glyph represents the direction and strength of the velocity at that point). Right now it spoofs the data for these glyphs, but this class will interface with the code for reading h5 data to accurately display the velocities at a given timestamp.
+ * the `LGlyphLayer` class renders a given set of particles as circular glyphs. These particles are advected according to an advection function, which in this implementation is spoofed. Like the EglyphLayer class, this layer will interact with the code for advecting particles according to the actual dataset, to accurately simulate its particles.
+
+The `LGlyphLayer` deserves some more explanation, as it depends on the `SpawnpointCallback` class to place particles in its dataset. The `SpawnpointCallback` makes use of the vtkCallbackCommand class and the vtk observer pattern to create new particles on mouseclick. It does so through a shared reference to the LGlyphLayer's `data` and `points` attributes, which the SpawnpointCallback then edits directlr. 
+
+The program also adds a second observer to the vtk pattern through the `TimerCallbackCommand`. This class subscribes to a vtkTimerEvent to manage the simulation of the program. To this end the TimerCallbackCommand has attributes for a timestep (dt) and current time (time). On every callback, the current time is updated according to the dt attribute, and this change is propagated to the layers containing the data by use of the program and layer's `updateData()` functions.
+
+
+## Location of data
+The data path is hardcoded such that the following tree structure is assumed:
+```
+data/
+  grid.h5
+  hydrodynamic_U.h5
+  hydrodynamic_V.h5
+interactive-track-and-trace/
+  opening-hdf5/
+     ...
+```
+
+## Compiling
+Let the current directory be the `src` directory. Run:
+```shell
+mkdir build
+cd build
+cmake ..
+make
+```
+

vtk/src/CMakeLists.txt

@@ -0,0 +1,76 @@
+cmake_minimum_required(VERSION 3.12 FATAL_ERROR)
+
+project(VtkBase)
+
+set(CMAKE_EXPORT_COMPILE_COMMANDS ON)
+
+find_package(VTK COMPONENTS 
+  CommonColor
+  CommonColor
+  CommonCore
+  CommonDataModel
+  FiltersGeneral
+  FiltersGeometry
+  FiltersProgrammable
+  FiltersSources
+  ImagingSources
+  InteractionStyle
+  IOImage
+  RenderingContextOpenGL2
+  RenderingCore
+  RenderingCore
+  RenderingFreeType
+  RenderingGL2PSOpenGL2
+  RenderingOpenGL2)
+
+
+if (NOT VTK_FOUND)
+  message(FATAL_ERROR "VtkBase: Unable to find the VTK build folder.")
+endif()
+
+# netcdf setup
+find_package(netCDF REQUIRED)
+
+add_executable(VtkBase MACOSX_BUNDLE main.cpp
+        layers/BackgroundImage.cpp
+        layers/BackgroundImage.h
+        layers/EGlyphLayer.cpp
+        layers/EGlyphLayer.h
+        layers/Layer.cpp
+        layers/Layer.h
+        layers/LGlyphLayer.cpp
+        layers/LGlyphLayer.h
+        Program.cpp
+        Program.h
+  commands/TimerCallbackCommand.h
+  commands/TimerCallbackCommand.cpp
+        commands/SpawnPointCallback.h
+        commands/SpawnPointCallback.cpp
+        CartographicTransformation.cpp
+)
+
+execute_process(
+    COMMAND nc-config --includedir
+    OUTPUT_VARIABLE NETCDF_INCLUDE_DIR
+    OUTPUT_STRIP_TRAILING_WHITESPACE
+)
+
+execute_process(
+    COMMAND ncxx4-config --libdir
+    OUTPUT_VARIABLE NETCDFCXX_LIB_DIR
+    OUTPUT_STRIP_TRAILING_WHITESPACE
+)
+
+target_include_directories(VtkBase PUBLIC ${netCDF_INCLUDE_DIR})
+
+find_library(NETCDF_LIB NAMES netcdf-cxx4 netcdf_c++4 PATHS ${NETCDFCXX_LIB_DIR} NO_DEFAULT_PATH)
+
+# Prevent a "command line is too long" failure in Windows.
+set(CMAKE_NINJA_FORCE_RESPONSE_FILE "ON" CACHE BOOL "Force Ninja to use response files.")
+target_link_libraries(VtkBase ${NETCDF_LIB} ${VTK_LIBRARIES})
+
+# vtk_module_autoinit is needed
+vtk_module_autoinit(
+  TARGETS VtkBase
+  MODULES ${VTK_LIBRARIES}
+)

vtk/src/CartographicTransformation.cpp

@@ -0,0 +1,46 @@
+#include "CartographicTransformation.h"
+#include <vtkMatrix4x4.h>
+#include <vtkTransform.h>
+#include <vtkTransformFilter.h>
+
+vtkSmartPointer<vtkCamera> createNormalisedCamera() {
+    vtkSmartPointer<vtkCamera> camera = vtkSmartPointer<vtkCamera>::New();
+    camera->ParallelProjectionOn();  // Enable parallel projection
+
+    camera->SetPosition(0, 0, 1000);  // Place the camera above the center
+    camera->SetFocalPoint(0, 0, 0);  // Look at the center
+    camera->SetViewUp(0, 1, 0); // Set the up vector to be along the Y-axis
+    camera->SetParallelScale(1); // x,y in [-1, 1]
+
+    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;
+
+    double eyeTransform[] = {
+            2/(XMax-XMin), 0, 0, -(XMax+XMin)/(XMax-XMin),
+            0, 2/(YMax-YMin), 0, -(YMax+YMin)/(YMax-YMin),
+            0, 0, 1, 0,
+            0, 0, 0, 1
+    };
+
+    auto matrix = vtkSmartPointer<vtkMatrix4x4>::New();
+    matrix->DeepCopy(eyeTransform);
+    return matrix;
+}
+
+// Assumes Normalised camera is used
+vtkSmartPointer<vtkTransformFilter> createCartographicTransformFilter() {
+    vtkNew<vtkTransform> transform;
+
+    transform->SetMatrix(getCartographicTransformMatrix());
+
+    vtkSmartPointer<vtkTransformFilter> transformFilter = vtkSmartPointer<vtkTransformFilter>::New();
+    transformFilter->SetTransform(transform);
+
+    return transformFilter;
+}

vtk/src/CartographicTransformation.h

@@ -0,0 +1,30 @@
+#include <vtkCamera.h>
+#include <vtkTransformFilter.h>
+
+#ifndef VTKBASE_NORMALISEDCARTOGRAPHICCAMERA_H
+#define VTKBASE_NORMALISEDCARTOGRAPHICCAMERA_H
+
+#endif //VTKBASE_NORMALISEDCARTOGRAPHICCAMERA_H
+
+/**
+ * Constructs a orthographically projected camera that looks at the square x,y in [-1, 1] with z = 0 and w = 1.
+ * The space [-1,1] x [-1,1] x {0} will be referred to as the normalised space.
+ * @return pointer to camera
+ */
+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();
+
+/**
+ * Convenience function that converts the 4x4 projection matrix into a vtkTransformFilter
+ * @return pointer to transform filter
+ */
+vtkSmartPointer<vtkTransformFilter> createCartographicTransformFilter();

vtk/src/Program.cpp

@@ -0,0 +1,86 @@
+#include <vtkRenderWindow.h>
+#include <vtkPointData.h>
+#include <vtkDoubleArray.h>
+#include <vtkGlyphSource2D.h>
+#include <vtkRegularPolygonSource.h>
+#include <vtkGlyph2D.h>
+#include <vtkActor2D.h>
+#include <vtkNamedColors.h>
+#include <vtkPolyDataMapper2D.h>
+#include <vtkPolyDataMapper.h>
+#include <vtkProperty.h>
+#include <vtkProperty2D.h>
+#include <vtkVertexGlyphFilter.h>
+#include <netcdf>
+#include <vtkArrowSource.h>
+#include <vtkNew.h>
+#include <vtkCallbackCommand.h>
+#include <vtkInteractorStyleUser.h>
+
+#include "Program.h"
+#include "commands/TimerCallbackCommand.h"
+#include "commands/SpawnPointCallback.h"
+
+void Program::setWinProperties() {
+  this->win->SetWindowName("Simulation");
+  this->win->SetSize(661, 661);
+  this->win->SetDesiredUpdateRate(60);
+
+  this->interact->SetRenderWindow(this->win);
+  this->interact->Initialize();
+
+  vtkNew<vtkInteractorStyleUser> style;
+  interact->SetInteractorStyle(style);
+}
+
+void Program::setupTimer() {
+  auto callback = vtkSmartPointer<TimerCallbackCommand>::New(this);
+  callback->SetClientData(this);
+  this->interact->AddObserver(vtkCommand::TimerEvent, callback);
+  this->interact->CreateRepeatingTimer(17); // 60 fps == 1000 / 60 == 16.7 ms per frame
+}
+
+Program::Program() {
+  this->win = vtkSmartPointer<vtkRenderWindow>::New();
+  this->interact = vtkSmartPointer<vtkRenderWindowInteractor>::New();
+
+  this->win->SetNumberOfLayers(0);
+  setWinProperties();
+  setupTimer();
+}
+
+
+void Program::addLayer(Layer *layer) {
+  this->layers.push_back(layer);
+  this->win->AddRenderer(layer->getLayer());
+  this->win->SetNumberOfLayers(this->win->GetNumberOfLayers() + 1);
+}
+
+void Program::removeLayer(Layer *layer) {
+  this->win->RemoveRenderer(layer->getLayer());
+
+  auto it = std::find(this->layers.begin(), this->layers.end(), layer);
+  if (it != this->layers.end()) {
+    this->layers.erase(it);
+    this->win->SetNumberOfLayers(this->win->GetNumberOfLayers() - 1);
+  }
+}
+
+void Program::updateData(int t) {
+  win->Render();
+  for (Layer *l: layers) {
+    l->updateData(t);
+  }
+}
+
+void Program::setupInteractions() {
+  for (Layer *l: layers) {
+    l->addObservers(interact);
+  }
+}
+
+void Program::render() {
+  setupInteractions();
+  win->Render();
+  interact->Start();
+}

vtk/src/Program.h

@@ -0,0 +1,68 @@
+#ifndef PROGRAM_H
+#define PROGRAM_H
+
+#include <vtkRenderWindow.h>
+#include <vtkRenderWindowInteractor.h>
+#include <vtkRenderer.h>
+
+#include "layers/Layer.h"
+#include "commands/SpawnPointCallback.h"
+
+/** This class manages the upper levels of the vtk pipeline; it has attributes for the vtkrenderWindow and a vector of Layers to represent a variable number of vtkRenderers.
+  * It can also set up a vtkTimer by connecting an instance of TimerCallbackCommand with its contained vtkRenderWindowInteractor.
+  */ 
+class Program {
+private:
+  /** This attribute models a variable number of vtkRenderers, managed through the abstract Layer class.
+    */ 
+  std::vector<Layer *> layers;
+
+  /** The window this program's layers render to.
+    */ 
+  vtkSmartPointer<vtkRenderWindow> win;
+
+  /** The interactor through which the layers can interact with the window.
+    */ 
+  vtkSmartPointer<vtkRenderWindowInteractor> interact;
+
+  /** This function sets some default properties on the vtkRenderWindow. Extracted to its' own function to keep the constructor from becoming cluttered.
+    */ 
+  void setWinProperties();
+
+  /** This function sets up and connects a TimerCallbackCommand with the program.
+    */ 
+  void setupTimer();
+
+  void setupInteractions();
+
+public:
+  /** Constructor.
+    */ 
+  Program();
+
+  /** 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.
+    * @param layer : pointer to the layer to add.
+    */ 
+  void addLayer(Layer *layer);
+
+  /** This function removes a given layer from the vtkRenderWindow and layers vector.
+    * If the given layer is not actually in the program, nothing happens.
+    * @param layer : the layer to removeLayer
+    */
+  void removeLayer(Layer *layer);
+
+  /** This function updates the data for the associated layers to the given timestamp.
+    * Also updates the renderWindow.
+    * @param t : the timestamp to update the data to.
+    */
+  void updateData(int t);
+
+  /**
+   * This function renders the vtkRenderWindow for the first time.
+   * Only call this function once!
+   */
+  void render();
+};
+
+#endif

vtk/src/commands/SpawnPointCallback.cpp

@@ -0,0 +1,77 @@
+#include "SpawnPointCallback.h"
+
+#include <vtkVertex.h>
+#include <vtkRenderer.h>
+#include <vtkRenderWindowInteractor.h>
+#include <vtkSmartPointer.h>
+#include <vtkCommand.h>
+#include <vtkRenderWindow.h>
+
+#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 SpawnPointCallback::Execute(vtkObject *caller, unsigned long evId, void *callData) {
+    // Note the use of reinterpret_cast to cast the caller to the expected type.
+    auto interactor = reinterpret_cast<vtkRenderWindowInteractor *>(caller);
+
+    if (evId == vtkCommand::LeftButtonPressEvent) {
+        dragging = true;
+    }
+    if (evId == vtkCommand::LeftButtonReleaseEvent) {
+        dragging = false;
+    }
+    if (!dragging) {
+        return;
+    }
+
+    int x, y;
+    interactor->GetEventPosition(x, y);
+
+    double worldPos[4] = {2, 0 ,0, 0};
+    double displayPos[3] = {static_cast<double>(x), static_cast<double>(y), 0.0};
+    ren->SetDisplayPoint(displayPos);
+    ren->DisplayToWorld();
+    ren->GetWorldPoint(worldPos);
+    inverseCartographicProjection->MultiplyPoint(worldPos, worldPos);
+    cout << "clicked on lon = " << worldPos[0] << " and lat = " << worldPos[1] << endl;
+
+    vtkIdType id = points->InsertNextPoint(worldPos[0], worldPos[1], 0);
+    data->SetPoints(points);
+
+    vtkSmartPointer<vtkVertex> vertex = vtkSmartPointer<vtkVertex>::New();
+    vertex->GetPointIds()->SetId(0, id);
+
+    vtkSmartPointer<vtkCellArray> vertices = vtkSmartPointer<vtkCellArray>::New();
+    vertices->InsertNextCell(vertex);
+    data->SetVerts(vertices);
+//    data->Modified(); // These might be needed im not sure.
+//    ren->GetRenderWindow()->Render();
+}
+
+
+SpawnPointCallback::SpawnPointCallback() : data(nullptr), points(nullptr), inverseCartographicProjection(nullptr) {
+    inverseCartographicProjection = getCartographicTransformMatrix();
+    inverseCartographicProjection->Invert();
+}
+
+SpawnPointCallback *SpawnPointCallback::New() {
+    return new SpawnPointCallback;
+}
+
+void SpawnPointCallback::setData(const vtkSmartPointer<vtkPolyData> &data) {
+    this->data = data;
+}
+
+void SpawnPointCallback::setPoints(const vtkSmartPointer<vtkPoints> &points) {
+    this->points = points;
+}
+
+void SpawnPointCallback::setRen(const vtkSmartPointer<vtkRenderer> &ren) {
+    this->ren = ren;
+}

vtk/src/commands/SpawnPointCallback.h

@@ -0,0 +1,33 @@
+#ifndef VTKBASE_SPAWNPOINTCALLBACK_H
+#define VTKBASE_SPAWNPOINTCALLBACK_H
+
+
+#include <vtkCallbackCommand.h>
+#include <vtkRenderWindowInteractor.h>
+#include <vtkPoints.h>
+#include <vtkPolyData.h>
+#include <vtkMatrix4x4.h>
+
+class SpawnPointCallback : public vtkCallbackCommand {
+
+public:
+    static SpawnPointCallback *New();
+    SpawnPointCallback();
+
+    void setPoints(const vtkSmartPointer<vtkPoints> &points);
+
+    void setData(const vtkSmartPointer<vtkPolyData> &data);
+
+    void setRen(const vtkSmartPointer<vtkRenderer> &ren);
+private:
+    vtkSmartPointer<vtkPolyData> data;
+    vtkSmartPointer<vtkPoints> points;
+    vtkSmartPointer<vtkRenderer> ren;
+    vtkSmartPointer<vtkMatrix4x4> inverseCartographicProjection;
+
+    void Execute(vtkObject *caller, unsigned long evId, void *callData) override;
+    bool dragging = false;
+};
+
+
+#endif //VTKBASE_SPAWNPOINTCALLBACK_H

vtk/src/commands/TimerCallbackCommand.cpp

@@ -0,0 +1,28 @@
+#include "TimerCallbackCommand.h"
+#include "../Program.h"
+
+
+TimerCallbackCommand::TimerCallbackCommand() : dt(3600), maxTime(3600*24*365), time(0) {}
+
+TimerCallbackCommand* TimerCallbackCommand::New(Program *program) {
+  TimerCallbackCommand *cb = new TimerCallbackCommand();
+  cb->setProgram(program);
+  return cb;
+}
+
+void TimerCallbackCommand::Execute(vtkObject* caller, unsigned long eventId, void* vtkNotUsed(callData)) {
+  this->time += this->dt;
+
+  if (this->time >= this->maxTime) {
+    return;
+    // TODO: how do we deal with reaching the end of the simulated dataset? Do we just stop simulating, loop back around? What about the location of the particles in this case? Just some ideas to consider, but we should iron this out pretty soon.
+  }
+
+  this->program->updateData(this->time);
+}
+
+
+
+void TimerCallbackCommand::setProgram(Program *program) {
+  this->program = program;
+}

vtk/src/commands/TimerCallbackCommand.h

@@ -0,0 +1,22 @@
+#ifndef TIMERCALLBACKCOMMAND_H
+#define TIMERCALLBACKCOMMAND_H
+
+#include <vtkCallbackCommand.h>
+#include "../Program.h"
+
+class TimerCallbackCommand : public vtkCallbackCommand {
+public:
+  TimerCallbackCommand();
+  static TimerCallbackCommand* New(Program *program);
+  void Execute(vtkObject* caller, unsigned long eventId, void* vtkNotUsed(callData)) override;
+
+  void setProgram(Program *program);
+
+private:
+  int time;
+  int dt;
+  int maxTime;
+  Program *program;
+};
+
+#endif

vtk/src/layers/BackgroundImage.cpp

@@ -0,0 +1,62 @@
+#include "BackgroundImage.h"
+#include <vtkCamera.h>
+#include <vtkImageActor.h>
+#include <vtkImageData.h>
+#include <vtkImageReader2.h>
+
+using std::string;
+
+BackgroundImage::BackgroundImage(string imagePath) : imagePath(imagePath) {
+  this->ren = vtkSmartPointer<vtkRenderer>::New();
+  this->ren->SetLayer(0);
+  this->ren->InteractiveOff();
+  updateImage();
+}
+
+
+void BackgroundImage::updateImage() {
+  vtkSmartPointer<vtkImageData> imageData;
+
+  vtkSmartPointer<vtkImageReader2> imageReader;
+  
+  imageReader.TakeReference(this->readerFactory->CreateImageReader2(this->imagePath.c_str()));
+  imageReader->SetFileName(this->imagePath.c_str());
+  imageReader->Update();
+  imageData = imageReader->GetOutput();
+
+  vtkNew<vtkImageActor> imageActor;
+  imageActor->SetInputData(imageData);
+
+  this->ren->AddActor(imageActor);
+
+
+  // camera stuff
+  // essentially sets the camera to the middle of the background, and points it at the background
+  // TODO: extract this to its own function, separate from the background class.
+  double origin[3], spacing[3];
+  int extent[6];
+  imageData->GetOrigin(origin);
+  imageData->GetSpacing(spacing);
+  imageData->GetExtent(extent);
+
+  vtkCamera *camera = this->ren->GetActiveCamera();
+  camera->ParallelProjectionOn();
+
+  double xc = origin[0] + 0.5 * (extent[0] + extent[1]) * spacing[0];
+  double yc = origin[1] + 0.5 * (extent[2] + extent[3]) * spacing[1];
+  double yd = (extent[3] - extent[2] + 1) * spacing[1];
+  double d = camera->GetDistance();
+  camera->SetParallelScale(0.5 * yd);
+  camera->SetFocalPoint(xc, yc, 0.0);
+  camera->SetPosition(xc, yc, d);
+}
+
+
+string BackgroundImage::getImagePath() {
+  return this->imagePath;
+}
+
+void BackgroundImage::setImagePath(string imagePath) {
+  this->imagePath = imagePath;
+  updateImage();
+}

vtk/src/layers/BackgroundImage.h

@@ -0,0 +1,37 @@
+#ifndef BACKGROUND_H
+#define BACKGROUND_H
+
+#include "Layer.h"
+#include <vtkImageReader2Factory.h>
+
+/** Implements the Layer class for the case of a background image.
+  * Specifically, reads a backgroundImage given by the imagePath attribute and puts it on layer 0.
+  */
+class BackgroundImage : public Layer {
+private:
+  std::string imagePath;
+  vtkSmartPointer<vtkImageReader2Factory> readerFactory;
+
+  /** This private function updates the background image using the imagePath attribute.
+   */
+  void updateImage();
+
+
+public:
+  /** Constructor.
+    * @param imagePath : String to the path of the image to use as background.
+    */
+  BackgroundImage(std::string imagePath);
+
+  /** Getter.
+    * @return the imagePath attribute.
+    */
+  std::string getImagePath();
+
+  /** Setter. Can be used to change the background image 
+    * @param imagePath : String to the path of the new image to use.
+    */
+  void setImagePath(std::string imagePath);
+};
+
+#endif

vtk/src/layers/EGlyphLayer.cpp

@@ -0,0 +1,124 @@
+#include "EGlyphLayer.h"
+#include <vtkPointData.h>
+#include <vtkDoubleArray.h>
+#include <vtkGlyphSource2D.h>
+#include <vtkRegularPolygonSource.h>
+#include <vtkGlyph2D.h>
+#include <vtkActor2D.h>
+#include <vtkNamedColors.h>
+#include <vtkPolyDataMapper2D.h>
+#include <vtkPolyDataMapper.h>
+#include <vtkProperty.h>
+#include <vtkProperty2D.h>
+#include <vtkVertexGlyphFilter.h>
+#include <netcdf>
+#include <vtkArrowSource.h>
+#include "../CartographicTransformation.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() {
+  this->ren = vtkSmartPointer<vtkRenderer>::New();
+  this->ren->SetLayer(1);
+  this->ren->InteractiveOff();
+
+  this->data = vtkSmartPointer<vtkPolyData>::New();
+  this->direction = vtkSmartPointer<vtkDoubleArray>::New();
+  this->direction->SetName("direction");
+  readCoordinates();
+}
+
+
+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->direction->SetNumberOfComponents(3);
+  this->direction->SetNumberOfTuples(numLats*numLons); 
+  points->Allocate(numLats*numLons);
+
+  auto camera = createNormalisedCamera();
+  ren->SetActiveCamera(camera);
+
+  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
+      points->InsertPoint(i++, lon, lat, 0);
+      // see also https://vtk.org/doc/nightly/html/classvtkPolyDataMapper2D.html
+    }
+  }
+  this->data->SetPoints(points);
+  this->data->GetPointData()->AddArray(this->direction);
+  this->data->GetPointData()->SetActiveVectors("direction");
+
+  vtkSmartPointer<vtkTransformFilter> transformFilter = createCartographicTransformFilter();
+  transformFilter->SetInputData(data);
+
+  vtkNew<vtkGlyphSource2D> arrowSource;
+  arrowSource->SetGlyphTypeToArrow();
+  arrowSource->SetScale(0.2); //TODO: set this properly
+  arrowSource->Update();
+
+  vtkNew<vtkGlyph2D> glyph2D;
+  glyph2D->SetSourceConnection(arrowSource->GetOutputPort());
+  glyph2D->SetInputConnection(transformFilter->GetOutputPort());
+  glyph2D->OrientOn();
+  glyph2D->ClampingOn();
+  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()->SetOpacity(0.2);
+
+  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.
+  }
+  this->direction->Modified();
+}

vtk/src/layers/EGlyphLayer.h

@@ -0,0 +1,35 @@
+#ifndef EGLYPHLAYER_H
+#define EGLYPHLAYER_H
+
+#include "Layer.h"
+#include <vtkPolyData.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.
+  */
+class EGlyphLayer : public Layer {
+private:
+  vtkSmartPointer<vtkPolyData> data;
+  vtkSmartPointer<vtkDoubleArray> direction;
+  int numLats;
+  int numLons;
+
+  /** This private function sets up the initial coordinates for the glyphs in the dataset.
+    * It also reads some initial data to actually display.
+    */
+  void readCoordinates();
+
+public:
+  /** Constructor.
+    */
+  EGlyphLayer();
+
+  /** updates the glyphs to reflect the given timestamp in the dataset.
+    * @param t : the time at which to fetch the data.
+    */
+  void updateData(int t);
+
+};
+
+
+#endif

vtk/src/layers/LGlyphLayer.cpp

@@ -0,0 +1,105 @@
+#include "LGlyphLayer.h"
+#include "../commands/SpawnPointCallback.h"
+#include <vtkActor2D.h>
+#include <vtkGlyph2D.h>
+#include <vtkGlyphSource2D.h>
+#include <vtkNamedColors.h>
+#include <vtkPolyDataMapper2D.h>
+#include <vtkProperty2D.h>
+#include <vtkVertexGlyphFilter.h>
+#include <vtkInteractorStyle.h>
+#include <vtkInteractorStyleUser.h>
+#include <vtkTransform.h>
+#include <vtkTransformFilter.h>
+#include <vtkPolyDataMapper.h>
+#include <vtkRenderWindow.h>
+#include <vtkCamera.h>
+
+#include "../CartographicTransformation.h"
+
+
+vtkSmartPointer<SpawnPointCallback> LGlyphLayer::createSpawnPointCallback() {
+    auto newPointCallBack = vtkSmartPointer<SpawnPointCallback>::New();
+    newPointCallBack->setData(data);
+    newPointCallBack->setPoints(points);
+    newPointCallBack->setRen(ren);
+    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.
+// Another challenge is the concept of beaching; dead particles must not be included in the advect function call (wasted computations), but they should not be outright deleted from the vtkPoints either (we still want to display them). Working Solution: have another array of ints in the vtkPolyData, which tracks for how many calls of UpdateData a given particle has not had its position changed. If this int reaches some treshold (5? 10? 3? needs some testing), exclude the particle from the advect call.
+//
+// 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);
+
+    this->points = vtkSmartPointer<vtkPoints>::New();
+    this->data = vtkSmartPointer<vtkPolyData>::New();
+    this->data->SetPoints(this->points);
+
+    auto camera = createNormalisedCamera();
+    ren->SetActiveCamera(camera);
+
+    auto transform = createCartographicTransformFilter();
+
+    vtkSmartPointer<vtkTransformFilter> transformFilter = createCartographicTransformFilter();
+    transformFilter->SetInputData(data);
+
+    vtkNew<vtkGlyphSource2D> circleSource;
+    circleSource->SetGlyphTypeToCircle();
+    circleSource->SetScale(0.05);
+    circleSource->Update();
+
+    vtkNew<vtkGlyph2D> glyph2D;
+    glyph2D->SetSourceConnection(circleSource->GetOutputPort());
+    glyph2D->SetInputConnection(transformFilter->GetOutputPort());
+    glyph2D->SetColorModeToColorByScalar();
+    glyph2D->Update();
+
+    vtkNew<vtkPolyDataMapper> mapper;
+    mapper->SetInputConnection(glyph2D->GetOutputPort());
+    mapper->Update();
+
+    vtkNew<vtkActor> actor;
+    actor->SetMapper(mapper);
+
+    this->ren->AddActor(actor);
+}
+
+// creates a few points so we can test the updateData function
+void LGlyphLayer::spoofPoints() {
+    this->points->InsertNextPoint(-4.125, 61.375 , 0);
+    this->points->InsertNextPoint(6.532949683882039, 53.24308582564463, 0); // Coordinates of Zernike
+    this->points->InsertNextPoint(5.315307819255385, 60.40001057122271, 0); // Coordinates of Bergen
+    this->points->InsertNextPoint( 6.646210231365825, 46.52346296009023, 0); // Coordinates of Lausanne
+    this->points->InsertNextPoint(-6.553894313570932, 62.39522131195857, 0); // Coordinates of the top of the Faroe islands
+
+    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);
+    }
+    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);
+}

vtk/src/layers/LGlyphLayer.h

@@ -0,0 +1,38 @@
+#ifndef LGLYPHLAYER_H
+#define LGLYPHLAYER_H
+
+#include "Layer.h"
+#include "../commands/SpawnPointCallback.h"
+#include <vtkPolyData.h>
+#include <vtkInteractorStyle.h>
+
+/** Implements the Layer class for the case of a Lagrangian visualization.
+  * Specifically, this class models the Lagrangian particles in the simulation using the 'glyph' mark and 'transparency' channel to denote age.
+  */
+class LGlyphLayer : public Layer {
+private:
+  vtkSmartPointer<vtkPoints> points;
+  vtkSmartPointer<vtkPolyData> data;
+
+
+
+public:
+  /** Constructor.
+    */
+  LGlyphLayer();
+
+  /** This function spoofs a few points in the dataset. Mostly used for testing.
+    */
+  void spoofPoints();
+
+  /** updates the glyphs to reflect the given timestamp in the dataset.
+    * @param t : the time at which to fetch the data.
+    */
+  void updateData(int t) override;
+
+  vtkSmartPointer<SpawnPointCallback> createSpawnPointCallback();
+
+  void addObservers(vtkSmartPointer<vtkRenderWindowInteractor> interactor) override;
+};
+
+#endif

vtk/src/layers/Layer.cpp

@@ -0,0 +1,17 @@
+#include "Layer.h"
+#include <vtkRenderWindow.h>
+#include <vtkRenderWindowInteractor.h>
+
+using std::string;
+
+vtkSmartPointer<vtkRenderer> Layer::getLayer() {
+  return this->ren;
+}
+
+void Layer::updateData(int t) {
+  // By default, do nothing
+}
+
+void Layer::addObservers(vtkSmartPointer<vtkRenderWindowInteractor> interactor) {
+    // By default, do nothing
+}

vtk/src/layers/Layer.h

@@ -0,0 +1,32 @@
+#ifndef LAYER_H
+#define LAYER_H
+
+#include <vtkInteractorStyle.h>
+#include <vtkRenderer.h>
+
+/** This class represents one abstract layer to be rendered to VTK.
+  * It exists to manage multiple different layers under the Program class.
+  */
+class Layer {
+protected:
+  vtkSmartPointer<vtkRenderer> ren;
+
+public:
+  /** gets the vtkRenderer to assign it to the vtkRenderWindow of the program class.
+   * @return pointer to the vtkRenderer of this class.
+   */
+  virtual vtkSmartPointer<vtkRenderer> getLayer();
+
+
+  /** updates the data in the layer to reflect the given timestamp.
+    * @param t : the timestamp which the data should reflect.
+    */
+  virtual void updateData(int t);
+
+  /** Adds observers to the renderWindowinteractor within which this layer is active.
+    * @param interactor : pointer to the interactor that observers can be added to.
+    */
+  virtual void addObservers(vtkSmartPointer<vtkRenderWindowInteractor> interactor);
+};
+
+#endif

vtk/src/main.cpp

@@ -0,0 +1,32 @@
+#include <netcdf>
+#include <vtkActor2D.h>
+#include <vtkNamedColors.h>
+#include <vtkPoints.h>
+#include <vtkPolyData.h>
+#include <vtkPolyDataMapper2D.h>
+#include <vtkProperty2D.h>
+#include <vtkRenderer.h>
+#include <vtkVertexGlyphFilter.h>
+
+#include "layers/BackgroundImage.h"
+#include "layers/EGlyphLayer.h"
+#include "layers/LGlyphLayer.h"
+#include "Program.h"
+
+using namespace std;
+
+
+int main() {
+  auto l = new LGlyphLayer();
+  l->spoofPoints();
+
+  Program *program = new Program();
+  program->addLayer(new BackgroundImage("../../../../data/map_661-661.png"));
+  program->addLayer(new EGlyphLayer());
+  program->addLayer(l);
+
+  program->render();
+
+  return EXIT_SUCCESS;
+}
+

vtk/src/script.py

@@ -0,0 +1,226 @@
+#!/usr/bin/env python
+
+import collections
+import json
+import os
+import re
+from pathlib import Path
+
+
+def get_program_parameters(argv):
+    import argparse
+    description = 'Generate a find_package(VTK COMPONENTS ...) command for CMake.'
+    epilogue = '''
+Uses modules.json and your source files to generate a
+  find_package(VTK COMPONENTS ...) command listing all the vtk modules
+  needed by the C++ source and header files in your code.
+
+Paths for more than one source path can be specified.
+
+Note than if there are spaces in the paths, enclose the path in quotes.
+
+If it is unable to find modules for your headers then
+  a list of these, along with the files they are in, is produced
+  so you can manually add the corresponding modules or rebuild VTK
+  to include the missing modules.
+
+You will need to manually add any third-party modules
+   (if used) to the find_package command.
+    '''
+    parser = argparse.ArgumentParser(description=description, epilog=epilogue,
+                                     formatter_class=argparse.RawTextHelpFormatter)
+    parser.add_argument('json', default=['modules.json'], help='The path to the VTK JSON file (modules.json).')
+    parser.add_argument('sources', nargs='+', help='The path to the source files.')
+    parser.add_argument('-f', '--file', help='The file name to write the output too.')
+    args = parser.parse_args()
+    return args.json, args.sources, args.file
+
+
+class Patterns:
+    header_pattern = re.compile(r'^#include *[<\"](\S+)[>\"]')
+    vtk_include_pattern = re.compile(r'^(vtk\S+)')
+    vtk_qt_include_pattern = re.compile(r'^(QVTK\S+)')
+
+
+def get_headers_modules(json_data):
+    """
+    From the parsed JSON data file make a dictionary whose key is the
+     header filename and value is the module.
+    :param json_data: The parsed JSON file modules.json.
+    :return:
+    """
+
+    # The headers should be unique to a module, however we will not assume this.
+    res = collections.defaultdict(set)
+    for k, v in json_data['modules'].items():
+        if 'headers' in v:
+            for k1 in v['headers']:
+                res[k1].add(k)
+    return res
+
+
+def get_vtk_components(jpath, paths):
+    """
+    Get the VTK components
+    :param jpath: The path to the JSON file.
+    :param paths: The C++ file paths.
+    :return:
+    """
+
+    with open(jpath) as data_file:
+        json_data = json.load(data_file)
+    vtk_headers_modules = get_headers_modules(json_data)
+
+    modules = set()
+    inc_no_mod = set()
+    inc_no_mod_headers = collections.defaultdict(set)
+    mod_implements = collections.defaultdict(set)
+    headers = collections.defaultdict(set)
+
+    for path in paths:
+        if path.is_file():
+            content = path.read_text().split('\n')
+            for line in content:
+                m = Patterns.header_pattern.match(line.strip())
+                if m:
+                    # We have a header name, split it from its path (if the path exists).
+                    header_parts = os.path.split(m.group(1))
+                    m = Patterns.vtk_include_pattern.match(header_parts[1])
+                    if m:
+                        headers[m.group(1)].add(path)
+                        continue
+                    m = Patterns.vtk_qt_include_pattern.match(header_parts[1])
+                    if m:
+                        headers[m.group(1)].add(path)
+    for incl in headers:
+        if incl in vtk_headers_modules:
+            m = vtk_headers_modules[incl]
+            for v in m:
+                modules.add(v)
+        else:
+            inc_no_mod.add(incl)
+            inc_no_mod_headers[incl] = headers[incl]
+
+    if headers:
+        for m in modules:
+            if not json_data['modules'][m]['implementable']:
+                continue
+            for i in json_data['modules']:
+                if i in modules:
+                    continue
+                if m in json_data['modules'][i]['implements']:
+                    # Suggest module i since it implements m
+                    mod_implements[i].add(m)
+
+    return modules, mod_implements, inc_no_mod, inc_no_mod_headers
+
+
+def disp_components(modules, module_implements):
+    """
+    For the found modules display them in a form that the user can
+     copy/paste into their CMakeLists.txt file.
+    :param modules: The modules.
+    :param module_implements: Modules implementing other modules.
+    :return:
+    """
+    res = ['find_package(VTK\n COMPONENTS']
+    for m in sorted(modules):
+        res.append('    {:s}'.format(m.split('::')[1]))
+    if module_implements:
+        keys = sorted(module_implements)
+        max_width = len(max(keys, key=len).split('::')[1])
+        comments = [
+            '    #',
+            '    # These modules are suggested since they implement an existing module.',
+            '    # You may need to uncomment one or more of these.',
+            '    # If vtkRenderWindow is used and you want to use OpenGL,',
+            '    #   you also need the RenderingOpenGL2 module.',
+            '    # If vtkRenderWindowInteractor is used,',
+            '    #    uncomment RenderingUI and possibly InteractionStyle.',
+            '    # If text rendering is used, uncomment RenderingFreeType',
+            '    #'
+        ]
+        res.extend(comments)
+        for key in keys:
+            res.append(
+                f'    # {key.split("::")[1]:<{max_width}} # implements {", ".join(sorted(module_implements[key]))}')
+    res.append(')\n')
+
+    return res
+
+
+def disp_missing_components(inc_no_mod, inc_no_mod_headers):
+    """
+    Display the headers along with the missing VTK modules.
+
+    :param inc_no_mod: Missing modules.
+    :param inc_no_mod_headers: Headers with missing modules.
+    :return:
+    """
+    if inc_no_mod:
+        res = [''
+               '*' * 64,
+               'You will need to manually add the modules that',
+               '  use these headers to the find_package command.',
+               'These could be external modules not in the modules.json file.',
+               'Or you may need to rebuild VTK to include the missing modules.',
+               '',
+               'Here are the vtk headers and corresponding files:']
+        sinmd = sorted(inc_no_mod)
+        for i in sinmd:
+            sh = sorted(list(inc_no_mod_headers[i]))
+            res.append(f'in {i}:')
+            for j in sh:
+                res.append(f'   {j}')
+        res.append('*' * 64)
+
+        return res
+    else:
+        return None
+
+
+def main(json_path, src_paths, ofn):
+    jpath = Path(json_path)
+    if jpath.is_dir():
+        jpath = jpath / 'modules.json'
+    if not jpath.is_file():
+        print(f'Non existent JSON file: {jpath}')
+        return
+
+    paths = list()
+    valid_ext = ['.h', '.hxx', '.cxx', '.cpp', '.txx']
+    path_list = list()
+    for fn in src_paths:
+        path = Path(fn)
+        if path.is_file() and path.suffix in valid_ext:
+            paths.append(path)
+        elif path.is_dir():
+            for e in valid_ext:
+                path_list += list(Path(fn).rglob(f'*{e}'))
+            program_path = Path(__file__)
+            for path in path_list:
+                if path.resolve() != program_path.resolve():
+                    paths.append(path)
+        else:
+            print(f'Non existent path: {path}')
+
+    modules, mod_implements, inc_no_mod, inc_no_mod_headers = get_vtk_components(jpath, paths)
+
+    res = '\n'.join(disp_components(modules, mod_implements))
+    if inc_no_mod:
+        res += '\n'.join(disp_missing_components(inc_no_mod, inc_no_mod_headers))
+
+    if ofn:
+        path = Path(ofn)
+        if path.suffix == '':
+            path = Path(ofn).with_suffix('.txt')
+        path.write_text(res)
+    else:
+        print(res)
+
+
+if __name__ == '__main__':
+    import sys
+
+    json_paths, src_paths, ofn = get_program_parameters(sys.argv)
+    main(json_paths, src_paths, ofn)