diff --git a/src/Cxx/Visualization/CurvatureBandsWithGlyphs.cxx b/src/Cxx/Visualization/CurvatureBandsWithGlyphs.cxx
index c4d18fe0d4913ffae5ff2bd3a8c27155d6b9140a..6f29ad65212252820d2006bd8d95dd9cca3383fa 100644
--- a/src/Cxx/Visualization/CurvatureBandsWithGlyphs.cxx
+++ b/src/Cxx/Visualization/CurvatureBandsWithGlyphs.cxx
@@ -2,7 +2,6 @@
#include <vtkArrowSource.h>
#include <vtkBandedPolyDataContourFilter.h>
#include <vtkCamera.h>
-#include <vtkCameraOrientationWidget.h>
#include <vtkCleanPolyData.h>
#include <vtkColorSeries.h>
#include <vtkColorTransferFunction.h>
@@ -45,6 +44,11 @@
#include <vtkVariantArray.h>
#include <vtkVersion.h>
+#if VTK_VERSION_NUMBER >= 90020210809ULL
+#define HAS_COW
+#include <vtkCameraOrientationWidget.h>
+#endif
+
#include <algorithm>
#include <array>
#include <cctype>
@@ -191,20 +195,25 @@ void PrintBandsFrequencies(std::map<int, std::vector<double>> const& bands,
int main(int, char*[])
{
// Get the surface
- //std::string desiredSurface = "Hills";
- //std::string desiredSurface = "ParametricTorus";
- //std::string desiredSurface = "Plane";
- std::string desiredSurface = "RandomHills";
- //std::string desiredSurface = "Sphere";
- //std::string desiredSurface = "Torus";
+ std::string desiredSurface{"RandomHills"};
+ // desiredSurface = "Hills";
+ // desiredSurface = "ParametricTorus";
+ // desiredSurface = "Plane";
+ // desiredSurface = "RandomHills";
+ // desiredSurface = "Sphere";
+ // desiredSurface = "Torus";
auto source = GetSource(desiredSurface);
// The length of the normal arrow glyphs.
auto scaleFactor = 1.0;
if (desiredSurface == "Hills")
- scaleFactor = 0.5;
+ {
+ scaleFactor = 0.5;
+ }
if (desiredSurface == "Sphere")
+ {
scaleFactor = 2.0;
+ }
std::cout << desiredSurface << std::endl;
auto gaussianCurvature = true;
@@ -270,19 +279,21 @@ int main(int, char*[])
lut->SetTableRange(scalarRangeCurvatures);
lut1->SetTableRange(scalarRangeElevation);
auto numberOfBands = lut->GetNumberOfTableValues();
- auto bands = GetBands(scalarRangeCurvatures, numberOfBands, 10);
+ auto precision = 10;
+ auto bands = GetBands(scalarRangeCurvatures, numberOfBands, precision, false);
+
if (desiredSurface == "RandomHills")
{
// These are my custom bands.
// Generated by first running:
- // bands = GetBands(scalarRangeCurvatures, numberOfBands, false);
+ // bands = GetBands(scalarRangeCurvatures, numberOfBands, precision, false);
// then:
// std::vector<int> freq = Frequencies(bands, src);
// PrintBandsFrequencies(bands, freq);
// Finally using the output to create this table:
// std::vector<std::array<double, 2>> myBands = {
// {-0.630, -0.190}, {-0.190, -0.043}, {-0.043, -0.0136},
- // {-0.0136, 0.0158}, {0.0158, 0.0452}, {0.0452, 0.0746},
+ // {-0.0136, 0.0158}, {0.0158, 0.0452}, {0.0452, 0.0746},12
// {0.0746, 0.104}, {0.104, 0.251}, {0.251, 1.131}};
// This demonstrates that the gaussian curvature of the surface
// is mostly planar with some hyperbolic regions (saddle points)
@@ -466,10 +477,12 @@ int main(int, char*[])
renWin->AddRenderer(ren);
// Important: The interactor must be set prior to enabling the widget.
iren->SetRenderWindow(renWin);
+#ifdef HAS_COW
vtkNew<vtkCameraOrientationWidget> camOrientManipulator;
camOrientManipulator->SetParentRenderer(ren);
// Enable the widget.
camOrientManipulator->On();
+#endif
// add actors
ren->AddViewProp(srcActor);
diff --git a/src/Cxx/Visualization/ElevationBandsWithGlyphs.cxx b/src/Cxx/Visualization/ElevationBandsWithGlyphs.cxx
index 82d040ea21ce1962651f260288ed9bd5f3da0624..0e465eb893588e12859cee750edafb483e0a79dd 100644
--- a/src/Cxx/Visualization/ElevationBandsWithGlyphs.cxx
+++ b/src/Cxx/Visualization/ElevationBandsWithGlyphs.cxx
@@ -2,11 +2,13 @@
#include <vtkArrowSource.h>
#include <vtkBandedPolyDataContourFilter.h>
#include <vtkCamera.h>
-#include <vtkCameraOrientationWidget.h>
#include <vtkCleanPolyData.h>
#include <vtkColorSeries.h>
#include <vtkColorTransferFunction.h>
+#include <vtkDelaunay2D.h>
+#include <vtkDoubleArray.h>
#include <vtkElevationFilter.h>
+#include <vtkFloatArray.h>
#include <vtkGlyph3D.h>
#include <vtkInteractorStyleTrackballCamera.h>
#include <vtkLookupTable.h>
@@ -20,19 +22,28 @@
#include <vtkPointData.h>
#include <vtkPolyData.h>
#include <vtkPolyDataMapper.h>
+#include <vtkPolyDataNormals.h>
#include <vtkProperty.h>
#include <vtkRenderWindow.h>
#include <vtkRenderWindowInteractor.h>
#include <vtkRenderer.h>
#include <vtkReverseSense.h>
#include <vtkScalarBarActor.h>
+#include <vtkSmartPointer.h>
#include <vtkSphereSource.h>
#include <vtkSuperquadricSource.h>
#include <vtkTextProperty.h>
#include <vtkTransform.h>
#include <vtkTransformPolyDataFilter.h>
#include <vtkTriangleFilter.h>
+#include <vtkVariant.h>
#include <vtkVariantArray.h>
+#include <vtkVersion.h>
+
+#if VTK_VERSION_NUMBER >= 90020210809ULL
+#define HAS_COW
+#include <vtkCameraOrientationWidget.h>
+#endif
#include <algorithm>
#include <array>
@@ -44,23 +55,52 @@
#include <iomanip>
#include <iostream>
#include <iterator>
+#include <numeric>
+#include <set>
#include <sstream>
#include <string>
#include <vector>
namespace {
+//! Generate elevations over the surface.
+/*!
+@param src - the vtkPolyData source.
+@return elev - the elevations.
+*/
+vtkSmartPointer<vtkPolyData> GetElevations(vtkPolyData* src);
-//! Divide a range into bands
+vtkSmartPointer<vtkPolyData> GetHills();
+vtkSmartPointer<vtkPolyData> GetParametricHills();
+vtkSmartPointer<vtkPolyData> GetParametricTorus();
+vtkSmartPointer<vtkPolyData> GetPlane();
+vtkSmartPointer<vtkPolyData> GetSphere();
+vtkSmartPointer<vtkPolyData> GetTorus();
+vtkSmartPointer<vtkPolyData> GetSource(std::string const& source);
+
+vtkSmartPointer<vtkColorSeries> GetColorSeries();
+
+vtkSmartPointer<vtkLookupTable> GetCategoricalLUT();
+
+vtkSmartPointer<vtkLookupTable> GetOrdinaLUT();
+
+vtkSmartPointer<vtkLookupTable> GetDivergingLUT();
+
+vtkSmartPointer<vtkLookupTable> ReverseLUT(vtkLookupTable* lut);
+
+//! Glyph the normals on the surface.
/*!
-@param dR - [min, max] the range that is to be covered by the bands.
-@param numberOfBands - the number of bands, a positive integer.
-@param nearestInteger - if True then [floor(min), ceil(max)] is used.
-@return A map consisting of the band inxex and [min, midpoint, max] for each
-band.
+@param src - the vtkPolyData source.
+#param scaleFactor - the scale factor for the glyphs.
+@param reverseNormals - if True the normals on the surface are reversed.
+@return The glyphs.
*/
-std::map<int, std::vector<double>> MakeBands(double const dR[2],
- int const& numberOfBands,
- bool const& nearestInteger);
+vtkNew<vtkGlyph3D> GetGlyphs(vtkPolyData* src, double const& scaleFactor = 1.0,
+ bool const& reverseNormals = false);
+
+std::map<int, std::vector<double>> GetBands(double const dR[2],
+ int const& numberOfBands,
+ int const& precision = 2,
+ bool const& nearestInteger = false);
//! Divide a range into custom bands
/*!
@@ -76,8 +116,8 @@ like this: [[r1, r2], [r2, r3], [r3, r4]...]
band.
*/
std::map<int, std::vector<double>>
-MakeCustomBands(double const dR[2], int const& numberOfBands,
- std::vector<std::array<double, 2>> const& myBands);
+GetCustomBands(double const dR[2], int const& numberOfBands,
+ std::vector<std::array<double, 2>> const& myBands);
//! Divide a range into integral bands
/*!
@@ -86,155 +126,76 @@ Divide a range into bands
@returnA map consisting of the band inxex and [min, midpoint, max] for each
band.
*/
-std::map<int, std::vector<double>> MakeIntegralBands(double const dR[2]);
-
-//! Generate elevations over the surface.
-/*!
-@param src - the vtkPolyData source.
-@return elev - the elevations.
-*/
-vtkSmartPointer<vtkPolyData> MakeElevations(vtkPolyData* src);
-
-//! Make a parametric hills surface as the source.
-/*!
-@return - vtkPolyData.
-*/
-vtkSmartPointer<vtkPolyData> MakeParametricHills();
-
-//! Make a parametric torus as the source.
-/*!
-@return - vtkPolyData.
-*/
-vtkSmartPointer<vtkPolyData> MakeParametricTorus();
-
-//! Make a Plane as the source.
-/*!
-@return - vtkPolyData.
-*/
-vtkSmartPointer<vtkPolyData> MakePlane();
-
-//! Make a MakeSphere as the source.
-/*!
-@return - vtkPolyData.
-*/
-vtkSmartPointer<vtkPolyData> MakeSphere();
-
-//! Make a torus as the source.
-/*!
-@return - vtkPolyData.
-*/
-vtkSmartPointer<vtkPolyData> MakeTorus();
-
-vtkSmartPointer<vtkColorSeries> GetColorSeries();
-
-vtkSmartPointer<vtkLookupTable> MakeCategoricalLUT();
-
-vtkSmartPointer<vtkLookupTable> MakeOrdinaLUT();
-
-vtkSmartPointer<vtkLookupTable> MakeDivergingLUT();
-
-vtkSmartPointer<vtkLookupTable> ReverseLUT(vtkLookupTable* lut);
+std::map<int, std::vector<double>> GetIntegralBands(double const dR[2]);
//! Count the number of scalars in each band.
-/*!
-@param bands - the bands.
-@param src - the vtkPolyData source.
-@return The frequencies of the scalars in each band.
-*/
-std::map<int, int> Frequencies(std::map<int, std::vector<double>>& bands,
- vtkPolyData* src);
-
-std::pair<int, int> AdjustFrequencyRanges(std::map<int, int>& freq);
-
-void PrintBands(std::map<int, std::vector<double>> const& bands);
-
-void PrintFrequencies(std::map<int, int> const& freq);
-
-void PrintBandsFrequencies(std::map<int, std::vector<double>> const& bands,
+/*
+ * The scalars used are the active scalars in the polydata.
+ *
+ * @param bands - the bands.
+ * @param src - the vtkPolyData source.
+ * @return The frequencies of the scalars in each band.
+ */
+std::map<int, int> GetFrequencies(std::map<int, std::vector<double>>& bands,
+ vtkPolyData* src);
+//!
+/*
+ * The bands and frequencies are adjusted so that the first and last
+ * frequencies in the range are non-zero.
+ * @param bands: The bands.
+ * @param freq: The frequencies.
+ */
+void AdjustFrequencyRanges(std::map<int, std::vector<double>>& bands,
std::map<int, int>& freq);
-//! Glyph the normals on the surface.
-/*!
-@param src - the vtkPolyData source.
-@param reverseNormals - if True the normals on the surface are reversed.
-@return The glyphs.
-*/
-vtkNew<vtkGlyph3D> MakeGlyphs(vtkPolyData* src, bool const& reverseNormals);
+void PrintBandsFrequencies(std::map<int, std::vector<double>> const& bands,
+ std::map<int, int>& freq, int const& precision = 2);
} // namespace
-//! Make and display the surface.
int main(int, char*[])
{
// Get the surface
- // std::string desiredSurface = "ParametricTorus";
- // std::string desiredSurface = "Plane";
- std::string desiredSurface = "RandomHills";
- // std::string desiredSurface = "Sphere";
- // std::string desiredSurface = "Torus";
- auto lcSurface = desiredSurface;
- std::transform(lcSurface.begin(), lcSurface.end(), lcSurface.begin(),
- [](char c) { return std::tolower(c); });
- std::map<std::string, int> availableSurfaces = {{"parametrictorus", 0},
- {"plane", 1},
- {"randomhills", 2},
- {"sphere", 3},
- {"torus", 4}};
- vtkSmartPointer<vtkPolyData> src;
- if (availableSurfaces.find(lcSurface) == availableSurfaces.end())
- {
- std::cout << "No surface specified." << std::endl;
- return EXIT_FAILURE;
- }
- switch (availableSurfaces[lcSurface])
- {
- case 0: {
- src = MakeParametricTorus();
- break;
- }
- case 1: {
- src = MakePlane();
- src = MakeElevations(src);
- break;
- }
- case 2: {
- src = MakeParametricHills();
- break;
- }
- case 3: {
- src = MakeSphere();
- src = MakeElevations(src);
- break;
- }
- case 4: {
- src = MakeTorus();
- src = MakeElevations(src);
- break;
- }
- default: {
- std::cout << "No surface specified." << std::endl;
- return EXIT_FAILURE;
+ std::string desiredSurface{"RandomHills"};
+ // desiredSurface = "Hills";
+ // desiredSurface = "ParametricTorus";
+ // desiredSurface = "Plane";
+ // desiredSurface = "RandomHills";
+ // desiredSurface = "Sphere";
+ // desiredSurface = "Torus";
+ auto source = GetSource(desiredSurface);
+
+ // The length of the normal arrow glyphs.
+ auto scaleFactor = 1.0;
+ if (desiredSurface == "Hills")
+ {
+ scaleFactor = 0.5;
}
+ if (desiredSurface == "Sphere")
+ {
+ scaleFactor = 2.0;
}
std::cout << desiredSurface << std::endl;
- src->GetPointData()->SetActiveScalars("Elevation");
- auto scalarRange = src->GetPointData()->GetScalars("Elevation")->GetRange();
+ source->GetPointData()->SetActiveScalars("Elevation");
+ auto scalarRange =
+ source->GetPointData()->GetScalars("Elevation")->GetRange();
- auto lut = MakeCategoricalLUT();
- auto lut1 = MakeOrdinaLUT();
+ auto lut = GetCategoricalLUT();
+ auto lut1 = GetOrdinaLUT();
lut->SetTableRange(scalarRange);
lut1->SetTableRange(scalarRange);
vtkIdType numberOfBands = lut->GetNumberOfTableValues();
- std::map<int, std::vector<double>> bands;
+ auto precision = 10;
+ auto bands = GetBands(scalarRange, numberOfBands, precision, false);
if (desiredSurface == "RandomHills")
{
// These are my custom bands.
// Generated by first running:
- // bands = MakeBands(scalarRange, numberOfBands, false);
+ // bands = GetBands(scalarRange, numberOfBands, precision, false);
// then:
- // std::vector<int> freq = Frequencies(bands, src);
+ // std::vector<int> freq = Frequencies(bands, source);
// PrintBandsFrequencies(bands, freq);
// Finally using the output to create this table:
std::vector<std::array<double, 2>> myBands = {
@@ -242,72 +203,27 @@ int main(int, char*[])
{4.0, 5.0}, {5.0, 6.0}, {6.0, 7.0}, {7.0, 8.0}};
// Comment this out if you want to see how allocating
// equally spaced bands works.
- bands = MakeCustomBands(scalarRange, numberOfBands, myBands);
- // bands = MakeBands(scalarRange, numberOfBands, false);
+ bands = GetCustomBands(scalarRange, numberOfBands, myBands);
+ // bands = GetBands(scalarRange, numberOfBands, precision, false);
// Adjust the number of table values
- lut->SetNumberOfTableValues(static_cast<vtkIdType>(bands.size()));
- lut1->SetNumberOfTableValues(static_cast<vtkIdType>(bands.size()));
- }
- else
- {
- bands = MakeBands(scalarRange, numberOfBands, false);
}
+ lut->SetNumberOfTableValues(static_cast<vtkIdType>(bands.size()));
+ lut1->SetNumberOfTableValues(static_cast<vtkIdType>(bands.size()));
// Let's do a frequency table.
- // The number of scalars in each band.
- std::map<int, int> freq = Frequencies(bands, src);
-
- std::vector<int> keys;
- for (int k = 0; k < freq.size(); ++k)
- {
- keys.push_back(k);
- }
- auto minKey = std::min_element(keys.begin(), keys.end());
- auto maxKey = std::max_element(keys.begin(), keys.end());
- if (lcSurface == "sphere")
- {
- freq[0] = 0;
- }
- auto freqRange = AdjustFrequencyRanges(freq);
- std::map<int, int>::iterator freqItr;
- freqItr = freq.find(freqRange.first);
- freq.erase(freq.begin(), freqItr);
- freqItr = ++freq.find(freqRange.second);
- freq.erase(freqItr, freq.end());
- std::map<int, std::vector<double>>::iterator bandItr;
- bandItr = bands.find(freqRange.first);
- bands.erase(bands.begin(), bandItr);
- bandItr = ++bands.find(freqRange.second);
- bands.erase(bandItr, bands.end());
- // Reindex freq and bands.
- std::map<int, int> adjFreq;
- int idx = 0;
- for (auto p : freq)
- {
- adjFreq[idx] = p.second;
- ++idx;
- }
- std::map<int, std::vector<double>> adjBands;
- idx = 0;
- for (auto p : bands)
- {
- adjBands[idx] = p.second;
- ++idx;
- }
- // PrintBandsFrequencies(bands, freq);
- PrintBandsFrequencies(adjBands, adjFreq);
- freq.erase(freq.begin(), freq.end());
- bands.erase(bands.begin(), bands.end());
+ auto freq = GetFrequencies(bands, source);
+ AdjustFrequencyRanges(bands, freq);
+ PrintBandsFrequencies(bands, freq);
- scalarRange[0] = adjBands.begin()->second[0];
- scalarRange[1] = std::prev(adjBands.end())->second[2];
+ scalarRange[0] = bands.begin()->second[0];
+ scalarRange[1] = std::prev(bands.end())->second[2];
lut->SetTableRange(scalarRange);
- lut->SetNumberOfTableValues(adjBands.size());
+ lut->SetNumberOfTableValues(bands.size());
// We will use the midpoint of the band as the label.
std::vector<std::string> labels;
- for (std::map<int, std::vector<double>>::const_iterator p = adjBands.begin();
- p != adjBands.end(); ++p)
+ for (std::map<int, std::vector<double>>::const_iterator p = bands.begin();
+ p != bands.end(); ++p)
{
std::ostringstream os;
os << std::fixed << std::setw(6) << std::setprecision(2) << p->second[1];
@@ -330,23 +246,21 @@ int main(int, char*[])
// Create the contour bands.
vtkNew<vtkBandedPolyDataContourFilter> bcf;
- bcf->SetInputData(src);
+ bcf->SetInputData(source);
// Use either the minimum or maximum value for each band.
int i = 0;
- for (std::map<int, std::vector<double>>::const_iterator p = adjBands.begin();
- p != adjBands.end(); ++p)
+ for (std::map<int, std::vector<double>>::const_iterator p = bands.begin();
+ p != bands.end(); ++p)
{
bcf->SetValue(i, p->second[2]);
++i;
}
-
// We will use an indexed lookup table.
bcf->SetScalarModeToIndex();
bcf->GenerateContourEdgesOn();
// Generate the glyphs on the original surface.
-
- auto glyph = MakeGlyphs(src, false);
+ auto glyph = GetGlyphs(source, scaleFactor, false);
// ------------------------------------------------------------
// Create the mappers and actors
@@ -423,10 +337,12 @@ int main(int, char*[])
renWin->AddRenderer(ren);
// Important: The interactor must be set prior to enabling the widget.
iren->SetRenderWindow(renWin);
+#ifdef HAS_COW
vtkNew<vtkCameraOrientationWidget> camOrientManipulator;
camOrientManipulator->SetParentRenderer(ren);
// Enable the widget.
camOrientManipulator->On();
+#endif
// add actors
ren->AddViewProp(srcActor);
@@ -455,126 +371,123 @@ int main(int, char*[])
}
namespace {
-
-std::map<int, std::vector<double>> MakeBands(double const dR[2],
- int const& numberOfBands,
- bool const& nearestInteger)
+vtkSmartPointer<vtkPolyData> GetElevations(vtkPolyData* src)
{
- std::map<int, std::vector<double>> bands;
- if ((dR[1] < dR[0]) || (numberOfBands <= 0))
- {
- return bands;
- }
- double x[2];
- for (int i = 0; i < 2; ++i)
- {
- x[i] = dR[i];
- }
- if (nearestInteger)
- {
- x[0] = std::floor(x[0]);
- x[1] = std::ceil(x[1]);
- }
- double dx = (x[1] - x[0]) / static_cast<double>(numberOfBands);
- std::vector<double> b;
- b.push_back(x[0]);
- b.push_back(x[0] + dx / 2.0);
- b.push_back(x[0] + dx);
- for (int i = 0; i < numberOfBands; ++i)
+ double bounds[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
+ src->GetBounds(bounds);
+ if (std::abs(bounds[2]) < 1.0e-8 && std::abs(bounds[3]) < 1.0e-8)
{
- bands[i] = b;
- for (std::vector<double>::iterator p = b.begin(); p != b.end(); ++p)
- {
- *p = *p + dx;
- }
+ bounds[3] = bounds[2] + 1;
}
- return bands;
+ vtkNew<vtkElevationFilter> elevFilter;
+ elevFilter->SetInputData(src);
+ elevFilter->SetLowPoint(0, bounds[2], 0);
+ elevFilter->SetHighPoint(0, bounds[3], 0);
+ elevFilter->SetScalarRange(bounds[2], bounds[3]);
+ elevFilter->Update();
+
+ return elevFilter->GetPolyDataOutput();
}
-std::map<int, std::vector<double>>
-MakeCustomBands(double const dR[2], int const& numberOfBands,
- std::vector<std::array<double, 2>> const& myBands)
+vtkSmartPointer<vtkPolyData> GetHills()
{
- std::map<int, std::vector<double>> bands;
- if ((dR[1] < dR[0]) || (numberOfBands <= 0))
+ // Create four hills on a plane.
+ // This will have regions of negative, zero and positive Gsaussian curvatures.
+
+ auto xRes = 50;
+ auto yRes = 50;
+ auto xMin = -5.0;
+ auto xMax = 5.0;
+ auto dx = (xMax - xMin) / (xRes - 1.0);
+ auto yMin = -5.0;
+ auto yMax = 5.0;
+ auto dy = (yMax - yMin) / (xRes - 1.0);
+
+ // Make a grid.
+ vtkNew<vtkPoints> points;
+ for (auto i = 0; i < xRes; ++i)
{
- return bands;
+ auto x = xMin + i * dx;
+ for (auto j = 0; j < yRes; ++j)
+ {
+ auto y = yMin + j * dy;
+ points->InsertNextPoint(x, y, 0);
+ }
}
- std::vector<std::array<double, 2>> x;
- std::copy(myBands.begin(), myBands.end(), std::back_inserter(x));
+ // Add the grid points to a polydata object.
+ vtkNew<vtkPolyData> plane;
+ plane->SetPoints(points);
- // Determine the index of the range minimum and range maximum.
- int idxMin = 0;
- for (auto idx = 0; idx < myBands.size(); ++idx)
+ // Triangulate the grid.
+ vtkNew<vtkDelaunay2D> delaunay;
+ delaunay->SetInputData(plane);
+ delaunay->Update();
+
+ auto polydata = delaunay->GetOutput();
+
+ vtkNew<vtkDoubleArray> elevation;
+ elevation->SetNumberOfTuples(points->GetNumberOfPoints());
+
+ // We define the parameters for the hills here.
+ // [[0: x0, 1: y0, 2: x variance, 3: y variance, 4: amplitude]...]
+ std::vector<std::array<double, 5>> hd{{-2.5, -2.5, 2.5, 6.5, 3.5},
+ {2.5, 2.5, 2.5, 2.5, 2},
+ {5.0, -2.5, 1.5, 1.5, 2.5},
+ {-5.0, 5, 2.5, 3.0, 3}};
+ std::array<double, 2> xx{0.0, 0.0};
+ for (auto i = 0; i < points->GetNumberOfPoints(); ++i)
{
- if (dR[0] < myBands[idx][1] && dR[0] >= myBands[idx][0])
+ auto x = polydata->GetPoint(i);
+ for (size_t j = 0; j < hd.size(); ++j)
{
- idxMin = idx;
- break;
+ xx[0] = std::pow(x[0] - hd[j][0] / hd[j][2], 2.0);
+ xx[1] = std::pow(x[1] - hd[j][1] / hd[j][3], 2.0);
+ x[2] += hd[j][4] * std::exp(-(xx[0] + xx[1]) / 2.0);
}
+ polydata->GetPoints()->SetPoint(i, x);
+ elevation->SetValue(i, x[2]);
}
- int idxMax = static_cast<int>(myBands.size()) - 1;
- for (auto idx = myBands.size() - 1; idx >= 0; --idx)
+
+ vtkNew<vtkFloatArray> textures;
+ textures->SetNumberOfComponents(2);
+ textures->SetNumberOfTuples(2 * polydata->GetNumberOfPoints());
+ textures->SetName("Textures");
+
+ for (auto i = 0; i < xRes; ++i)
{
- if (dR[1] < myBands[idx][1] && dR[1] >= myBands[idx][0])
+ float tc[2];
+ tc[0] = i / (xRes - 1.0);
+ for (auto j = 0; j < yRes; ++j)
{
- idxMax = static_cast<int>(idx);
- break;
+ // tc[1] = 1.0 - j / (yRes - 1.0);
+ tc[1] = j / (yRes - 1.0);
+ textures->SetTuple(static_cast<vtkIdType>(i) * yRes + j, tc);
}
}
- // Set the minimum to match the range minimum.
- x[idxMin][0] = dR[0];
- x[idxMax][1] = dR[1];
- for (int i = idxMin; i < idxMax + 1; ++i)
- {
- std::vector<double> b(3);
- b[0] = x[i][0];
- b[1] = x[i][0] + (x[i][1] - x[i][0]) / 2.0;
- b[2] = x[i][1];
- bands[i] = b;
- }
- return bands;
-}
+ polydata->GetPointData()->SetScalars(elevation);
+ polydata->GetPointData()->GetScalars()->SetName("Elevation");
+ polydata->GetPointData()->SetTCoords(textures);
-std::map<int, std::vector<double>> MakeIntegralBands(double const dR[2])
-{
- std::map<int, std::vector<double>> bands;
- if (dR[1] < dR[0])
- {
- return bands;
- }
- double x[2];
- for (int i = 0; i < 2; ++i)
- {
- x[i] = dR[i];
- }
- x[0] = std::floor(x[0]);
- x[1] = std::ceil(x[1]);
- int numberOfBands = static_cast<int>(std::abs(x[1]) + std::abs(x[0]));
- return MakeBands(x, numberOfBands, false);
-}
+ vtkNew<vtkPolyDataNormals> normals;
+ normals->SetInputData(polydata);
+ normals->SetInputData(polydata);
+ normals->SetFeatureAngle(30);
+ normals->SplittingOff();
-vtkSmartPointer<vtkPolyData> MakeElevations(vtkPolyData* src)
-{
- double bounds[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
- src->GetBounds(bounds);
- if (std::abs(bounds[2]) < 1.0e-8 && std::abs(bounds[3]) < 1.0e-8)
- {
- bounds[3] = bounds[2] + 1;
- }
- vtkNew<vtkElevationFilter> elevFilter;
- elevFilter->SetInputData(src);
- elevFilter->SetLowPoint(0, bounds[2], 0);
- elevFilter->SetHighPoint(0, bounds[3], 0);
- elevFilter->SetScalarRange(bounds[2], bounds[3]);
- elevFilter->Update();
+ vtkNew<vtkTransform> tr1;
+ tr1->RotateX(-90);
- return elevFilter->GetPolyDataOutput();
+ vtkNew<vtkTransformPolyDataFilter> tf1;
+ tf1->SetInputConnection(normals->GetOutputPort());
+ tf1->SetTransform(tr1);
+ tf1->Update();
+
+ return tf1->GetOutput();
}
-vtkSmartPointer<vtkPolyData> MakeParametricHills()
+vtkSmartPointer<vtkPolyData> GetParametricHills()
{
vtkNew<vtkParametricRandomHills> fn;
fn->AllowRandomGenerationOn();
@@ -606,7 +519,7 @@ vtkSmartPointer<vtkPolyData> MakeParametricHills()
return transformFilter->GetOutput();
}
-vtkSmartPointer<vtkPolyData> MakeParametricTorus()
+vtkSmartPointer<vtkPolyData> GetParametricTorus()
{
vtkNew<vtkParametricTorus> fn;
fn->SetRingRadius(5);
@@ -636,7 +549,7 @@ vtkSmartPointer<vtkPolyData> MakeParametricTorus()
return transformFilter->GetOutput();
}
-vtkSmartPointer<vtkPolyData> MakePlane()
+vtkSmartPointer<vtkPolyData> GetPlane()
{
vtkNew<vtkPlaneSource> source;
source->SetOrigin(-10.0, -10.0, 0.0);
@@ -670,7 +583,7 @@ vtkSmartPointer<vtkPolyData> MakePlane()
return cleaner->GetOutput();
}
-vtkSmartPointer<vtkPolyData> MakeSphere()
+vtkSmartPointer<vtkPolyData> GetSphere()
{
vtkNew<vtkSphereSource> source;
source->SetCenter(0.0, 0.0, 0.0);
@@ -682,7 +595,7 @@ vtkSmartPointer<vtkPolyData> MakeSphere()
return source->GetOutput();
}
-vtkSmartPointer<vtkPolyData> MakeTorus()
+vtkSmartPointer<vtkPolyData> GetTorus()
{
vtkNew<vtkSuperquadricSource> source;
source->SetCenter(0.0, 0.0, 0.0);
@@ -711,6 +624,44 @@ vtkSmartPointer<vtkPolyData> MakeTorus()
return cleaner->GetOutput();
}
+vtkSmartPointer<vtkPolyData> GetSource(std::string const& source)
+{
+ std::string surface = source;
+ std::transform(surface.begin(), surface.end(), surface.begin(),
+ [](unsigned char c) { return std::tolower(c); });
+ std::map<std::string, int> available_surfaces = {
+ {"hills", 0}, {"parametrictorus", 1}, {"plane", 2},
+ {"randomhills", 3}, {"sphere", 4}, {"torus", 5}};
+ if (available_surfaces.find(surface) == available_surfaces.end())
+ {
+ std::cout << "The surface is not available." << std::endl;
+ std::cout << "Using RandomHills instead." << std::endl;
+ surface = "randomhills";
+ }
+ switch (available_surfaces[surface])
+ {
+ case 0:
+ return GetHills();
+ break;
+ case 1:
+ return GetParametricTorus();
+ break;
+ case 2:
+ return GetElevations(GetPlane());
+ break;
+ case 3:
+ return GetParametricHills();
+ break;
+ case 4:
+ return GetElevations(GetSphere());
+ break;
+ case 5:
+ return GetElevations(GetTorus());
+ break;
+ }
+ return GetParametricHills();
+}
+
vtkSmartPointer<vtkColorSeries> GetColorSeries()
{
@@ -729,7 +680,7 @@ vtkSmartPointer<vtkColorSeries> GetColorSeries()
return colorSeries;
}
-vtkSmartPointer<vtkLookupTable> MakeCategoricalLUT()
+vtkSmartPointer<vtkLookupTable> GetCategoricalLUT()
{
vtkSmartPointer<vtkColorSeries> colorSeries = GetColorSeries();
// Make the lookup table.
@@ -740,7 +691,7 @@ vtkSmartPointer<vtkLookupTable> MakeCategoricalLUT()
return lut;
}
-vtkSmartPointer<vtkLookupTable> MakeOrdinaLUT()
+vtkSmartPointer<vtkLookupTable> GetOrdinaLUT()
{
vtkSmartPointer<vtkColorSeries> colorSeries = GetColorSeries();
// Make the lookup table.
@@ -764,7 +715,7 @@ vtkSmartPointer<vtkLookupTable> MakeOrdinaLUT()
*
*/
// clang-format on
-vtkSmartPointer<vtkLookupTable> MakeDivergingLUT()
+vtkSmartPointer<vtkLookupTable> GetDivergingLUT()
{
vtkNew<vtkColorTransferFunction> ctf;
@@ -816,11 +767,169 @@ vtkSmartPointer<vtkLookupTable> ReverseLUT(vtkLookupTable* lut)
return lutr;
}
-std::map<int, int> Frequencies(std::map<int, std::vector<double>>& bands,
- vtkPolyData* src)
+vtkNew<vtkGlyph3D> GetGlyphs(vtkPolyData* src, double const& scaleFactor,
+ bool const& reverseNormals)
+{
+ // Sometimes the contouring algorithm can create a volume whose gradient
+ // vector and ordering of polygon(using the right hand rule) are
+ // inconsistent. vtkReverseSense cures this problem.
+ vtkNew<vtkReverseSense> reverse;
+ vtkNew<vtkMaskPoints> maskPts;
+ maskPts->SetOnRatio(5);
+ maskPts->RandomModeOn();
+ if (reverseNormals)
+ {
+ reverse->SetInputData(src);
+ reverse->ReverseCellsOn();
+ reverse->ReverseNormalsOn();
+ maskPts->SetInputConnection(reverse->GetOutputPort());
+ }
+ else
+ {
+ maskPts->SetInputData(src);
+ }
+
+ // Source for the glyph filter
+ vtkNew<vtkArrowSource> arrow;
+ arrow->SetTipResolution(16);
+ arrow->SetTipLength(0.3);
+ arrow->SetTipRadius(0.1);
+
+ vtkNew<vtkGlyph3D> glyph;
+ glyph->SetSourceConnection(arrow->GetOutputPort());
+ glyph->SetInputConnection(maskPts->GetOutputPort());
+ glyph->SetVectorModeToUseNormal();
+ glyph->SetScaleFactor(scaleFactor);
+ glyph->SetColorModeToColorByVector();
+ glyph->SetScaleModeToScaleByVector();
+ glyph->OrientOn();
+ glyph->Update();
+
+ return glyph;
+}
+
+std::map<int, std::vector<double>> GetBands(double const dR[2],
+ int const& numberOfBands,
+ int const& precision,
+ bool const& nearestInteger)
+{
+ auto prec = abs(precision);
+ prec = (prec > 14) ? 14 : prec;
+
+ auto RoundOff = [&prec](const double& x) {
+ auto pow_10 = std::pow(10.0, prec);
+ return std::round(x * pow_10) / pow_10;
+ };
+
+ std::map<int, std::vector<double>> bands;
+ if ((dR[1] < dR[0]) || (numberOfBands <= 0))
+ {
+ return bands;
+ }
+ double x[2];
+ for (int i = 0; i < 2; ++i)
+ {
+ x[i] = dR[i];
+ }
+ if (nearestInteger)
+ {
+ x[0] = std::floor(x[0]);
+ x[1] = std::ceil(x[1]);
+ }
+ double dx = (x[1] - x[0]) / static_cast<double>(numberOfBands);
+ std::vector<double> b;
+ b.push_back(x[0]);
+ b.push_back(x[0] + dx / 2.0);
+ b.push_back(x[0] + dx);
+ for (int i = 0; i < numberOfBands; ++i)
+ {
+ if (i == 0)
+ {
+ for (std::vector<double>::iterator p = b.begin(); p != b.end(); ++p)
+ {
+ *p = RoundOff(*p);
+ }
+ b[0] = x[0];
+ }
+ bands[i] = b;
+ for (std::vector<double>::iterator p = b.begin(); p != b.end(); ++p)
+ {
+ *p = RoundOff(*p + dx);
+ }
+ }
+ return bands;
+}
+
+std::map<int, std::vector<double>>
+GetCustomBands(double const dR[2], int const& numberOfBands,
+ std::vector<std::array<double, 2>> const& myBands)
+{
+ std::map<int, std::vector<double>> bands;
+ if ((dR[1] < dR[0]) || (numberOfBands <= 0))
+ {
+ return bands;
+ }
+
+ std::vector<std::array<double, 2>> x;
+ std::copy(myBands.begin(), myBands.end(), std::back_inserter(x));
+
+ // Determine the index of the range minimum and range maximum.
+ int idxMin = 0;
+ for (auto idx = 0; idx < static_cast<int>(myBands.size()); ++idx)
+ {
+ if (dR[0] < myBands[idx][1] && dR[0] >= myBands[idx][0])
+ {
+ idxMin = idx;
+ break;
+ }
+ }
+ int idxMax = static_cast<int>(myBands.size()) - 1;
+ for (int idx = static_cast<int>(myBands.size()) - 1; idx >= 0; --idx)
+ {
+ if (dR[1] < myBands[idx][1] && dR[1] >= myBands[idx][0])
+ {
+ idxMax = static_cast<int>(idx);
+ break;
+ }
+ }
+
+ // Set the minimum to match the range minimum.
+ x[idxMin][0] = dR[0];
+ x[idxMax][1] = dR[1];
+ for (int i = idxMin; i < idxMax + 1; ++i)
+ {
+ std::vector<double> b(3);
+ b[0] = x[i][0];
+ b[1] = x[i][0] + (x[i][1] - x[i][0]) / 2.0;
+ b[2] = x[i][1];
+ bands[i] = b;
+ }
+ return bands;
+}
+
+std::map<int, std::vector<double>> GetIntegralBands(double const dR[2])
+{
+ std::map<int, std::vector<double>> bands;
+ if (dR[1] < dR[0])
+ {
+ return bands;
+ }
+ double x[2];
+ for (int i = 0; i < 2; ++i)
+ {
+ x[i] = dR[i];
+ }
+ x[0] = std::floor(x[0]);
+ x[1] = std::ceil(x[1]);
+ int numberOfBands = static_cast<int>(std::abs(x[1]) + std::abs(x[0]));
+ return GetBands(x, numberOfBands, false);
+}
+
+std::map<int, int> GetFrequencies(std::map<int, std::vector<double>>& bands,
+ vtkPolyData* src)
{
std::map<int, int> freq;
- for (int i = 0; i < bands.size(); ++i)
+ for (auto i = 0; i < static_cast<int>(bands.size()); ++i)
{
freq[i] = 0;
}
@@ -828,7 +937,7 @@ std::map<int, int> Frequencies(std::map<int, std::vector<double>>& bands,
for (int i = 0; i < tuples; ++i)
{
double* x = src->GetPointData()->GetScalars()->GetTuple(i);
- for (int j = 0; j < bands.size(); ++j)
+ for (auto j = 0; j < static_cast<int>(bands.size()); ++j)
{
if (*x <= bands[j][2])
{
@@ -840,10 +949,12 @@ std::map<int, int> Frequencies(std::map<int, std::vector<double>>& bands,
return freq;
}
-std::pair<int, int> AdjustFrequencyRanges(std::map<int, int>& freq)
+void AdjustFrequencyRanges(std::map<int, std::vector<double>>& bands,
+ std::map<int, int>& freq)
{
+ // Get the indices of the first and last non-zero elements.
auto first = 0;
- for (auto i = 0; i < freq.size(); ++i)
+ for (auto i = 0; i < static_cast<int>(freq.size()); ++i)
{
if (freq[i] != 0)
{
@@ -858,7 +969,7 @@ std::pair<int, int> AdjustFrequencyRanges(std::map<int, int>& freq)
}
std::reverse(keys.begin(), keys.end());
auto last = keys[0];
- for (auto i = 0; i < keys.size(); ++i)
+ for (size_t i = 0; i < keys.size(); ++i)
{
if (freq[keys[i]] != 0)
{
@@ -866,76 +977,56 @@ std::pair<int, int> AdjustFrequencyRanges(std::map<int, int>& freq)
break;
}
}
- return std::pair<int, int>(first, last);
-}
-
-void PrintBands(std::map<int, std::vector<double>> const& bands)
-{
- std::ostringstream os;
- os << "Bands:\n";
- size_t idx = 0;
- for (std::map<int, std::vector<double>>::const_iterator p = bands.begin();
- p != bands.end(); ++p)
+ // Now adjust the ranges.
+ std::map<int, int>::iterator freqItr;
+ freqItr = freq.find(first);
+ freq.erase(freq.begin(), freqItr);
+ freqItr = ++freq.find(last);
+ freq.erase(freqItr, freq.end());
+ std::map<int, std::vector<double>>::iterator bandItr;
+ bandItr = bands.find(first);
+ bands.erase(bands.begin(), bandItr);
+ bandItr = ++bands.find(last);
+ bands.erase(bandItr, bands.end());
+ // Reindex freq and bands.
+ std::map<int, int> adjFreq;
+ int idx = 0;
+ for (auto p : freq)
{
- for (std::vector<double>::const_iterator q = p->second.begin();
- q != p->second.end(); ++q)
- {
- if (q == p->second.begin())
- {
- os << std::setw(4) << idx << " [";
- }
- if (q == std::prev(p->second.end()))
- {
- os << std::fixed << std::setw(6) << std::setprecision(3) << *q << "]\n";
- }
- else
- {
- os << std::fixed << std::setw(6) << std::setprecision(3) << *q << ", ";
- }
- }
+ adjFreq[idx] = p.second;
++idx;
}
- std::cout << os.str() << std::endl;
-}
-
-void PrintFrequencies(std::map<int, int> const& freq)
-{
- std::ostringstream os;
- int i = 0;
- for (std::map<int, int>::const_iterator p = freq.begin(); p != freq.end();
- ++p)
+ std::map<int, std::vector<double>> adjBands;
+ idx = 0;
+ for (auto p : bands)
{
- if (i == 0)
- {
- os << "Frequencies: [";
- }
- if (p == std::prev(freq.end()))
- {
- os << i << ": " << p->second << "]\n";
- }
- else
- {
- os << i << ": " << p->second << ", ";
- }
- ++i;
+ adjBands[idx] = p.second;
+ ++idx;
}
- std::cout << os.str() << endl;
+ bands = adjBands;
+ freq = adjFreq;
}
void PrintBandsFrequencies(std::map<int, std::vector<double>> const& bands,
- std::map<int, int>& freq)
+ std::map<int, int>& freq, int const& precision)
{
+ auto prec = abs(precision);
+ prec = (prec > 14) ? 14 : prec;
+
if (bands.size() != freq.size())
{
std::cout << "Bands and frequencies must be the same size." << std::endl;
return;
}
std::ostringstream os;
- os << "Bands & frequencies:\n";
+ os << "Bands & Frequencies:\n";
size_t idx = 0;
+ auto total = 0;
+ auto width = prec + 6;
for (std::map<int, std::vector<double>>::const_iterator p = bands.begin();
p != bands.end(); ++p)
{
+ total += freq[p->first];
for (std::vector<double>::const_iterator q = p->second.begin();
q != p->second.end(); ++q)
{
@@ -945,57 +1036,21 @@ void PrintBandsFrequencies(std::map<int, std::vector<double>> const& bands,
}
if (q == std::prev(p->second.end()))
{
- os << std::fixed << std::setw(6) << std::setprecision(3) << *q
- << "]: " << std::setw(6) << freq[p->first] << "\n";
+ os << std::fixed << std::setw(width) << std::setprecision(prec) << *q
+ << "]: " << std::setw(8) << freq[p->first] << "\n";
}
else
{
- os << std::fixed << std::setw(6) << std::setprecision(3) << *q << ", ";
+ os << std::fixed << std::setw(width) << std::setprecision(prec) << *q
+ << ", ";
}
}
++idx;
}
+ width = 3 * width + 13;
+ os << std::left << std::setw(width) << "Total" << std::right << std::setw(8)
+ << total << std::endl;
std::cout << os.str() << endl;
}
-vtkNew<vtkGlyph3D> MakeGlyphs(vtkPolyData* src, bool const& reverseNormals)
-{
- // Sometimes the contouring algorithm can create a volume whose gradient
- // vector and ordering of polygon(using the right hand rule) are
- // inconsistent. vtkReverseSense cures this problem.
- vtkNew<vtkReverseSense> reverse;
- vtkNew<vtkMaskPoints> maskPts;
- maskPts->SetOnRatio(5);
- maskPts->RandomModeOn();
- if (reverseNormals)
- {
- reverse->SetInputData(src);
- reverse->ReverseCellsOn();
- reverse->ReverseNormalsOn();
- maskPts->SetInputConnection(reverse->GetOutputPort());
- }
- else
- {
- maskPts->SetInputData(src);
- }
-
- // Source for the glyph filter
- vtkNew<vtkArrowSource> arrow;
- arrow->SetTipResolution(16);
- arrow->SetTipLength(0.3);
- arrow->SetTipRadius(0.1);
-
- vtkNew<vtkGlyph3D> glyph;
- glyph->SetSourceConnection(arrow->GetOutputPort());
- glyph->SetInputConnection(maskPts->GetOutputPort());
- glyph->SetVectorModeToUseNormal();
- glyph->SetScaleFactor(1.0);
- glyph->SetColorModeToColorByVector();
- glyph->SetScaleModeToScaleByVector();
- glyph->OrientOn();
- glyph->Update();
-
- return glyph;
-}
-
} // namespace
diff --git a/src/Python/Utilities/VTKImportsForPython.py b/src/Python/Utilities/VTKImportsForPython.py
index 05ae68119552510f4c893e159b54ce9440032706..6653b8f926db6cdf2b7adf135aa8f24feb72eed2 100755
--- a/src/Python/Utilities/VTKImportsForPython.py
+++ b/src/Python/Utilities/VTKImportsForPython.py
@@ -131,7 +131,7 @@ def main(json_path, src_paths, ofn):
if path.is_file() and path.suffix == '.py':
classes_constants.update(get_classes_constants(path))
elif path.is_dir():
- paths = list(Path(fn).rglob('*.py'))
+ paths = list(Path(fn).glob('*.py'))
program_path = Path(__file__)
for path in paths:
if path.resolve() != program_path.resolve():
@@ -165,7 +165,9 @@ def main(json_path, src_paths, ofn):
res = format_imports(imports)
if ofn:
- path = Path(ofn).with_suffix('.txt')
+ path = Path(ofn)
+ if path.suffix == '':
+ path = Path(ofn).with_suffix('.txt')
path.write_text('\n'.join(res))
else:
print('\n'.join(res))
diff --git a/src/Python/Visualization/CurvatureBandsWithGlyphs.py b/src/Python/Visualization/CurvatureBandsWithGlyphs.py
index 1839b195770820803304e240d9127e3801d089cd..e61b63e81aaf0ab93caf66b4ffa552417880f713 100755
--- a/src/Python/Visualization/CurvatureBandsWithGlyphs.py
+++ b/src/Python/Visualization/CurvatureBandsWithGlyphs.py
@@ -1,7 +1,6 @@
#!/usr/bin/env python
import math
-from collections import OrderedDict
import numpy as np
from vtk.util import numpy_support
@@ -22,7 +21,8 @@ from vtkmodules.vtkCommonCore import (
vtkLookupTable,
vtkPoints,
vtkVariant,
- vtkVariantArray
+ vtkVariantArray,
+ vtkVersion
)
from vtkmodules.vtkCommonDataModel import vtkPolyData
from vtkmodules.vtkCommonTransforms import vtkTransform
@@ -162,7 +162,7 @@ def main(argv):
# Let's do a frequency table.
# The number of scalars in each band.
- freq = frequencies(bands, cc.GetOutput())
+ freq = get_frequencies(bands, cc.GetOutput())
bands, freq = adjust_ranges(bands, freq)
print_bands_frequencies(bands, freq)
@@ -293,10 +293,12 @@ def main(argv):
ren_win.AddRenderer(ren)
# Important: The interactor must be set prior to enabling the widget.
iren.SetRenderWindow(ren_win)
- cam_orient_manipulator = vtkCameraOrientationWidget()
- cam_orient_manipulator.SetParentRenderer(ren)
- # Enable the widget.
- cam_orient_manipulator.On()
+ if vtk_version_ok(9, 0, 20210718):
+ cam_orient_manipulator = vtkCameraOrientationWidget()
+ cam_orient_manipulator = vtkCameraOrientationWidget()
+ cam_orient_manipulator.SetParentRenderer(ren)
+ # Enable the widget.
+ cam_orient_manipulator.On()
# add actors
ren.AddViewProp(src_actor)
@@ -321,83 +323,179 @@ def main(argv):
iren.Start()
-def get_custom_bands(d_r, number_of_bands, my_bands):
+def vtk_version_ok(major, minor, build):
"""
- Divide a range into custom bands.
-
- You need to specify each band as an list [r1, r2] where r1 < r2 and
- append these to a list.
- The list should ultimately look
- like this: [[r1, r2], [r2, r3], [r3, r4]...]
+ Check the VTK version.
- :param: d_r - [min, max] the range that is to be covered by the bands.
- :param: number_of_bands - the number of bands, a positive integer.
- :return: A dictionary consisting of band number and [min, midpoint, max] for each band.
+ :param major: Requested major version.
+ :param minor: Requested minor version.
+ :param build: Requested build version.
+ :return: True if the requested VTK version is >= the actual VTK version.
"""
- bands = dict()
- if (d_r[1] < d_r[0]) or (number_of_bands <= 0):
- return bands
- x = my_bands
- # Determine the index of the range minimum and range maximum.
- idx_min = 0
- for idx in range(0, len(my_bands)):
- if my_bands[idx][1] > d_r[0] >= my_bands[idx][0]:
- idx_min = idx
- break
+ requested_version = (100 * int(major) + int(minor)) * 100000000 + int(build)
+ ver = vtkVersion()
+ actual_version = (100 * ver.GetVTKMajorVersion() + ver.GetVTKMinorVersion()) \
+ * 100000000 + ver.GetVTKBuildVersion()
+ if actual_version >= requested_version:
+ return True
+ else:
+ return False
- idx_max = len(my_bands) - 1
- for idx in range(len(my_bands) - 1, -1, -1):
- if my_bands[idx][1] > d_r[1] >= my_bands[idx][0]:
- idx_max = idx
- break
- # Set the minimum to match the range minimum.
- x[idx_min][0] = d_r[0]
- x[idx_max][1] = d_r[1]
- x = x[idx_min: idx_max + 1]
- for idx, e in enumerate(x):
- bands[idx] = [e[0], e[0] + (e[1] - e[0]) / 2, e[1]]
- return bands
+def adjust_edge_curvatures(source, curvature_name, epsilon=1.0e-08):
+ """
+ This function adjusts curvatures along the edges of the surface by replacing
+ the value with the average value of the curvatures of points in the neighborhood.
+ Remember to update the vtkCurvatures object before calling this.
-def frequencies(bands, src):
- """
- Count the number of scalars in each band.
- :param: bands - the bands.
- :param: src - the vtkPolyData source.
- :return: The frequencies of the scalars in each band.
+ :param source: A vtkPolyData object corresponding to the vtkCurvatures object.
+ :param curvature_name: The name of the curvature, 'Gauss_Curvature' or 'Mean_Curvature'.
+ :param epsilon: Absolute curvature values less than this will be set to zero.
+ :return:
"""
- freq = OrderedDict()
- for i in range(len(bands)):
- freq[i] = 0
- tuples = src.GetPointData().GetScalars().GetNumberOfTuples()
- for i in range(tuples):
- x = src.GetPointData().GetScalars().GetTuple1(i)
- for j in range(len(bands)):
- if x <= bands[j][2]:
- freq[j] = freq[j] + 1
- break
- return freq
+ def point_neighbourhood(pt_id):
+ """
+ Find the ids of the neighbours of pt_id.
+
+ :param pt_id: The point id.
+ :return: The neighbour ids.
+ """
+ """
+ Extract the topological neighbors for point pId. In two steps:
+ 1) source.GetPointCells(pt_id, cell_ids)
+ 2) source.GetCellPoints(cell_id, cell_point_ids) for all cell_id in cell_ids
+ """
+ cell_ids = vtkIdList()
+ source.GetPointCells(pt_id, cell_ids)
+ neighbour = set()
+ for cell_idx in range(0, cell_ids.GetNumberOfIds()):
+ cell_id = cell_ids.GetId(cell_idx)
+ cell_point_ids = vtkIdList()
+ source.GetCellPoints(cell_id, cell_point_ids)
+ for cell_pt_idx in range(0, cell_point_ids.GetNumberOfIds()):
+ neighbour.add(cell_point_ids.GetId(cell_pt_idx))
+ return neighbour
+
+ def compute_distance(pt_id_a, pt_id_b):
+ """
+ Compute the distance between two points given their ids.
+
+ :param pt_id_a:
+ :param pt_id_b:
+ :return:
+ """
+ pt_a = np.array(source.GetPoint(pt_id_a))
+ pt_b = np.array(source.GetPoint(pt_id_b))
+ return np.linalg.norm(pt_a - pt_b)
+
+ # Get the active scalars
+ source.GetPointData().SetActiveScalars(curvature_name)
+ np_source = dsa.WrapDataObject(source)
+ curvatures = np_source.PointData[curvature_name]
+
+ # Get the boundary point IDs.
+ array_name = 'ids'
+ id_filter = vtkIdFilter()
+ id_filter.SetInputData(source)
+ id_filter.SetPointIds(True)
+ id_filter.SetCellIds(False)
+ id_filter.SetPointIdsArrayName(array_name)
+ id_filter.SetCellIdsArrayName(array_name)
+ id_filter.Update()
+
+ edges = vtkFeatureEdges()
+ edges.SetInputConnection(id_filter.GetOutputPort())
+ edges.BoundaryEdgesOn()
+ edges.ManifoldEdgesOff()
+ edges.NonManifoldEdgesOff()
+ edges.FeatureEdgesOff()
+ edges.Update()
+
+ edge_array = edges.GetOutput().GetPointData().GetArray(array_name)
+ boundary_ids = []
+ for i in range(edges.GetOutput().GetNumberOfPoints()):
+ boundary_ids.append(edge_array.GetValue(i))
+ # Remove duplicate Ids.
+ p_ids_set = set(boundary_ids)
-def adjust_frequency_ranges(freq):
+ # Iterate over the edge points and compute the curvature as the weighted
+ # average of the neighbours.
+ count_invalid = 0
+ for p_id in boundary_ids:
+ p_ids_neighbors = point_neighbourhood(p_id)
+ # Keep only interior points.
+ p_ids_neighbors -= p_ids_set
+ # Compute distances and extract curvature values.
+ curvs = [curvatures[p_id_n] for p_id_n in p_ids_neighbors]
+ dists = [compute_distance(p_id_n, p_id) for p_id_n in p_ids_neighbors]
+ curvs = np.array(curvs)
+ dists = np.array(dists)
+ curvs = curvs[dists > 0]
+ dists = dists[dists > 0]
+ if len(curvs) > 0:
+ weights = 1 / np.array(dists)
+ weights /= weights.sum()
+ new_curv = np.dot(curvs, weights)
+ else:
+ # Corner case.
+ count_invalid += 1
+ # Assuming the curvature of the point is planar.
+ new_curv = 0.0
+ # Set the new curvature value.
+ curvatures[p_id] = new_curv
+
+ # Set small values to zero.
+ if epsilon != 0.0:
+ curvatures = np.where(abs(curvatures) < epsilon, 0, curvatures)
+ # Curvatures is now an ndarray
+ curv = numpy_support.numpy_to_vtk(num_array=curvatures.ravel(),
+ deep=True,
+ array_type=VTK_DOUBLE)
+ curv.SetName(curvature_name)
+ source.GetPointData().RemoveArray(curvature_name)
+ source.GetPointData().AddArray(curv)
+ source.GetPointData().SetActiveScalars(curvature_name)
+
+
+def constrain_curvatures(source, curvature_name, lower_bound=0.0, upper_bound=0.0):
"""
- Get the indices of the first and last non-zero elements.
- :param freq: The frequency dictionary.
- :return: The indices of the first and last non-zero elements.
+ This function constrains curvatures to the range [lower_bound ... upper_bound].
+
+ Remember to update the vtkCurvatures object before calling this.
+
+ :param source: A vtkPolyData object corresponding to the vtkCurvatures object.
+ :param curvature_name: The name of the curvature, 'Gauss_Curvature' or 'Mean_Curvature'.
+ :param lower_bound: The lower bound.
+ :param upper_bound: The upper bound.
+ :return:
"""
- first = 0
- for k, v in freq.items():
- if v != 0:
- first = k
- break
- rev_keys = list(freq.keys())[::-1]
- last = rev_keys[0]
- for idx in list(freq.keys())[::-1]:
- if freq[idx] != 0:
- last = idx
- break
- return first, last
+
+ bounds = list()
+ if lower_bound < upper_bound:
+ bounds.append(lower_bound)
+ bounds.append(upper_bound)
+ else:
+ bounds.append(upper_bound)
+ bounds.append(lower_bound)
+
+ # Get the active scalars
+ source.GetPointData().SetActiveScalars(curvature_name)
+ np_source = dsa.WrapDataObject(source)
+ curvatures = np_source.PointData[curvature_name]
+
+ # Set upper and lower bounds.
+ curvatures = np.where(curvatures < bounds[0], bounds[0], curvatures)
+ curvatures = np.where(curvatures > bounds[1], bounds[1], curvatures)
+ # Curvatures is now an ndarray
+ curv = numpy_support.numpy_to_vtk(num_array=curvatures.ravel(),
+ deep=True,
+ array_type=VTK_DOUBLE)
+ curv.SetName(curvature_name)
+ source.GetPointData().RemoveArray(curvature_name)
+ source.GetPointData().AddArray(curv)
+ source.GetPointData().SetActiveScalars(curvature_name)
def get_elevations(src):
@@ -651,160 +749,24 @@ def get_torus():
return cleaner.GetOutput()
-def adjust_edge_curvatures(source, curvature_name, epsilon=1.0e-08):
- """
- This function adjusts curvatures along the edges of the surface by replacing
- the value with the average value of the curvatures of points in the neighborhood.
-
- Remember to update the vtkCurvatures object before calling this.
-
- :param source: A vtkPolyData object corresponding to the vtkCurvatures object.
- :param curvature_name: The name of the curvature, 'Gauss_Curvature' or 'Mean_Curvature'.
- :param epsilon: Absolute curvature values less than this will be set to zero.
- :return:
- """
-
- def point_neighbourhood(pt_id):
- """
- Find the ids of the neighbours of pt_id.
-
- :param pt_id: The point id.
- :return: The neighbour ids.
- """
- """
- Extract the topological neighbors for point pId. In two steps:
- 1) source.GetPointCells(pt_id, cell_ids)
- 2) source.GetCellPoints(cell_id, cell_point_ids) for all cell_id in cell_ids
- """
- cell_ids = vtkIdList()
- source.GetPointCells(pt_id, cell_ids)
- neighbour = set()
- for cell_idx in range(0, cell_ids.GetNumberOfIds()):
- cell_id = cell_ids.GetId(cell_idx)
- cell_point_ids = vtkIdList()
- source.GetCellPoints(cell_id, cell_point_ids)
- for cell_pt_idx in range(0, cell_point_ids.GetNumberOfIds()):
- neighbour.add(cell_point_ids.GetId(cell_pt_idx))
- return neighbour
-
- def compute_distance(pt_id_a, pt_id_b):
- """
- Compute the distance between two points given their ids.
-
- :param pt_id_a:
- :param pt_id_b:
- :return:
- """
- pt_a = np.array(source.GetPoint(pt_id_a))
- pt_b = np.array(source.GetPoint(pt_id_b))
- return np.linalg.norm(pt_a - pt_b)
-
- # Get the active scalars
- source.GetPointData().SetActiveScalars(curvature_name)
- np_source = dsa.WrapDataObject(source)
- curvatures = np_source.PointData[curvature_name]
-
- # Get the boundary point IDs.
- array_name = 'ids'
- id_filter = vtkIdFilter()
- id_filter.SetInputData(source)
- id_filter.SetPointIds(True)
- id_filter.SetCellIds(False)
- id_filter.SetPointIdsArrayName(array_name)
- id_filter.SetCellIdsArrayName(array_name)
- id_filter.Update()
-
- edges = vtkFeatureEdges()
- edges.SetInputConnection(id_filter.GetOutputPort())
- edges.BoundaryEdgesOn()
- edges.ManifoldEdgesOff()
- edges.NonManifoldEdgesOff()
- edges.FeatureEdgesOff()
- edges.Update()
-
- edge_array = edges.GetOutput().GetPointData().GetArray(array_name)
- boundary_ids = []
- for i in range(edges.GetOutput().GetNumberOfPoints()):
- boundary_ids.append(edge_array.GetValue(i))
- # Remove duplicate Ids.
- p_ids_set = set(boundary_ids)
-
- # Iterate over the edge points and compute the curvature as the weighted
- # average of the neighbours.
- count_invalid = 0
- for p_id in boundary_ids:
- p_ids_neighbors = point_neighbourhood(p_id)
- # Keep only interior points.
- p_ids_neighbors -= p_ids_set
- # Compute distances and extract curvature values.
- curvs = [curvatures[p_id_n] for p_id_n in p_ids_neighbors]
- dists = [compute_distance(p_id_n, p_id) for p_id_n in p_ids_neighbors]
- curvs = np.array(curvs)
- dists = np.array(dists)
- curvs = curvs[dists > 0]
- dists = dists[dists > 0]
- if len(curvs) > 0:
- weights = 1 / np.array(dists)
- weights /= weights.sum()
- new_curv = np.dot(curvs, weights)
- else:
- # Corner case.
- count_invalid += 1
- # Assuming the curvature of the point is planar.
- new_curv = 0.0
- # Set the new curvature value.
- curvatures[p_id] = new_curv
-
- # Set small values to zero.
- if epsilon != 0.0:
- curvatures = np.where(abs(curvatures) < epsilon, 0, curvatures)
- # Curvatures is now an ndarray
- curv = numpy_support.numpy_to_vtk(num_array=curvatures.ravel(),
- deep=True,
- array_type=VTK_DOUBLE)
- curv.SetName(curvature_name)
- source.GetPointData().RemoveArray(curvature_name)
- source.GetPointData().AddArray(curv)
- source.GetPointData().SetActiveScalars(curvature_name)
-
-
-def constrain_curvatures(source, curvature_name, lower_bound=0.0, upper_bound=0.0):
- """
- This function constrains curvatures to the range [lower_bound ... upper_bound].
-
- Remember to update the vtkCurvatures object before calling this.
-
- :param source: A vtkPolyData object corresponding to the vtkCurvatures object.
- :param curvature_name: The name of the curvature, 'Gauss_Curvature' or 'Mean_Curvature'.
- :param lower_bound: The lower bound.
- :param upper_bound: The upper bound.
- :return:
- """
-
- bounds = list()
- if lower_bound < upper_bound:
- bounds.append(lower_bound)
- bounds.append(upper_bound)
- else:
- bounds.append(upper_bound)
- bounds.append(lower_bound)
-
- # Get the active scalars
- source.GetPointData().SetActiveScalars(curvature_name)
- np_source = dsa.WrapDataObject(source)
- curvatures = np_source.PointData[curvature_name]
-
- # Set upper and lower bounds.
- curvatures = np.where(curvatures < bounds[0], bounds[0], curvatures)
- curvatures = np.where(curvatures > bounds[1], bounds[1], curvatures)
- # Curvatures is now an ndarray
- curv = numpy_support.numpy_to_vtk(num_array=curvatures.ravel(),
- deep=True,
- array_type=VTK_DOUBLE)
- curv.SetName(curvature_name)
- source.GetPointData().RemoveArray(curvature_name)
- source.GetPointData().AddArray(curv)
- source.GetPointData().SetActiveScalars(curvature_name)
+def get_source(source):
+ surface = source.lower()
+ available_surfaces = ['hills', 'parametrictorus', 'plane', 'randomhills', 'sphere', 'torus']
+ if surface not in available_surfaces:
+ return None
+ elif surface == 'hills':
+ return get_hills()
+ elif surface == 'parametrictorus':
+ return get_parametric_torus()
+ elif surface == 'plane':
+ return get_elevations(get_plane())
+ elif surface == 'randomhills':
+ return get_parametric_hills()
+ elif surface == 'sphere':
+ return get_elevations(get_sphere())
+ elif surface == 'torus':
+ return get_elevations(get_torus())
+ return None
def get_color_series():
@@ -950,26 +912,6 @@ def get_glyphs(src, scale_factor=1.0, reverse_normals=False):
return glyph
-def get_source(source):
- surface = source.lower()
- available_surfaces = ['hills', 'parametrictorus', 'plane', 'randomhills', 'sphere', 'torus']
- if surface not in available_surfaces:
- return None
- elif surface == 'hills':
- return get_hills()
- elif surface == 'parametrictorus':
- return get_parametric_torus()
- elif surface == 'plane':
- return get_elevations(get_plane())
- elif surface == 'randomhills':
- return get_parametric_hills()
- elif surface == 'sphere':
- return get_elevations(get_sphere())
- elif surface == 'torus':
- return get_elevations(get_torus())
- return None
-
-
def get_bands(d_r, number_of_bands, precision=2, nearest_integer=False):
"""
Divide a range into bands
@@ -1003,6 +945,45 @@ def get_bands(d_r, number_of_bands, precision=2, nearest_integer=False):
return bands
+def get_custom_bands(d_r, number_of_bands, my_bands):
+ """
+ Divide a range into custom bands.
+
+ You need to specify each band as an list [r1, r2] where r1 < r2 and
+ append these to a list.
+ The list should ultimately look
+ like this: [[r1, r2], [r2, r3], [r3, r4]...]
+
+ :param: d_r - [min, max] the range that is to be covered by the bands.
+ :param: number_of_bands - the number of bands, a positive integer.
+ :return: A dictionary consisting of band number and [min, midpoint, max] for each band.
+ """
+ bands = dict()
+ if (d_r[1] < d_r[0]) or (number_of_bands <= 0):
+ return bands
+ x = my_bands
+ # Determine the index of the range minimum and range maximum.
+ idx_min = 0
+ for idx in range(0, len(my_bands)):
+ if my_bands[idx][1] > d_r[0] >= my_bands[idx][0]:
+ idx_min = idx
+ break
+
+ idx_max = len(my_bands) - 1
+ for idx in range(len(my_bands) - 1, -1, -1):
+ if my_bands[idx][1] > d_r[1] >= my_bands[idx][0]:
+ idx_max = idx
+ break
+
+ # Set the minimum to match the range minimum.
+ x[idx_min][0] = d_r[0]
+ x[idx_max][1] = d_r[1]
+ x = x[idx_min: idx_max + 1]
+ for idx, e in enumerate(x):
+ bands[idx] = [e[0], e[0] + (e[1] - e[0]) / 2, e[1]]
+ return bands
+
+
def get_frequencies(bands, src):
"""
Count the number of scalars in each band.
diff --git a/src/Python/Visualization/ElevationBandsWithGlyphs.py b/src/Python/Visualization/ElevationBandsWithGlyphs.py
index 8dd6bfb75e64380b0d1a79bb52ab2f68d6b5c96f..594efa4057dbe851eded1f3b42f1d9264bcde34a 100755
--- a/src/Python/Visualization/ElevationBandsWithGlyphs.py
+++ b/src/Python/Visualization/ElevationBandsWithGlyphs.py
@@ -1,57 +1,122 @@
#!/usr/bin/env python
import math
-from collections import OrderedDict
-import vtk
+# noinspection PyUnresolvedReferences
+import vtkmodules.vtkRenderingOpenGL2
+from vtkmodules.vtkCommonColor import (
+ vtkColorSeries,
+ vtkNamedColors
+)
+from vtkmodules.vtkCommonComputationalGeometry import (
+ vtkParametricRandomHills,
+ vtkParametricTorus
+)
+from vtkmodules.vtkCommonCore import (
+ vtkDoubleArray,
+ vtkFloatArray,
+ vtkLookupTable,
+ vtkPoints,
+ vtkVariant,
+ vtkVariantArray,
+ vtkVersion
+)
+from vtkmodules.vtkCommonDataModel import vtkPolyData
+from vtkmodules.vtkCommonTransforms import vtkTransform
+from vtkmodules.vtkFiltersCore import (
+ vtkCleanPolyData,
+ vtkDelaunay2D,
+ vtkElevationFilter,
+ vtkGlyph3D,
+ vtkMaskPoints,
+ vtkPolyDataNormals,
+ vtkReverseSense,
+ vtkTriangleFilter
+)
+from vtkmodules.vtkFiltersGeneral import (
+ vtkTransformPolyDataFilter
+)
+from vtkmodules.vtkFiltersModeling import vtkBandedPolyDataContourFilter
+from vtkmodules.vtkFiltersSources import (
+ vtkArrowSource,
+ vtkParametricFunctionSource,
+ vtkPlaneSource,
+ vtkSphereSource,
+ vtkSuperquadricSource
+)
+from vtkmodules.vtkInteractionStyle import vtkInteractorStyleTrackballCamera
+from vtkmodules.vtkInteractionWidgets import vtkCameraOrientationWidget
+from vtkmodules.vtkRenderingAnnotation import vtkScalarBarActor
+from vtkmodules.vtkRenderingCore import (
+ vtkActor,
+ vtkColorTransferFunction,
+ vtkPolyDataMapper,
+ vtkRenderWindow,
+ vtkRenderWindowInteractor,
+ vtkRenderer
+)
+
+
+def vtk_version_ok(major, minor, build):
+ """
+ Check the VTK version.
+
+ :param major: Requested major version.
+ :param minor: Requested minor version.
+ :param build: Requested build version.
+ :return: True if the requested VTK version is >= the actual VTK version.
+ """
+ requested_version = (100 * int(major) + int(minor)) * 100000000 + int(build)
+ ver = vtkVersion()
+ actual_version = (100 * ver.GetVTKMajorVersion() + ver.GetVTKMinorVersion()) \
+ * 100000000 + ver.GetVTKBuildVersion()
+ if actual_version >= requested_version:
+ return True
+ else:
+ return False
-def main():
+def main(argv):
# ------------------------------------------------------------
# Create the surface, lookup tables, contour filter etc.
# ------------------------------------------------------------
- # desired_surface = 'ParametricTorus'
- # desired_surface = 'Plane'
+ # desired_surface = 'Hills'
+ # desired_surface = 'ParametricTorus'
+ # desired_surface = 'Plane'
desired_surface = 'RandomHills'
- # desired_surface = 'Sphere'
- # desired_surface = 'Torus'
- surface = desired_surface.lower()
- available_surfaces = ['parametrictorus', 'plane', 'randomhills', 'sphere', 'torus']
- if surface not in available_surfaces:
- print('No surface specified.')
- return
- if surface == 'parametrictorus':
- src = make_parametric_torus()
- elif surface == 'plane':
- src = make_elevations(make_plane())
- elif surface == 'randomhills':
- src = make_parametric_hills()
- elif surface == 'sphere':
- src = make_elevations(make_sphere())
- elif surface == 'torus':
- src = make_elevations(make_torus())
- else:
- print('No surface specified.')
+ # desired_surface = 'Sphere'
+ # desired_surface = 'Torus'
+ source = get_source(desired_surface)
+ if not source:
+ print('The surface is not available.')
return
+
+ # The length of the normal arrow glyphs.
+ scale_factor = 1.0
+ if desired_surface == 'Hills':
+ scale_factor = 0.5
+ elif desired_surface == 'Sphere':
+ scale_factor = 2.0
print(desired_surface)
- src.GetPointData().SetActiveScalars('Elevation')
- scalar_range = src.GetPointData().GetScalars('Elevation').GetRange()
+ source.GetPointData().SetActiveScalars('Elevation')
+ scalar_range = source.GetPointData().GetScalars('Elevation').GetRange()
- lut = make_categorical_lut()
- lut1 = make_ordinal_lut()
+ lut = get_categorical_lut()
+ lut1 = get_ordinal_lut()
lut.SetTableRange(scalar_range)
lut1.SetTableRange(scalar_range)
number_of_bands = lut.GetNumberOfTableValues()
lut.SetNumberOfTableValues(number_of_bands)
- bands = make_bands(scalar_range, number_of_bands, False)
+ precision = 10
+ bands = get_bands(scalar_range, number_of_bands, precision, False)
- if surface == 'randomhills':
+ if desired_surface == 'RandomHills':
# These are my custom bands.
# Generated by first running:
- # bands = make_bands(scalar_range, number_of_bands, False)
+ # bands = get_bands(scalar_range, number_of_bands, precision, False)
# then:
- # freq = frequencies(bands, src)
+ # freq = get_frequencies(bands, source)
# print_bands_frequencies(bands, freq)
# Finally using the output to create this table:
my_bands = [
@@ -60,52 +125,29 @@ def main():
[6.0, 7.0], [7.0, 8.0]]
# Comment this out if you want to see how allocating
# equally spaced bands works.
- bands = make_custom_bands(scalar_range, number_of_bands, my_bands)
- # bands = make_bands(scalar_range, number_of_bands, False)
- # Adjust the number of table values
- number_of_bands = len(bands)
- lut.SetNumberOfTableValues(number_of_bands)
- lut1.SetNumberOfTableValues(number_of_bands)
+ bands = get_custom_bands(scalar_range, number_of_bands, my_bands)
+ # bands = get_bands(scalar_range, number_of_bands, precision, False)
# Let's do a frequency table.
# The number of scalars in each band.
- freq = frequencies(bands, src)
+ freq = get_frequencies(bands, source)
+ bands, freq = adjust_ranges(bands, freq)
+ print_bands_frequencies(bands, freq)
- min_key = min(freq.keys())
- max_key = max(freq.keys())
- first, last = adjust_frequency_ranges(freq)
- for idx in range(min_key, first):
- freq.pop(idx)
- bands.pop(idx)
- for idx in range(last + 1, max_key + 1):
- freq.popitem()
- bands.popitem()
- old_keys = freq.keys()
- adj_freq = OrderedDict()
- adj_bands = OrderedDict()
-
- for idx, k in enumerate(old_keys):
- adj_freq[idx] = freq[k]
- adj_bands[idx] = bands[k]
- # print_bands_frequencies(bands, freq)
- print_bands_frequencies(adj_bands, adj_freq)
-
- min_key = min(adj_freq.keys())
- max_key = max(adj_freq.keys())
- scalar_range_curvatures = (adj_bands[min_key][0], adj_bands[max_key][2])
- lut.SetTableRange(scalar_range_curvatures)
- lut.SetNumberOfTableValues(len(adj_bands))
- lut1.SetNumberOfTableValues(len(adj_bands))
+ scalar_range = (bands[0][0], bands[len(bands) - 1][2])
+ lut.SetTableRange(scalar_range)
+ lut.SetNumberOfTableValues(len(bands))
+ lut1.SetNumberOfTableValues(len(bands))
# We will use the midpoint of the band as the label.
labels = []
- for i in range(len(adj_bands)):
- labels.append('{:4.2f}'.format(adj_bands[i][1]))
+ for k in bands:
+ labels.append('{:4.2f}'.format(bands[k][1]))
# Annotate
- values = vtk.vtkVariantArray()
+ values = vtkVariantArray()
for i in range(len(labels)):
- values.InsertNextValue(vtk.vtkVariant(labels[i]))
+ values.InsertNextValue(vtkVariant(labels[i]))
for i in range(values.GetNumberOfTuples()):
lut.SetAnnotation(i, values.GetValue(i).ToString())
@@ -113,47 +155,47 @@ def main():
lutr = reverse_lut(lut)
# Create the contour bands.
- bcf = vtk.vtkBandedPolyDataContourFilter()
- bcf.SetInputData(src)
+ bcf = vtkBandedPolyDataContourFilter()
+ bcf.SetInputData(source)
# Use either the minimum or maximum value for each band.
- for i in range(len(adj_bands)):
- bcf.SetValue(i, adj_bands[i][2])
+ for i in range(len(bands)):
+ bcf.SetValue(i, bands[i][2])
# We will use an indexed lookup table.
bcf.SetScalarModeToIndex()
bcf.GenerateContourEdgesOn()
# Generate the glyphs on the original surface.
- glyph = make_glyphs(src, False)
+ glyph = get_glyphs(source, scale_factor, False)
# ------------------------------------------------------------
# Create the mappers and actors
# ------------------------------------------------------------
- colors = vtk.vtkNamedColors()
+ colors = vtkNamedColors()
# Set the background color.
colors.SetColor('BkgColor', [179, 204, 255, 255])
colors.SetColor("ParaViewBkg", [82, 87, 110, 255])
- src_mapper = vtk.vtkPolyDataMapper()
+ src_mapper = vtkPolyDataMapper()
src_mapper.SetInputConnection(bcf.GetOutputPort())
src_mapper.SetScalarRange(scalar_range)
src_mapper.SetLookupTable(lut)
src_mapper.SetScalarModeToUseCellData()
- src_actor = vtk.vtkActor()
+ src_actor = vtkActor()
src_actor.SetMapper(src_mapper)
# Create contour edges
- edge_mapper = vtk.vtkPolyDataMapper()
+ edge_mapper = vtkPolyDataMapper()
edge_mapper.SetInputData(bcf.GetContourEdgesOutput())
edge_mapper.SetResolveCoincidentTopologyToPolygonOffset()
- edge_actor = vtk.vtkActor()
+ edge_actor = vtkActor()
edge_actor.SetMapper(edge_mapper)
edge_actor.GetProperty().SetColor(colors.GetColor3d('Black'))
- glyph_mapper = vtk.vtkPolyDataMapper()
+ glyph_mapper = vtkPolyDataMapper()
glyph_mapper.SetInputConnection(glyph.GetOutputPort())
glyph_mapper.SetScalarModeToUsePointFieldData()
glyph_mapper.SetColorModeToMapScalars()
@@ -163,14 +205,14 @@ def main():
glyph_mapper.SetLookupTable(lut1)
glyph_mapper.SetScalarRange(scalar_range)
- glyph_actor = vtk.vtkActor()
+ glyph_actor = vtkActor()
glyph_actor.SetMapper(glyph_mapper)
window_width = 800
window_height = 800
# Add a scalar bar.
- scalar_bar = vtk.vtkScalarBarActor()
+ scalar_bar = vtkScalarBarActor()
# This LUT puts the lowest value at the top of the scalar bar.
# scalar_bar->SetLookupTable(lut);
# Use this LUT if you want the highest value at the top.
@@ -190,19 +232,21 @@ def main():
# ------------------------------------------------------------
# Create the RenderWindow, Renderer and Interactor
# ------------------------------------------------------------
- ren = vtk.vtkRenderer()
- ren_win = vtk.vtkRenderWindow()
- iren = vtk.vtkRenderWindowInteractor()
- style = vtk.vtkInteractorStyleTrackballCamera()
+ ren = vtkRenderer()
+ ren_win = vtkRenderWindow()
+ iren = vtkRenderWindowInteractor()
+ style = vtkInteractorStyleTrackballCamera()
iren.SetInteractorStyle(style)
ren_win.AddRenderer(ren)
# Important: The interactor must be set prior to enabling the widget.
iren.SetRenderWindow(ren_win)
- cam_orient_manipulator = vtk.vtkCameraOrientationWidget()
- cam_orient_manipulator.SetParentRenderer(ren)
- # Enable the widget.
- cam_orient_manipulator.On()
+ if vtk_version_ok(9, 0, 20210718):
+ cam_orient_manipulator = vtkCameraOrientationWidget()
+ cam_orient_manipulator = vtkCameraOrientationWidget()
+ cam_orient_manipulator.SetParentRenderer(ren)
+ # Enable the widget.
+ cam_orient_manipulator.On()
# add actors
ren.AddViewProp(src_actor)
@@ -226,111 +270,7 @@ def main():
iren.Start()
-def make_bands(d_r, number_of_bands, nearest_integer):
- """
- Divide a range into bands
- :param: d_r - [min, max] the range that is to be covered by the bands.
- :param: number_of_bands - the number of bands, a positive integer.
- :param: nearest_integer - if True then [floor(min), ceil(max)] is used.
- :return: A dictionary consisting of the band number and [min, midpoint, max] for each band.
- """
- bands = OrderedDict()
- if (d_r[1] < d_r[0]) or (number_of_bands <= 0):
- return bands
- x = list(d_r)
- if nearest_integer:
- x[0] = math.floor(x[0])
- x[1] = math.ceil(x[1])
- dx = (x[1] - x[0]) / float(number_of_bands)
- b = [x[0], x[0] + dx / 2.0, x[0] + dx]
- i = 0
- while i < number_of_bands:
- bands[i] = b
- b = [b[0] + dx, b[1] + dx, b[2] + dx]
- i += 1
- return bands
-
-
-def make_custom_bands(d_r, number_of_bands, my_bands):
- """
- Divide a range into custom bands.
-
- You need to specify each band as an list [r1, r2] where r1 < r2 and
- append these to a list.
- The list should ultimately look
- like this: [[r1, r2], [r2, r3], [r3, r4]...]
-
- :param: d_r - [min, max] the range that is to be covered by the bands.
- :param: number_of_bands - the number of bands, a positive integer.
- :return: A dixtionary consisting of band number and [min, midpoint, max] for each band.
- """
- bands = OrderedDict()
- if (d_r[1] < d_r[0]) or (number_of_bands <= 0):
- return bands
- x = my_bands
- # Determine the index of the range minimum and range maximum.
- idx_min = 0
- for idx in range(0, len(my_bands)):
- if my_bands[idx][1] > d_r[0] >= my_bands[idx][0]:
- idx_min = idx
- break
-
- idx_max = len(my_bands) - 1
- for idx in range(len(my_bands) - 1, -1, -1):
- if my_bands[idx][1] > d_r[1] >= my_bands[idx][0]:
- idx_max = idx
- break
-
- # Set the minimum to match the range minimum.
- x[idx_min][0] = d_r[0]
- x[idx_max][1] = d_r[1]
- x = x[idx_min: idx_max + 1]
- for idx, e in enumerate(x):
- bands[idx] = [e[0], e[0] + (e[1] - e[0]) / 2, e[1]]
- return bands
-
-
-def frequencies(bands, src):
- """
- Count the number of scalars in each band.
- :param: bands - the bands.
- :param: src - the vtkPolyData source.
- :return: The frequencies of the scalars in each band.
- """
- freq = OrderedDict()
- for i in range(len(bands)):
- freq[i] = 0
- tuples = src.GetPointData().GetScalars().GetNumberOfTuples()
- for i in range(tuples):
- x = src.GetPointData().GetScalars().GetTuple1(i)
- for j in range(len(bands)):
- if x <= bands[j][2]:
- freq[j] = freq[j] + 1
- break
- return freq
-
-
-def adjust_frequency_ranges(freq):
- """
- Get the indices of the first and last non-zero elements.
- :param freq: The frequency dictionary.
- :return: The indices of the first and last non-zero elements.
- """
- first = 0
- for k, v in freq.items():
- if v != 0:
- first = k
- break
- rev_keys = list(freq.keys())[::-1]
- last = rev_keys[0]
- for idx in list(freq.keys())[::-1]:
- if freq[idx] != 0:
- last = idx
- break
- return first, last
-
-
-def make_elevations(src):
+def get_elevations(src):
"""
Generate elevations over the surface.
:param: src - the vtkPolyData source.
@@ -340,7 +280,7 @@ def make_elevations(src):
src.GetBounds(bounds)
if abs(bounds[2]) < 1.0e-8 and abs(bounds[3]) < 1.0e-8:
bounds[3] = bounds[2] + 1
- elev_filter = vtk.vtkElevationFilter()
+ elev_filter = vtkElevationFilter()
elev_filter.SetInputData(src)
elev_filter.SetLowPoint(0, bounds[2], 0)
elev_filter.SetHighPoint(0, bounds[3], 0)
@@ -349,12 +289,94 @@ def make_elevations(src):
return elev_filter.GetPolyDataOutput()
-def make_parametric_hills():
+def get_hills():
+ # Create four hills on a plane.
+ # This will have regions of negative, zero and positive Gsaussian curvatures.
+
+ x_res = 50
+ y_res = 50
+ x_min = -5.0
+ x_max = 5.0
+ dx = (x_max - x_min) / (x_res - 1)
+ y_min = -5.0
+ y_max = 5.0
+ dy = (y_max - y_min) / (x_res - 1)
+
+ # Make a grid.
+ points = vtkPoints()
+ for i in range(0, x_res):
+ x = x_min + i * dx
+ for j in range(0, y_res):
+ y = y_min + j * dy
+ points.InsertNextPoint(x, y, 0)
+
+ # Add the grid points to a polydata object.
+ plane = vtkPolyData()
+ plane.SetPoints(points)
+
+ # Triangulate the grid.
+ delaunay = vtkDelaunay2D()
+ delaunay.SetInputData(plane)
+ delaunay.Update()
+
+ polydata = delaunay.GetOutput()
+
+ elevation = vtkDoubleArray()
+ elevation.SetNumberOfTuples(points.GetNumberOfPoints())
+
+ # We define the parameters for the hills here.
+ # [[0: x0, 1: y0, 2: x variance, 3: y variance, 4: amplitude]...]
+ hd = [[-2.5, -2.5, 2.5, 6.5, 3.5], [2.5, 2.5, 2.5, 2.5, 2],
+ [5.0, -2.5, 1.5, 1.5, 2.5], [-5.0, 5, 2.5, 3.0, 3]]
+ xx = [0.0] * 2
+ for i in range(0, points.GetNumberOfPoints()):
+ x = list(polydata.GetPoint(i))
+ for j in range(0, len(hd)):
+ xx[0] = (x[0] - hd[j][0] / hd[j][2]) ** 2.0
+ xx[1] = (x[1] - hd[j][1] / hd[j][3]) ** 2.0
+ x[2] += hd[j][4] * math.exp(-(xx[0] + xx[1]) / 2.0)
+ polydata.GetPoints().SetPoint(i, x)
+ elevation.SetValue(i, x[2])
+
+ textures = vtkFloatArray()
+ textures.SetNumberOfComponents(2)
+ textures.SetNumberOfTuples(2 * polydata.GetNumberOfPoints())
+ textures.SetName("Textures")
+
+ for i in range(0, x_res):
+ tc = [i / (x_res - 1.0), 0.0]
+ for j in range(0, y_res):
+ # tc[1] = 1.0 - j / (y_res - 1.0)
+ tc[1] = j / (y_res - 1.0)
+ textures.SetTuple(i * y_res + j, tc)
+
+ polydata.GetPointData().SetScalars(elevation)
+ polydata.GetPointData().GetScalars().SetName("Elevation")
+ polydata.GetPointData().SetTCoords(textures)
+
+ normals = vtkPolyDataNormals()
+ normals.SetInputData(polydata)
+ normals.SetInputData(polydata)
+ normals.SetFeatureAngle(30)
+ normals.SplittingOff()
+
+ tr1 = vtkTransform()
+ tr1.RotateX(-90)
+
+ tf1 = vtkTransformPolyDataFilter()
+ tf1.SetInputConnection(normals.GetOutputPort())
+ tf1.SetTransform(tr1)
+ tf1.Update()
+
+ return tf1.GetOutput()
+
+
+def get_parametric_hills():
"""
Make a parametric hills surface as the source.
:return: vtkPolyData with normal and scalar data.
"""
- fn = vtk.vtkParametricRandomHills()
+ fn = vtkParametricRandomHills()
fn.AllowRandomGenerationOn()
fn.SetRandomSeed(1)
fn.SetNumberOfHills(30)
@@ -363,7 +385,7 @@ def make_parametric_hills():
# if fn.GetClassName() == 'vtkParametricRandomHills':
# fn.ClockwiseOrderingOff()
- source = vtk.vtkParametricFunctionSource()
+ source = vtkParametricFunctionSource()
source.SetParametricFunction(fn)
source.SetUResolution(50)
source.SetVResolution(50)
@@ -375,10 +397,10 @@ def make_parametric_hills():
# Rename the scalars to 'Elevation' since we are using the Z-scalars as elevations.
source.GetOutput().GetPointData().GetScalars().SetName('Elevation')
- transform = vtk.vtkTransform()
+ transform = vtkTransform()
transform.Translate(0.0, 5.0, 15.0)
transform.RotateX(-90.0)
- transform_filter = vtk.vtkTransformPolyDataFilter()
+ transform_filter = vtkTransformPolyDataFilter()
transform_filter.SetInputConnection(source.GetOutputPort())
transform_filter.SetTransform(transform)
transform_filter.Update()
@@ -386,17 +408,17 @@ def make_parametric_hills():
return transform_filter.GetOutput()
-def make_parametric_torus():
+def get_parametric_torus():
"""
Make a parametric torus as the source.
:return: vtkPolyData with normal and scalar data.
"""
- fn = vtk.vtkParametricTorus()
+ fn = vtkParametricTorus()
fn.SetRingRadius(5)
fn.SetCrossSectionRadius(2)
- source = vtk.vtkParametricFunctionSource()
+ source = vtkParametricFunctionSource()
source.SetParametricFunction(fn)
source.SetUResolution(50)
source.SetVResolution(50)
@@ -409,9 +431,9 @@ def make_parametric_torus():
# Rename the scalars to 'Elevation' since we are using the Z-scalars as elevations.
source.GetOutput().GetPointData().GetScalars().SetName('Elevation')
- transform = vtk.vtkTransform()
+ transform = vtkTransform()
transform.RotateX(-90.0)
- transform_filter = vtk.vtkTransformPolyDataFilter()
+ transform_filter = vtkTransformPolyDataFilter()
transform_filter.SetInputConnection(source.GetOutputPort())
transform_filter.SetTransform(transform)
transform_filter.Update()
@@ -419,13 +441,13 @@ def make_parametric_torus():
return transform_filter.GetOutput()
-def make_plane():
+def get_plane():
"""
Make a plane as the source.
:return: vtkPolyData with normal and scalar data.
"""
- source = vtk.vtkPlaneSource()
+ source = vtkPlaneSource()
source.SetOrigin(-10.0, -10.0, 0.0)
source.SetPoint2(-10.0, 10.0, 0.0)
source.SetPoint1(10.0, -10.0, 0.0)
@@ -433,10 +455,10 @@ def make_plane():
source.SetYResolution(20)
source.Update()
- transform = vtk.vtkTransform()
+ transform = vtkTransform()
transform.Translate(0.0, 0.0, 0.0)
transform.RotateX(-90.0)
- transform_filter = vtk.vtkTransformPolyDataFilter()
+ transform_filter = vtkTransformPolyDataFilter()
transform_filter.SetInputConnection(source.GetOutputPort())
transform_filter.SetTransform(transform)
transform_filter.Update()
@@ -444,12 +466,12 @@ def make_plane():
# We have a m x n array of quadrilaterals arranged as a regular tiling in a
# plane. So pass it through a triangle filter since the curvature filter only
# operates on polys.
- tri = vtk.vtkTriangleFilter()
+ tri = vtkTriangleFilter()
tri.SetInputConnection(transform_filter.GetOutputPort())
# Pass it though a CleanPolyDataFilter and merge any points which
# are coincident, or very close
- cleaner = vtk.vtkCleanPolyData()
+ cleaner = vtkCleanPolyData()
cleaner.SetInputConnection(tri.GetOutputPort())
cleaner.SetTolerance(0.005)
cleaner.Update()
@@ -457,8 +479,8 @@ def make_plane():
return cleaner.GetOutput()
-def make_sphere():
- source = vtk.vtkSphereSource()
+def get_sphere():
+ source = vtkSphereSource()
source.SetCenter(0.0, 0.0, 0.0)
source.SetRadius(10.0)
source.SetThetaResolution(32)
@@ -468,12 +490,12 @@ def make_sphere():
return source.GetOutput()
-def make_torus():
+def get_torus():
"""
Make a torus as the source.
:return: vtkPolyData with normal and scalar data.
"""
- source = vtk.vtkSuperquadricSource()
+ source = vtkSuperquadricSource()
source.SetCenter(0.0, 0.0, 0.0)
source.SetScale(1.0, 1.0, 1.0)
source.SetPhiResolution(64)
@@ -485,13 +507,13 @@ def make_torus():
# The quadric is made of strips, so pass it through a triangle filter as
# the curvature filter only operates on polys
- tri = vtk.vtkTriangleFilter()
+ tri = vtkTriangleFilter()
tri.SetInputConnection(source.GetOutputPort())
# The quadric has nasty discontinuities from the way the edges are generated
# so let's pass it though a CleanPolyDataFilter and merge any points which
# are coincident, or very close
- cleaner = vtk.vtkCleanPolyData()
+ cleaner = vtkCleanPolyData()
cleaner.SetInputConnection(tri.GetOutputPort())
cleaner.SetTolerance(0.005)
cleaner.Update()
@@ -499,80 +521,28 @@ def make_torus():
return cleaner.GetOutput()
-def clipper(src, dx, dy, dz):
- """
- Clip a vtkPolyData source.
- A cube is made whose size corresponds the the bounds of the source.
- Then each side is shrunk by the appropriate dx, dy or dz. After
- this operation the source is clipped by the cube.
- :param: src - the vtkPolyData source
- :param: dx - the amount to clip in the x-direction
- :param: dy - the amount to clip in the y-direction
- :param: dz - the amount to clip in the z-direction
- :return: vtkPolyData.
- """
- bounds = [0, 0, 0, 0, 0, 0]
- src.GetBounds(bounds)
-
- plane1 = vtk.vtkPlane()
- plane1.SetOrigin(bounds[0] + dx, 0, 0)
- plane1.SetNormal(1, 0, 0)
-
- plane2 = vtk.vtkPlane()
- plane2.SetOrigin(bounds[1] - dx, 0, 0)
- plane2.SetNormal(-1, 0, 0)
-
- plane3 = vtk.vtkPlane()
- plane3.SetOrigin(0, bounds[2] + dy, 0)
- plane3.SetNormal(0, 1, 0)
-
- plane4 = vtk.vtkPlane()
- plane4.SetOrigin(0, bounds[3] - dy, 0)
- plane4.SetNormal(0, -1, 0)
-
- plane5 = vtk.vtkPlane()
- plane5.SetOrigin(0, 0, bounds[4] + dz)
- plane5.SetNormal(0, 0, 1)
-
- plane6 = vtk.vtkPlane()
- plane6.SetOrigin(0, 0, bounds[5] - dz)
- plane6.SetNormal(0, 0, -1)
-
- clip_function = vtk.vtkImplicitBoolean()
- clip_function.SetOperationTypeToUnion()
- clip_function.AddFunction(plane1)
- clip_function.AddFunction(plane2)
- clip_function.AddFunction(plane3)
- clip_function.AddFunction(plane4)
- clip_function.AddFunction(plane5)
- clip_function.AddFunction(plane6)
-
- # Clip it.
- pd_clipper = vtk.vtkClipPolyData()
- pd_clipper.SetClipFunction(clip_function)
- pd_clipper.SetInputData(src)
- pd_clipper.GenerateClipScalarsOff()
- pd_clipper.GenerateClippedOutputOff()
- # pd_clipper.GenerateClippedOutputOn()
- pd_clipper.Update()
- return pd_clipper.GetOutput()
-
-
-def calculate_curvatures(src):
- """
- The source must be triangulated.
- :param: src - the source.
- :return: vtkPolyData with normal and scalar data representing curvatures.
- """
- curvature = vtk.vtkCurvatures()
- curvature.SetCurvatureTypeToGaussian()
- curvature.SetInputData(src)
- curvature.Update()
- return curvature.GetOutput()
+def get_source(source):
+ surface = source.lower()
+ available_surfaces = ['hills', 'parametrictorus', 'plane', 'randomhills', 'sphere', 'torus']
+ if surface not in available_surfaces:
+ return None
+ elif surface == 'hills':
+ return get_hills()
+ elif surface == 'parametrictorus':
+ return get_parametric_torus()
+ elif surface == 'plane':
+ return get_elevations(get_plane())
+ elif surface == 'randomhills':
+ return get_parametric_hills()
+ elif surface == 'sphere':
+ return get_elevations(get_sphere())
+ elif surface == 'torus':
+ return get_elevations(get_torus())
+ return None
def get_color_series():
- color_series = vtk.vtkColorSeries()
+ color_series = vtkColorSeries()
# Select a color scheme.
# color_series_enum = color_series.BREWER_DIVERGING_BROWN_BLUE_GREEN_9
# color_series_enum = color_series.BREWER_DIVERGING_SPECTRAL_10
@@ -586,39 +556,71 @@ def get_color_series():
return color_series
-def make_categorical_lut():
+def get_categorical_lut():
"""
Make a lookup table using vtkColorSeries.
:return: An indexed (categorical) lookup table.
"""
color_series = get_color_series()
# Make the lookup table.
- lut = vtk.vtkLookupTable()
+ lut = vtkLookupTable()
color_series.BuildLookupTable(lut, color_series.CATEGORICAL)
lut.SetNanColor(0, 0, 0, 1)
return lut
-def make_ordinal_lut():
+def get_ordinal_lut():
"""
Make a lookup table using vtkColorSeries.
:return: An ordinal (not indexed) lookup table.
"""
color_series = get_color_series()
# Make the lookup table.
- lut = vtk.vtkLookupTable()
+ lut = vtkLookupTable()
color_series.BuildLookupTable(lut, color_series.ORDINAL)
lut.SetNanColor(0, 0, 0, 1)
return lut
+def get_diverging_lut():
+ """
+ See: [Diverging Color Maps for Scientific Visualization](https://www.kennethmoreland.com/color-maps/)
+ start point midPoint end point
+ cool to warm: 0.230, 0.299, 0.754 0.865, 0.865, 0.865 0.706, 0.016, 0.150
+ purple to orange: 0.436, 0.308, 0.631 0.865, 0.865, 0.865 0.759, 0.334, 0.046
+ green to purple: 0.085, 0.532, 0.201 0.865, 0.865, 0.865 0.436, 0.308, 0.631
+ blue to brown: 0.217, 0.525, 0.910 0.865, 0.865, 0.865 0.677, 0.492, 0.093
+ green to red: 0.085, 0.532, 0.201 0.865, 0.865, 0.865 0.758, 0.214, 0.233
+
+ :return:
+ """
+ ctf = vtkColorTransferFunction()
+ ctf.SetColorSpaceToDiverging()
+ # Cool to warm.
+ ctf.AddRGBPoint(0.0, 0.085, 0.532, 0.201)
+ ctf.AddRGBPoint(0.5, 0.865, 0.865, 0.865)
+ ctf.AddRGBPoint(1.0, 0.758, 0.214, 0.233)
+
+ table_size = 256
+ lut = vtkLookupTable()
+ lut.SetNumberOfTableValues(table_size)
+ lut.Build()
+
+ for i in range(0, table_size):
+ rgba = list(ctf.GetColor(float(i) / table_size))
+ rgba.append(1)
+ lut.SetTableValue(i, rgba)
+
+ return lut
+
+
def reverse_lut(lut):
"""
Create a lookup table with the colors reversed.
:param: lut - An indexed lookup table.
:return: The reversed indexed lookup table.
"""
- lutr = vtk.vtkLookupTable()
+ lutr = vtkLookupTable()
lutr.DeepCopy(lut)
t = lut.GetNumberOfTableValues() - 1
rev_range = reversed(list(range(t + 1)))
@@ -635,7 +637,7 @@ def reverse_lut(lut):
return lutr
-def make_glyphs(src, reverse_normals):
+def get_glyphs(src, scale_factor=1.0, reverse_normals=False):
"""
Glyph the normals on the surface.
@@ -650,10 +652,10 @@ def make_glyphs(src, reverse_normals):
# Sometimes the contouring algorithm can create a volume whose gradient
# vector and ordering of polygon (using the right hand rule) are
# inconsistent. vtkReverseSense cures this problem.
- reverse = vtk.vtkReverseSense()
+ reverse = vtkReverseSense()
# Choose a random subset of points.
- mask_pts = vtk.vtkMaskPoints()
+ mask_pts = vtkMaskPoints()
mask_pts.SetOnRatio(5)
mask_pts.RandomModeOn()
if reverse_normals:
@@ -665,16 +667,16 @@ def make_glyphs(src, reverse_normals):
mask_pts.SetInputData(src)
# Source for the glyph filter
- arrow = vtk.vtkArrowSource()
+ arrow = vtkArrowSource()
arrow.SetTipResolution(16)
arrow.SetTipLength(0.3)
arrow.SetTipRadius(0.1)
- glyph = vtk.vtkGlyph3D()
+ glyph = vtkGlyph3D()
glyph.SetSourceConnection(arrow.GetOutputPort())
glyph.SetInputConnection(mask_pts.GetOutputPort())
glyph.SetVectorModeToUseNormal()
- glyph.SetScaleFactor(1.0)
+ glyph.SetScaleFactor(scale_factor)
glyph.SetColorModeToColorByVector()
glyph.SetScaleModeToScaleByVector()
glyph.OrientOn()
@@ -682,46 +684,166 @@ def make_glyphs(src, reverse_normals):
return glyph
-def print_bands(bands):
- s = f'Bands:\n'
- for k, v in bands.items():
- for j, q in enumerate(v):
- if j == 0:
- s += f'{k:4d} ['
- if j == len(v) - 1:
- s += f'{q:8.3f}]\n'
- else:
- s += f'{q:8.3f}, '
- print(s)
-
+def get_bands(d_r, number_of_bands, precision=2, nearest_integer=False):
+ """
+ Divide a range into bands
+ :param: d_r - [min, max] the range that is to be covered by the bands.
+ :param: number_of_bands - The number of bands, a positive integer.
+ :param: precision - The decimal precision of the bounds.
+ :param: nearest_integer - If True then [floor(min), ceil(max)] is used.
+ :return: A dictionary consisting of the band number and [min, midpoint, max] for each band.
+ """
+ prec = abs(precision)
+ if prec > 14:
+ prec = 14
-def print_frequencies(freq):
- s = ''
- for i, p in freq.items():
+ bands = dict()
+ if (d_r[1] < d_r[0]) or (number_of_bands <= 0):
+ return bands
+ x = list(d_r)
+ if nearest_integer:
+ x[0] = math.floor(x[0])
+ x[1] = math.ceil(x[1])
+ dx = (x[1] - x[0]) / float(number_of_bands)
+ b = [x[0], x[0] + dx / 2.0, x[0] + dx]
+ i = 0
+ while i < number_of_bands:
+ b = list(map(lambda ele_b: round(ele_b, prec), b))
if i == 0:
- s += f'Frequencies: ['
- if i == len(freq) - 1:
- s += f'{i}: {p} ]'
- else:
- s += f'{i}: {p}, '
- print(s)
+ b[0] = x[0]
+ bands[i] = b
+ b = [b[0] + dx, b[1] + dx, b[2] + dx]
+ i += 1
+ return bands
-def print_bands_frequencies(bands, freq):
+def get_custom_bands(d_r, number_of_bands, my_bands):
+ """
+ Divide a range into custom bands.
+
+ You need to specify each band as an list [r1, r2] where r1 < r2 and
+ append these to a list.
+ The list should ultimately look
+ like this: [[r1, r2], [r2, r3], [r3, r4]...]
+
+ :param: d_r - [min, max] the range that is to be covered by the bands.
+ :param: number_of_bands - the number of bands, a positive integer.
+ :return: A dictionary consisting of band number and [min, midpoint, max] for each band.
+ """
+ bands = dict()
+ if (d_r[1] < d_r[0]) or (number_of_bands <= 0):
+ return bands
+ x = my_bands
+ # Determine the index of the range minimum and range maximum.
+ idx_min = 0
+ for idx in range(0, len(my_bands)):
+ if my_bands[idx][1] > d_r[0] >= my_bands[idx][0]:
+ idx_min = idx
+ break
+
+ idx_max = len(my_bands) - 1
+ for idx in range(len(my_bands) - 1, -1, -1):
+ if my_bands[idx][1] > d_r[1] >= my_bands[idx][0]:
+ idx_max = idx
+ break
+
+ # Set the minimum to match the range minimum.
+ x[idx_min][0] = d_r[0]
+ x[idx_max][1] = d_r[1]
+ x = x[idx_min: idx_max + 1]
+ for idx, e in enumerate(x):
+ bands[idx] = [e[0], e[0] + (e[1] - e[0]) / 2, e[1]]
+ return bands
+
+
+def get_frequencies(bands, src):
+ """
+ Count the number of scalars in each band.
+ The scalars used are the active scalars in the polydata.
+
+ :param: bands - The bands.
+ :param: src - The vtkPolyData source.
+ :return: The frequencies of the scalars in each band.
+ """
+ freq = dict()
+ for i in range(len(bands)):
+ freq[i] = 0
+ tuples = src.GetPointData().GetScalars().GetNumberOfTuples()
+ for i in range(tuples):
+ x = src.GetPointData().GetScalars().GetTuple1(i)
+ for j in range(len(bands)):
+ if x <= bands[j][2]:
+ freq[j] += 1
+ break
+ return freq
+
+
+def adjust_ranges(bands, freq):
+ """
+ The bands and frequencies are adjusted so that the first and last
+ frequencies in the range are non-zero.
+ :param bands: The bands dictionary.
+ :param freq: The frequency dictionary.
+ :return: Adjusted bands and frequencies.
+ """
+ # Get the indices of the first and last non-zero elements.
+ first = 0
+ for k, v in freq.items():
+ if v != 0:
+ first = k
+ break
+ rev_keys = list(freq.keys())[::-1]
+ last = rev_keys[0]
+ for idx in list(freq.keys())[::-1]:
+ if freq[idx] != 0:
+ last = idx
+ break
+ # Now adjust the ranges.
+ min_key = min(freq.keys())
+ max_key = max(freq.keys())
+ for idx in range(min_key, first):
+ freq.pop(idx)
+ bands.pop(idx)
+ for idx in range(last + 1, max_key + 1):
+ freq.popitem()
+ bands.popitem()
+ old_keys = freq.keys()
+ adj_freq = dict()
+ adj_bands = dict()
+
+ for idx, k in enumerate(old_keys):
+ adj_freq[idx] = freq[k]
+ adj_bands[idx] = bands[k]
+
+ return adj_bands, adj_freq
+
+
+def print_bands_frequencies(bands, freq, precision=2):
+ prec = abs(precision)
+ if prec > 14:
+ prec = 14
+
if len(bands) != len(freq):
- print('Bands and frequencies must be the same size.')
+ print('Bands and Frequencies must be the same size.')
return
- s = f'Bands & frequencies:\n'
+ s = f'Bands & Frequencies:\n'
+ total = 0
+ width = prec + 6
for k, v in bands.items():
+ total += freq[k]
for j, q in enumerate(v):
if j == 0:
s += f'{k:4d} ['
if j == len(v) - 1:
- s += f'{q:8.3f}]: {freq[k]:6d}\n'
+ s += f'{q:{width}.{prec}f}]: {freq[k]:8d}\n'
else:
- s += f'{q:8.3f}, '
+ s += f'{q:{width}.{prec}f}, '
+ width = 3 * width + 13
+ s += f'{"Total":{width}s}{total:8d}\n'
print(s)
if __name__ == '__main__':
- main()
+ import sys
+
+ main(sys.argv)