Merge pull request #30 from MakeNEnjoy/robin-projection

Robin projection

particle-track-and-trace/src/CMakeLists.txt

@@ -8,6 +8,7 @@ set(CMAKE_CXX_STANDARD_REQUIRED ON)
 set(CMAKE_EXPORT_COMPILE_COMMANDS ON)
 
 find_package(VTK COMPONENTS
+  GeovisCore
   CommonColor
   CommonColor
   CommonCore

particle-track-and-trace/src/CartographicTransformation.cpp

@@ -2,6 +2,9 @@
 #include <vtkMatrix4x4.h>
 #include <vtkTransform.h>
 #include <vtkTransformFilter.h>
+#include <vtkGeoProjection.h>
+#include <vtkGeoTransform.h>
+#include <vtkGeneralTransform.h>
 
 vtkSmartPointer<vtkCamera> createNormalisedCamera() {
   vtkSmartPointer<vtkCamera> camera = vtkSmartPointer<vtkCamera>::New();
@@ -15,32 +18,75 @@ vtkSmartPointer<vtkCamera> createNormalisedCamera() {
   return camera;
 }
 
-vtkSmartPointer<vtkMatrix4x4> getCartographicTransformMatrix(const std::shared_ptr<UVGrid> uvGrid) {
+// Assumes Normalised camera is used
+vtkSmartPointer<vtkTransformFilter> createCartographicTransformFilter(const std::shared_ptr<UVGrid> uvGrid) {
+  auto proj = vtkSmartPointer<vtkGeoProjection>::New();
+  proj->SetName("merc");
+
+  auto geoTransform = vtkSmartPointer<vtkGeoTransform>::New();
+  geoTransform->SetDestinationProjection(proj);
+
   const double XMin = uvGrid->lonMin();
   const double XMax = uvGrid->lonMax();
   const double YMin = uvGrid->latMin();
   const double YMax = uvGrid->latMax();
 
-    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
-    };
+  double bottomLeft[3] = {XMin, YMin, 0};
+  double topRight[3] = {XMax, YMax, 0};
+  geoTransform->TransformPoint(bottomLeft, bottomLeft);
+  geoTransform->TransformPoint(topRight, topRight);
 
-    auto matrix = vtkSmartPointer<vtkMatrix4x4>::New();
-    matrix->DeepCopy(eyeTransform);
-    return matrix;
-}
+  double width = topRight[0] - bottomLeft[0];
+  double height = topRight[1] - bottomLeft[1];
 
-// Assumes Normalised camera is used
-vtkSmartPointer<vtkTransformFilter> createCartographicTransformFilter(const std::shared_ptr<UVGrid> uvGrid) {
-    vtkNew<vtkTransform> transform;
+  auto scaleIntoNormalisedSpace = vtkSmartPointer<vtkTransform>::New();
+  scaleIntoNormalisedSpace->Scale(2 / (width), 2 / (height), 1);
+  scaleIntoNormalisedSpace->Translate(-(bottomLeft[0] + topRight[0]) / 2, -(bottomLeft[1] + topRight[1]) / 2, 0);
 
-    transform->SetMatrix(getCartographicTransformMatrix(uvGrid));
+  auto totalProjection = vtkSmartPointer<vtkGeneralTransform>::New();
+  totalProjection->Identity();
+  totalProjection->Concatenate(scaleIntoNormalisedSpace);
+  totalProjection->Concatenate(geoTransform);
 
   vtkSmartPointer<vtkTransformFilter> transformFilter = vtkSmartPointer<vtkTransformFilter>::New();
-    transformFilter->SetTransform(transform);
+  transformFilter->SetTransform(totalProjection);
+
+  return transformFilter;
+}
+
+vtkSmartPointer<vtkTransformFilter> createInverseCartographicTransformFilter(const std::shared_ptr<UVGrid> uvGrid) {
+  auto proj = vtkSmartPointer<vtkGeoProjection>::New();
+  proj->SetName("merc");
+
+  auto geoTransform = vtkSmartPointer<vtkGeoTransform>::New();
+  geoTransform->SetDestinationProjection(proj);
+
+  const double XMin = uvGrid->lonMin();
+  const double XMax = uvGrid->lonMax();
+  const double YMin = uvGrid->latMin();
+  const double YMax = uvGrid->latMax();
+
+  double bottomLeft[3] = {XMin, YMin, 0};
+  double topRight[3] = {XMax, YMax, 0};
+  geoTransform->TransformPoint(bottomLeft, bottomLeft);
+  geoTransform->TransformPoint(topRight, topRight);
+  geoTransform->Inverse();
+
+  double width = topRight[0] - bottomLeft[0];
+  double height = topRight[1] - bottomLeft[1];
+
+  auto scaleIntoNormalisedSpace = vtkSmartPointer<vtkTransform>::New();
+  scaleIntoNormalisedSpace->Scale(2 / (width), 2 / (height), 1);
+  scaleIntoNormalisedSpace->Translate(-(bottomLeft[0] + topRight[0]) / 2, -(bottomLeft[1] + topRight[1]) / 2, 0);
+  scaleIntoNormalisedSpace->Inverse();
+
+  auto totalProjection = vtkSmartPointer<vtkGeneralTransform>::New();
+  totalProjection->Identity();
+  totalProjection->Concatenate(geoTransform);
+  totalProjection->Concatenate(scaleIntoNormalisedSpace);
+
+  vtkSmartPointer<vtkTransformFilter> transformFilter = vtkSmartPointer<vtkTransformFilter>::New();
+  transformFilter->SetTransform(totalProjection);
 
   return transformFilter;
 }

particle-track-and-trace/src/CartographicTransformation.h

@@ -15,17 +15,20 @@
 vtkSmartPointer<vtkCamera> createNormalisedCamera();
 
 /**
- * Constructs a 4x4 projection matrix that maps homogenious (longitude, latitude, 0, 1) points
- * to the normalised space.
- * 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(const std::shared_ptr<UVGrid> uvGrid);
-
-/**
- * Convenience function that converts the 4x4 projection matrix into a vtkTransformFilter
+ *  Creates a mercator projection and then scales into normalised space.
+ *  @param uvGrid grid that contains latitudes and longitudes
  *  @return pointer to transform filter
  */
 vtkSmartPointer<vtkTransformFilter> createCartographicTransformFilter(const std::shared_ptr<UVGrid> uvGrid);
+
+/**
+ * Creates an inverse mercator projection and then scales into normalised space.
+ *
+ * Note that it is not possible to call the Inverse() member function on the transform of the
+ * createCartographicTransformFilter. It produces the wrong output,
+ * I believe that this is a bug in VTK, I might make an issue about this sometime.
+ * @param uvGrid grid that contains latitudes and longitudes
+ * @return pointer to transform filter
+ */
+vtkSmartPointer<vtkTransformFilter> createInverseCartographicTransformFilter(const std::shared_ptr<UVGrid> uvGrid);
 #endif //NORMALISEDCARTOGRAPHICCAMERA_H

particle-track-and-trace/src/commands/SpawnPointCallback.cpp

@@ -17,7 +17,6 @@ void convertDisplayToWorld(vtkRenderer *renderer, int x, int y, double *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);
 
@@ -39,7 +38,8 @@ void SpawnPointCallback::Execute(vtkObject *caller, unsigned long evId, void *ca
     ren->SetDisplayPoint(displayPos);
     ren->DisplayToWorld();
     ren->GetWorldPoint(worldPos);
-    inverseCartographicProjection->MultiplyPoint(worldPos, worldPos);
+    cout << "clicked on " << worldPos[1] << ", " << worldPos[0] << endl;
+    inverseCartographicProjection->TransformPoint(worldPos, worldPos);
 
     points->InsertNextPoint(worldPos[0], worldPos[1], 0);
     this->particlesBeached->InsertNextValue(0);
@@ -75,8 +75,7 @@ void SpawnPointCallback::setRen(const vtkSmartPointer<vtkRenderer> &ren) {
 
 void SpawnPointCallback::setUVGrid(const std::shared_ptr<UVGrid> &uvGrid) {
   this->uvGrid = uvGrid;
-  inverseCartographicProjection = getCartographicTransformMatrix(uvGrid);
-  inverseCartographicProjection->Invert();
+  inverseCartographicProjection = createInverseCartographicTransformFilter(uvGrid)->GetTransform();
 }
 
 void SpawnPointCallback::setBeached(const vtkSmartPointer<vtkIntArray> &ints) {

particle-track-and-trace/src/commands/SpawnPointCallback.h

@@ -1,13 +1,13 @@
 #ifndef SPAWNPOINTCALLBACK_H
 #define SPAWNPOINTCALLBACK_H
 
-
 #include <memory>
 #include <vtkCallbackCommand.h>
 #include <vtkRenderWindowInteractor.h>
 #include <vtkPoints.h>
 #include <vtkPolyData.h>
-#include <vtkMatrix4x4.h>
+#include <vtkAbstractTransform.h>
+
 #include "../advection/UVGrid.h"
 
 class SpawnPointCallback : public vtkCallbackCommand {
@@ -33,7 +33,7 @@ private:
   vtkSmartPointer<vtkRenderer> ren;
   vtkSmartPointer<vtkIntArray> particlesBeached;
   std::shared_ptr<UVGrid> uvGrid;
-  vtkSmartPointer<vtkMatrix4x4> inverseCartographicProjection;
+  vtkSmartPointer<vtkAbstractTransform> inverseCartographicProjection;
 
   void Execute(vtkObject *caller, unsigned long evId, void *callData) override;
 

particle-track-and-trace/src/layers/EGlyphLayer.cpp

@@ -12,10 +12,13 @@
 #include <vtkProperty2D.h>
 #include <vtkVertexGlyphFilter.h>
 #include <vtkArrowSource.h>
+
 #include "../CartographicTransformation.h"
 #include "../advection/readdata.h"
 #include "../advection/interpolate.h"
 
+#include "vtkTransform.h"
+
 using namespace std;
 
 EGlyphLayer::EGlyphLayer(std::shared_ptr<UVGrid> uvGrid) {
@@ -41,25 +44,29 @@ void EGlyphLayer::readCoordinates() {
   this->direction->SetNumberOfTuples(numLats * numLons);
   points->Allocate(numLats * numLons);
 
+  vtkSmartPointer<vtkTransformFilter> filter = createCartographicTransformFilter(uvGrid);
+  auto transform = filter->GetTransform();
+
   int i = 0;
   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);
+      direction->SetTuple3(i, u/2, v/2, 0);
+      // pre-transform the points, so we don't have to include the cartographic transform while running the program.
+      double out[3] = {lon, lat, 0};
+      transform->TransformPoint(out, out);
+      points->InsertPoint(i++, out[0], out[1], 0);
       lonIndex++;
     }
     latIndex++;
   }
+
   this->data->SetPoints(points);
   this->data->GetPointData()->AddArray(this->direction);
   this->data->GetPointData()->SetActiveVectors("direction");
 
-  vtkSmartPointer<vtkTransformFilter> transformFilter = createCartographicTransformFilter(uvGrid);
-  transformFilter->SetInputData(data);
-
   vtkNew<vtkGlyphSource2D> arrowSource;
   arrowSource->SetGlyphTypeToArrow();
   arrowSource->SetScale(0.2); //TODO: set this properly
@@ -67,7 +74,7 @@ void EGlyphLayer::readCoordinates() {
 
   vtkNew<vtkGlyph2D> glyph2D;
   glyph2D->SetSourceConnection(arrowSource->GetOutputPort());
-  glyph2D->SetInputConnection(transformFilter->GetOutputPort());
+  glyph2D->SetInputData(data);
   glyph2D->OrientOn();
   glyph2D->ClampingOn();
   glyph2D->SetScaleModeToScaleByVector();
@@ -94,7 +101,7 @@ void EGlyphLayer::updateData(int t) {
     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);
+      this->direction->SetTuple3(i, u/2, v/2, 0);
       i++;
     }
   }

particle-track-and-trace/src/layers/LGlyphLayer.cpp

@@ -95,7 +95,6 @@ LGlyphLayer::LGlyphLayer(std::shared_ptr<UVGrid> uvGrid, std::unique_ptr<Advecti
   this->ren->AddActor(actor);
 }
 
-
 // creates a few points so we can test the updateData function
 void LGlyphLayer::spoofPoints() {
     // auto id =this->points->InsertNextPoint(6.532949683882039, 53.24308582564463, 0); // Coordinates of Zernike

particle-track-and-trace/src/layers/LGlyphLayer.h

@@ -25,11 +25,13 @@ public:
     */
   LGlyphLayer(std::shared_ptr<UVGrid> uvGrid, std::unique_ptr<AdvectionKernel> advectionKernel);
 
-  /** This function spoofs a few points in the dataset. Mostly used for testing.
+  /**
+   * 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.
+  /**
+   * 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;

particle-track-and-trace/src/main.cpp

@@ -29,6 +29,7 @@ int main() {
   cout << "Starting vtk..." << endl;
 
   auto l = new LGlyphLayer(uvGrid, std::move(kernelRK4BoundaryChecked));
+  l->spoofPoints();
 
   unique_ptr<Program> program = make_unique<Program>(DT);
   program->addLayer(new BackgroundImage("../../../../data/map_661-661.png"));