kind of working but where did the euler layer go?

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,93 @@ vtkSmartPointer<vtkCamera> createNormalisedCamera() {
     return camera;
 }
 
-vtkSmartPointer<vtkMatrix4x4> getCartographicTransformMatrix(const std::shared_ptr<UVGrid> uvGrid) {
+//vtkSmartPointer<vtkMatrix4x4> getCartographicTransformMatrix(const std::shared_ptr<UVGrid> uvGrid) {
-    const double XMin = uvGrid->lons.front();
+//    const double XMin = uvGrid->lons.front();
-    const double XMax = uvGrid->lons.back();
+//    const double XMax = uvGrid->lons.back();
-    const double YMin = uvGrid->lats.front();
+//    const double YMin = uvGrid->lats.front();
-    const double YMax = uvGrid->lats.back();
+//    const double YMax = uvGrid->lats.back();
-
+//
-    double eyeTransform[] = {
+//    double eyeTransform[] = {
-            2/(XMax-XMin), 0, 0, -(XMax+XMin)/(XMax-XMin),
+//            2/(XMax-XMin), 0, 0, -(XMax+XMin)/(XMax-XMin),
-            0, 2/(YMax-YMin), 0, -(YMax+YMin)/(YMax-YMin),
+//            0, 2/(YMax-YMin), 0, -(YMax+YMin)/(YMax-YMin),
-            0, 0, 1, 0,
+//            0, 0, 1, 0,
-            0, 0, 0, 1
+//            0, 0, 0, 1
-    };
+//    };
-
+//
-    auto matrix = vtkSmartPointer<vtkMatrix4x4>::New();
+//    auto matrix = vtkSmartPointer<vtkMatrix4x4>::New();
-    matrix->DeepCopy(eyeTransform);
+//    matrix->DeepCopy(eyeTransform);
-    return matrix;
+//    return matrix;
-}
+//}
 
 // Assumes Normalised camera is used
 vtkSmartPointer<vtkTransformFilter> createCartographicTransformFilter(const std::shared_ptr<UVGrid> uvGrid) {
-    vtkNew<vtkTransform> transform;
+  auto proj = vtkSmartPointer<vtkGeoProjection>::New();
+  proj->SetName("merc");
 
-    transform->SetMatrix(getCartographicTransformMatrix(uvGrid));
+  auto geoTransform = vtkSmartPointer<vtkGeoTransform>::New();
+  geoTransform->SetDestinationProjection(proj);
 
-    vtkSmartPointer<vtkTransformFilter> transformFilter = vtkSmartPointer<vtkTransformFilter>::New();
+  const double XMin = -15.875;
-    transformFilter->SetTransform(transform);
+  const double XMax = 12.875;
+  const double YMin = 46.125;
+  const double YMax = 62.625;
 
-    return transformFilter;
+  double bottomLeft[3] = {XMin, YMin, 0};
+  double topRight[3] = {XMax, YMax, 0};
+  geoTransform->TransformPoint(bottomLeft, bottomLeft);
+  geoTransform->TransformPoint(topRight, topRight);
+
+  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);
+
+  auto totalProjection = vtkSmartPointer<vtkGeneralTransform>::New();
+  totalProjection->Identity();
+  totalProjection->Concatenate(scaleIntoNormalisedSpace);
+  totalProjection->Concatenate(geoTransform);
+
+  vtkSmartPointer<vtkTransformFilter> transformFilter = vtkSmartPointer<vtkTransformFilter>::New();
+  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 = -15.875;
+  const double XMax = 12.875;
+  const double YMin = 46.125;
+  const double YMax = 62.625;
+
+  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

