diff --git a/particle-track-and-trace/src/CMakeLists.txt b/particle-track-and-trace/src/CMakeLists.txt index 9364e77..ac28701 100644 --- a/particle-track-and-trace/src/CMakeLists.txt +++ b/particle-track-and-trace/src/CMakeLists.txt @@ -7,7 +7,8 @@ set(CMAKE_CXX_STANDARD_REQUIRED ON) set(CMAKE_EXPORT_COMPILE_COMMANDS ON) -find_package(VTK COMPONENTS +find_package(VTK COMPONENTS + GeovisCore CommonColor CommonColor CommonCore diff --git a/particle-track-and-trace/src/CartographicTransformation.cpp b/particle-track-and-trace/src/CartographicTransformation.cpp index 5518d61..ee07e53 100644 --- a/particle-track-and-trace/src/CartographicTransformation.cpp +++ b/particle-track-and-trace/src/CartographicTransformation.cpp @@ -2,6 +2,9 @@ #include #include #include +#include +#include +#include vtkSmartPointer createNormalisedCamera() { vtkSmartPointer camera = vtkSmartPointer::New(); @@ -15,32 +18,93 @@ vtkSmartPointer createNormalisedCamera() { return camera; } -vtkSmartPointer getCartographicTransformMatrix(const std::shared_ptr 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), - 0, 2/(YMax-YMin), 0, -(YMax+YMin)/(YMax-YMin), - 0, 0, 1, 0, - 0, 0, 0, 1 - }; - - auto matrix = vtkSmartPointer::New(); - matrix->DeepCopy(eyeTransform); - return matrix; -} +//vtkSmartPointer getCartographicTransformMatrix(const std::shared_ptr 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), +// 0, 2/(YMax-YMin), 0, -(YMax+YMin)/(YMax-YMin), +// 0, 0, 1, 0, +// 0, 0, 0, 1 +// }; +// +// auto matrix = vtkSmartPointer::New(); +// matrix->DeepCopy(eyeTransform); +// return matrix; +//} // Assumes Normalised camera is used vtkSmartPointer createCartographicTransformFilter(const std::shared_ptr uvGrid) { - vtkNew transform; + auto proj = vtkSmartPointer::New(); + proj->SetName("merc"); - transform->SetMatrix(getCartographicTransformMatrix(uvGrid)); + auto geoTransform = vtkSmartPointer::New(); + geoTransform->SetDestinationProjection(proj); - vtkSmartPointer transformFilter = vtkSmartPointer::New(); - transformFilter->SetTransform(transform); + const double XMin = -15.875; + 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::New(); + scaleIntoNormalisedSpace->Scale(2/(width), 2/(height), 1); + scaleIntoNormalisedSpace->Translate(-(bottomLeft[0]+topRight[0])/2, -(bottomLeft[1] + topRight[1])/2, 0); + + auto totalProjection = vtkSmartPointer::New(); + totalProjection->Identity(); + totalProjection->Concatenate(scaleIntoNormalisedSpace); + totalProjection->Concatenate(geoTransform); + + vtkSmartPointer transformFilter = vtkSmartPointer::New(); + transformFilter->SetTransform(totalProjection); + + return transformFilter; +} + +vtkSmartPointer createInverseCartographicTransformFilter(const std::shared_ptr uvGrid) { + auto proj = vtkSmartPointer::New(); + proj->SetName("merc"); + + auto geoTransform = vtkSmartPointer::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::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::New(); + totalProjection->Identity(); + totalProjection->Concatenate(geoTransform); + totalProjection->Concatenate(scaleIntoNormalisedSpace); + + vtkSmartPointer transformFilter = vtkSmartPointer::New(); + transformFilter->SetTransform(totalProjection); + + return transformFilter; } diff --git a/particle-track-and-trace/src/CartographicTransformation.h b/particle-track-and-trace/src/CartographicTransformation.h index e59be8d..9b265e5 100644 --- a/particle-track-and-trace/src/CartographicTransformation.h +++ b/particle-track-and-trace/src/CartographicTransformation.h @@ -28,4 +28,6 @@ vtkSmartPointer getCartographicTransformMatrix(const std::shared_p * @return pointer to transform filter */ vtkSmartPointer createCartographicTransformFilter(const std::shared_ptr uvGrid); + +vtkSmartPointer createInverseCartographicTransformFilter(const std::shared_ptr uvGrid); #endif //NORMALISEDCARTOGRAPHICCAMERA_H diff --git a/particle-track-and-trace/src/commands/SpawnPointCallback.cpp b/particle-track-and-trace/src/commands/SpawnPointCallback.cpp index 3253e23..3327443 100644 --- a/particle-track-and-trace/src/commands/SpawnPointCallback.cpp +++ b/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 lon = " << worldPos[0] << " and lat = " << worldPos[1] << endl; + cout << "clicked on " << worldPos[1] << ", " << worldPos[0] << 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 &ren) { void SpawnPointCallback::setUVGrid(const std::shared_ptr &uvGrid) { this->uvGrid = uvGrid; - inverseCartographicProjection = getCartographicTransformMatrix(uvGrid); - inverseCartographicProjection->Invert(); + inverseCartographicProjection = createInverseCartographicTransformFilter(uvGrid)->GetTransform(); } diff --git a/particle-track-and-trace/src/commands/SpawnPointCallback.h b/particle-track-and-trace/src/commands/SpawnPointCallback.h index b7f7682..0896796 100644 --- a/particle-track-and-trace/src/commands/SpawnPointCallback.h +++ b/particle-track-and-trace/src/commands/SpawnPointCallback.h @@ -1,13 +1,13 @@ #ifndef SPAWNPOINTCALLBACK_H #define SPAWNPOINTCALLBACK_H - #include #include #include #include #include -#include +#include + #include "../advection/UVGrid.h" class SpawnPointCallback : public vtkCallbackCommand { @@ -30,7 +30,7 @@ private: vtkSmartPointer points; vtkSmartPointer ren; std::shared_ptr uvGrid; - vtkSmartPointer inverseCartographicProjection; + vtkSmartPointer inverseCartographicProjection; void Execute(vtkObject *caller, unsigned long evId, void *callData) override; diff --git a/particle-track-and-trace/src/layers/LGlyphLayer.cpp b/particle-track-and-trace/src/layers/LGlyphLayer.cpp index 3a4077a..60b12db 100644 --- a/particle-track-and-trace/src/layers/LGlyphLayer.cpp +++ b/particle-track-and-trace/src/layers/LGlyphLayer.cpp @@ -59,7 +59,7 @@ LGlyphLayer::LGlyphLayer(std::shared_ptr uvGrid, std::unique_ptr mapper; mapper->SetInputConnection(glyph2D->GetOutputPort()); mapper->Update(); - + vtkNew actor; actor->SetMapper(mapper); @@ -68,17 +68,11 @@ LGlyphLayer::LGlyphLayer(std::shared_ptr uvGrid, std::unique_ptrpoints->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 - - 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->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(); } @@ -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); } diff --git a/particle-track-and-trace/src/main.cpp b/particle-track-and-trace/src/main.cpp index f6a794e..c974467 100644 --- a/particle-track-and-trace/src/main.cpp +++ b/particle-track-and-trace/src/main.cpp @@ -17,6 +17,12 @@ #include "advection/kernel/RK4AdvectionKernel.h" #include "advection/kernel/SnapBoundaryConditionKernel.h" +#include +#include +#include +#include +#include + 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 = make_unique(DT); program->addLayer(new BackgroundImage("../../../../data/map_661-661.png")); @@ -37,6 +44,47 @@ int main() { program->render(); +// auto proj = vtkSmartPointer::New(); proj->SetName("merc"); auto geoTransform = vtkSmartPointer::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::New(); +// scaleIntoNormalisedSpace->Scale(2/(width), 2/(height), 1); +// scaleIntoNormalisedSpace->Translate(-(bottomLeft[0]+topRight[0])/2, -(bottomLeft[1] + topRight[1])/2, 0); +// +// auto totalProjection = vtkSmartPointer::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; }