diff --git a/particle-track-and-trace/src/CMakeLists.txt b/particle-track-and-trace/src/CMakeLists.txt index cf6696a..dabecca 100644 --- a/particle-track-and-trace/src/CMakeLists.txt +++ b/particle-track-and-trace/src/CMakeLists.txt @@ -48,8 +48,6 @@ add_executable(ParticleTrackTrace MACOSX_BUNDLE main.cpp layers/Layer.h layers/LGlyphLayer.cpp layers/LGlyphLayer.h - layers/LColLayer.cpp - layers/LColLayer.h Program.cpp Program.h commands/TimerCallbackCommand.h diff --git a/particle-track-and-trace/src/layers/EColLayer.cpp b/particle-track-and-trace/src/layers/EColLayer.cpp index e6c5325..c7df049 100644 --- a/particle-track-and-trace/src/layers/EColLayer.cpp +++ b/particle-track-and-trace/src/layers/EColLayer.cpp @@ -2,6 +2,7 @@ #include #include #include +#include #include #include #include @@ -24,8 +25,9 @@ #include "../CartographicTransformation.h" -using namespace std; +using std::numbers::pi; +using namespace std; EColLayer::EColLayer(std::shared_ptr uvGrid) { this->ren = vtkSmartPointer::New(); @@ -41,9 +43,53 @@ EColLayer::EColLayer(std::shared_ptr uvGrid) { this->strength->SetNumberOfComponents(1); this->strength->SetNumberOfTuples((numLats-1)*(numLons-1)); + this->direction = vtkSmartPointer::New(); + this->direction->SetName("direction"); + this->direction->SetNumberOfComponents(1); + this->direction->SetNumberOfTuples((numLats-1)*(numLons-1)); + readCoordinates(); } +/** + * Sets a given rgba colour to a range of values [start, end] in the lut. + * @param lut : lookuptable to operate on. + * @ param start : starting index of range to assign + * @ param end: ending index of range to assign + * @param r : red value [0,1] + * @param g : green value [0,1] + * @param n : blue value [0,1] + * @param a : alpha value [0,1] + */ +void setLutRange(vtkSmartPointer lut, int start, int end, double r, double g, double b, double a) { + for (int i=start; i <= end; i++) { + lut->SetTableValue(i, r, g, b, a); + } +} + +// builds a 4-way lookuptable, used to encode the directional component +vtkSmartPointer buildLutDirs() { + vtkNew lut; + lut->SetNumberOfColors(360); + lut->SetTableRange(0, 359); + lut->SetScaleToLinear(); + lut->Build(); + //currently builds a corkO cyclic colour map, divided into 8 colours (see https://www.fabiocrameri.ch/cycliccolourmaps/) + setLutRange(lut, 000, 020, 0.247, 0.243, 0.227, 1); + setLutRange(lut, 021, 060, 0.243, 0.267, 0.365, 1); + setLutRange(lut, 061, 100, 0.318, 0.416, 0.557, 1); + setLutRange(lut, 101, 140, 0.518, 0.620, 0.729, 1); + setLutRange(lut, 141, 180, 0.667, 0.757, 0.773, 1); + setLutRange(lut, 181, 220, 0.631, 0.769, 0.651, 1); + setLutRange(lut, 221, 260, 0.451, 0.639, 0.435, 1); + setLutRange(lut, 261, 300, 0.298, 0.431, 0.224, 1); + setLutRange(lut, 301, 340, 0.263, 0.310, 0.173, 1); + setLutRange(lut, 341, 359, 0.247, 0.243, 0.227, 1); + + lut->SetNanColor(0,0,0,0); + + return lut; +} // TODO: Bit of a superfunction here; can do with some refactoring. void EColLayer::readCoordinates() { @@ -64,27 +110,26 @@ void EColLayer::readCoordinates() { for (auto lat : uvGrid->lats) { double out[3] = {lon, lat, 0}; transform->TransformPoint(out, out); - // cout << "inserting point (" << pointId << "): " << out[0] << " " << out[1] << endl; points->InsertPoint(pointId++, out[0], out[1], 0); + // logic for adding cells if (latIndex > 0 and lonIndex > 0 ) { int idx = latIndex+lonIndex*numLats; - // cout << idx << " at " << lonIndex << " " << latIndex << endl; vtkNew l; l->SetNumberOfIds(4); l->SetId(0, idx-1); l->SetId(1, idx-numLats-1); l->SetId(2, idx-numLats); l->SetId(3, idx); - // cout << "inserting cell: " << idx-1 << " " << idx-numLats-1 << " " << idx- numLats<< " " << idx << endl; + double coords[12]; points->GetPoint(idx-1, coords); points->GetPoint(idx-numLats-1, coords+3); points->GetPoint(idx-numLats, coords+6); points->GetPoint(idx, coords+9); - // cout << "Inserting cell with points at: (" << coords[0] << " " << coords[1] << "), (" << coords[3] << " " << coords[4] << "), (" << coords[6] << " " << coords[7] << "), (" << coords[9] << " " << coords[10] << ")" << endl; data->InsertNextCell(VTK_QUAD, l); + // ltake the average of the four surrounding points as the cell's velocity. double u=0, v=0; for (int j=0; j < 2; j++ ) { for (int k=0; k < 2; k++ ) { @@ -95,7 +140,8 @@ void EColLayer::readCoordinates() { } u /= 4; v /= 4; - this->strength->SetTuple1(cellId++, std::sqrt(u*u + v*v)); + this->strength->SetTuple1(cellId, std::sqrt(u*u + v*v)); + this->direction->SetTuple1(cellId++, atan(u/v)*180/pi); } latIndex++; } @@ -103,11 +149,15 @@ void EColLayer::readCoordinates() { } data->GetCellData()->AddArray(this->strength); - data->GetCellData()->SetActiveScalars("strength"); + data->GetCellData()->AddArray(this->direction); + // data->GetCellData()->SetActiveScalars("strength"); vtkNew(mapper); mapper->SetInputData(data); + mapper->SetLookupTable(buildLutDirs()); + mapper->UseLookupTableScalarRangeOn(); mapper->Update(); + data->GetCellData()->SetActiveScalars("direction"); vtkNew actor; actor->SetMapper(mapper); @@ -137,7 +187,8 @@ void EColLayer::updateData(int t) { } u /= 4; v /= 4; - this->strength->SetTuple1(i++, std::sqrt(u*u + v*v)); + this->strength->SetTuple1(i, std::sqrt(u*u + v*v)); + this->direction->SetTuple1(i++, atan(u/v)*180/pi); } } this->strength->Modified(); diff --git a/particle-track-and-trace/src/layers/EColLayer.h b/particle-track-and-trace/src/layers/EColLayer.h index 4934728..1097d21 100644 --- a/particle-track-and-trace/src/layers/EColLayer.h +++ b/particle-track-and-trace/src/layers/EColLayer.h @@ -15,6 +15,7 @@ class EColLayer : public Layer { private: vtkSmartPointer strength; + vtkSmartPointer direction; std::shared_ptr uvGrid; int numLats; int numLons; diff --git a/particle-track-and-trace/src/main.cpp b/particle-track-and-trace/src/main.cpp index 00cc38d..7d768d5 100644 --- a/particle-track-and-trace/src/main.cpp +++ b/particle-track-and-trace/src/main.cpp @@ -35,7 +35,7 @@ int main() { l->setDt(DT); unique_ptr program = make_unique(DT); - program->addLayer(new BackgroundImage(dataPath + "/map_661-661.png")); + program->addLayer(new BackgroundImage(dataPath + "/map_qgis_1035.png")); // program->addLayer(new EGlyphLayer(uvGrid)); program->addLayer(new EColLayer(uvGrid)); program->addLayer(l);