@@ -28,4 +28,6 @@ vtkSmartPointer<vtkMatrix4x4> getCartographicTransformMatrix(const std::shared_p
  * @return pointer to transform filter
  */
 vtkSmartPointer<vtkTransformFilter> createCartographicTransformFilter(const std::shared_ptr<UVGrid> uvGrid);
+
+vtkSmartPointer<vtkTransformFilter> createInverseCartographicTransformFilter(const std::shared_ptr<UVGrid> uvGrid);
 #endif //NORMALISEDCARTOGRAPHICCAMERA_H

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

@@ -38,8 +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;
-    // cout << "clicked on lon = " << worldPos[0] << " and lat = " << worldPos[1] << endl;
+    inverseCartographicProjection->TransformPoint(worldPos, worldPos);
 
     vtkIdType id = points->InsertNextPoint(worldPos[0], worldPos[1], 0);
     data->SetPoints(points);
@@ -77,6 +77,5 @@ void SpawnPointCallback::setRen(const vtkSmartPointer<vtkRenderer> &ren) {
 
 void SpawnPointCallback::setUVGrid(const std::shared_ptr<UVGrid> &uvGrid) {
   this->uvGrid = uvGrid;
-  inverseCartographicProjection = getCartographicTransformMatrix(uvGrid);
+  inverseCartographicProjection = createInverseCartographicTransformFilter(uvGrid)->GetTransform();
-  inverseCartographicProjection->Invert();
 }

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 {
@@ -30,7 +30,7 @@ private:
   vtkSmartPointer<vtkPoints> points;
   vtkSmartPointer<vtkRenderer> ren;
   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/LGlyphLayer.cpp

@@ -68,17 +68,11 @@ LGlyphLayer::LGlyphLayer(std::shared_ptr<UVGrid> uvGrid, std::unique_ptr<Advecti
 
 // 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(6.532949683882039, 53.24308582564463, 0); // Coordinates of Zernike
+  this->points->InsertNextPoint(5.315307819255385, 60.40001057122271, 0); // Coordinates of Bergen
-    // 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.646210231365825, 46.52346296009023, 0); // Coordinates of Lausanne
+  this->points->InsertNextPoint(-6.553894313570932, 62.39522131195857,
-    // this->points->InsertNextPoint(-6.553894313570932, 62.39522131195857, 0); // Coordinates of the top of the Faroe islands
+                                0); // Coordinates of the top of the Faroe islands
-
-  for (int i=0; i < 330; i+=5) {
-    for (int j=0; j < 330; j+=5) {
-    this->points->InsertNextPoint(-15.875+(12.875+15.875)/330*j, 46.125+(62.625-46.125)/330*i, 0);
-    }
-  }
 
   this->points->Modified();
 }
@@ -89,7 +83,7 @@ void LGlyphLayer::updateData(int t) {
   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);
+      std::tie(point[1], point[0]) = advector->advect(t, point[1], point[0], (t - lastT) / SUPERSAMPLINGRATE);
     }
     this->points->SetPoint(n, point[0], point[1], 0);
   }

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

@@ -17,6 +17,12 @@
 #include "advection/kernel/RK4AdvectionKernel.h"
 #include "advection/kernel/SnapBoundaryConditionKernel.h"
 
+#include <vtkTransform.h>
+#include <vtkTransformFilter.h>
+#include <vtkGeoProjection.h>
+#include <vtkGeoTransform.h>
+#include <vtkGeneralTransform.h>
+
 using namespace std;
 
 #define DT 60 * 60 // 60 sec/min * 60 mins
@@ -29,6 +35,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"));
@@ -37,6 +44,47 @@ int main() {
 
   program->render();
 
+//  auto proj = vtkSmartPointer<vtkGeoProjection>::New(); proj->SetName("merc"); auto geoTransform = vtkSmartPointer<vtkGeoTransform>::New(); geoTransform->SetDestinationProjection(proj);
+//  const double XMin = -15.875;
+//  const double XMax = 12.875;
+//  const double YMin = 46.125;
+//  const double YMax = 62.625;
+//
+//  double bottomLeft[3] = {XMin, YMin, 0};
+//  double topRight[3] = {XMax, YMax, 0};
+//  geoTransform->TransformPoint(bottomLeft, bottomLeft);
+//  geoTransform->TransformPoint(topRight, topRight);
+//
+//  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);
+//
+//  auto totalProjection = vtkSmartPointer<vtkGeneralTransform>::New();
+//  totalProjection->PostMultiply();
+//  totalProjection->Identity();
+//  totalProjection->Concatenate(geoTransform);
+//  totalProjection->Concatenate(scaleIntoNormalisedSpace);
+//
+//  double in[3] = {4.846871030623073, 52.364810061968335, 0};
+//  geoTransform->TransformPoint(in, in);
+//  cout << "in[3] = {" << in[0] << "," << in[1] << "," << in[2] << "}" << endl;
+//  scaleIntoNormalisedSpace->TransformPoint(in, in);
+//  cout << "in[3] = {" << in[0] << "," << in[1] << "," << in[2] << "}" << endl;
+//  scaleIntoNormalisedSpace->Inverse();
+//  scaleIntoNormalisedSpace->TransformPoint(in, in);
+//  cout << "in[3] = {" << in[0] << "," << in[1] << "," << in[2] << "}" << endl;
+//  geoTransform->Inverse();
+//  geoTransform->TransformPoint(in, in);
+//  cout << "in[3] = {" << in[0] << "," << in[1] << "," << in[2] << "}" << endl;
+////  totalProjection->TransformPoint(in, in);
+////  cout << "in[3] = {" << in[0] << "," << in[1] << "," << in[2] << "}" << endl;
+////  totalProjection->Inverse();
+////  totalProjection->TransformPoint(in, in);
+////  cout << "in[3] = {" << in[0] << "," << in[1] << "," << in[2] << "}" << endl;
+
   return EXIT_SUCCESS;
 }