From 49601a2f465d6b3b4ba6d6cec4265af0fb7efd73 Mon Sep 17 00:00:00 2001 From: amaclean Date: Tue, 9 Jun 2026 09:18:37 +1000 Subject: [PATCH 1/5] Using an namedtuple instead of a dataclass, also added a glyph type field --- .../CurvatureBandsWithGlyphs.cxx | 3 ++ .../ElevationBandsWithGlyphs.cxx | 2 + .../Visualization/CurvatureBandsWithGlyphs.py | 39 +++++++------------ .../Visualization/ElevationBandsWithGlyphs.py | 20 ++++------ 4 files changed, 25 insertions(+), 39 deletions(-) diff --git a/src/Cxx/Visualization/CurvatureBandsWithGlyphs.cxx b/src/Cxx/Visualization/CurvatureBandsWithGlyphs.cxx index 64f3f467fad..6c5317a3441 100644 --- a/src/Cxx/Visualization/CurvatureBandsWithGlyphs.cxx +++ b/src/Cxx/Visualization/CurvatureBandsWithGlyphs.cxx @@ -77,6 +77,7 @@ vtkNew ReverseLUT(vtkLookupTable* lut); struct BandedGlyphs { + std::string glyphType; vtkSmartPointer glyphs; vtkSmartPointer bcf; vtkSmartPointer lut; @@ -1060,6 +1061,7 @@ BandedGlyphs GetCurvatureGlyphs(std::string const& surfaceName, auto glyphs = GetGlyphs(source, scaleFactor, false); BandedGlyphs bg; + bg.glyphType = curvature; bg.glyphs = glyphs; bg.bcf = bcf; bg.lut = lut; @@ -1189,6 +1191,7 @@ BandedGlyphs GetElevationGlyphs(std::string const& surfaceName, auto glyphs = GetGlyphs(source, scaleFactor, false); BandedGlyphs bg; + bg.glyphType = "Elevation"; bg.glyphs = glyphs; bg.bcf = bcf; bg.lut = lut; diff --git a/src/Cxx/Visualization/ElevationBandsWithGlyphs.cxx b/src/Cxx/Visualization/ElevationBandsWithGlyphs.cxx index dedca86b48b..e04e3d4615e 100644 --- a/src/Cxx/Visualization/ElevationBandsWithGlyphs.cxx +++ b/src/Cxx/Visualization/ElevationBandsWithGlyphs.cxx @@ -68,6 +68,7 @@ vtkSmartPointer ReverseLUT(vtkLookupTable* lut); struct BandedGlyphs { + std::string glyphType; vtkSmartPointer glyphs; vtkSmartPointer bcf; vtkSmartPointer lut; @@ -653,6 +654,7 @@ BandedGlyphs GetElevationGlyphs(std::string const& surfaceName, auto glyphs = GetGlyphs(source, scaleFactor, false); BandedGlyphs bg; + bg.glyphType = "Elevation"; bg.glyphs = glyphs; bg.bcf = bcf; bg.lut = lut; diff --git a/src/PythonicAPI/Visualization/CurvatureBandsWithGlyphs.py b/src/PythonicAPI/Visualization/CurvatureBandsWithGlyphs.py index 99cc5d54f33..d1eed81cdd6 100644 --- a/src/PythonicAPI/Visualization/CurvatureBandsWithGlyphs.py +++ b/src/PythonicAPI/Visualization/CurvatureBandsWithGlyphs.py @@ -2,6 +2,7 @@ import copy import math +from collections import namedtuple from dataclasses import dataclass import numpy as np @@ -983,7 +984,12 @@ def get_curvature_glyphs(surface_name, source, curvature, glyphs = get_glyphs(source, scale_factor, reverse_normals=False) - bg = CurvatureBandedGlyphs + BandedGlyphs = namedtuple('BandedGlyphs', + ['glyph_type', 'glyphs', 'bcf', + 'lut', 'lutr', 'lut1', 'lut1r', + 'scalar_range', 'labels']) + bg = BandedGlyphs + bg.glyph_type = curvature bg.glyphs = glyphs bg.bcf = bcf bg.lut = lut @@ -996,18 +1002,6 @@ def get_curvature_glyphs(surface_name, source, curvature, return bg -@dataclass -class CurvatureBandedGlyphs: - glyphs: vtkGlyph3D - bcf: vtkBandedPolyDataContourFilter - lut: vtkLookupTable - lutr = vtkLookupTable - lut1: vtkLookupTable - lut1r: vtkLookupTable - scalar_range: tuple - labels: list - - def get_elevation_glyphs(surface_name, source, precision, frequency_table=False, nearest_integer=False): """ @@ -1104,7 +1098,12 @@ def get_elevation_glyphs(surface_name, source, glyphs = get_glyphs(source, scale_factor, reverse_normals=False) - bg = ElevationBandedGlyphs + BandedGlyphs = namedtuple('BandedGlyphs', + ['glyph_type', 'glyphs', 'bcf', + 'lut', 'lutr', 'lut1', 'lut1r', + 'scalar_range', 'labels']) + bg = BandedGlyphs + bg.glyph_type = 'Elevation' bg.glyphs = glyphs bg.bcf = bcf bg.lut = lut @@ -1117,18 +1116,6 @@ def get_elevation_glyphs(surface_name, source, return bg -@dataclass -class ElevationBandedGlyphs: - glyphs: vtkGlyph3D - bcf: vtkBandedPolyDataContourFilter - lut: vtkLookupTable - lutr = vtkLookupTable - lut1: vtkLookupTable - lut1r: vtkLookupTable - scalar_range: tuple - labels: list - - class ScalarBarProperties: """ The properties needed for scalar bars. diff --git a/src/PythonicAPI/Visualization/ElevationBandsWithGlyphs.py b/src/PythonicAPI/Visualization/ElevationBandsWithGlyphs.py index 4dd06aa7430..dd822a8efab 100644 --- a/src/PythonicAPI/Visualization/ElevationBandsWithGlyphs.py +++ b/src/PythonicAPI/Visualization/ElevationBandsWithGlyphs.py @@ -2,6 +2,7 @@ import copy import math +from collections import namedtuple from dataclasses import dataclass # noinspection PyUnresolvedReferences @@ -751,7 +752,12 @@ def get_elevation_glyphs(surface_name, source, glyphs = get_glyphs(source, scale_factor, reverse_normals=False) - bg = ElevationBandedGlyphs + BandedGlyphs = namedtuple('BandedGlyphs', + ['glyph_type', 'glyphs', 'bcf', + 'lut', 'lutr', 'lut1', 'lut1r', + 'scalar_range', 'labels']) + bg = BandedGlyphs + bg.glyph_type = 'Elevation' bg.glyphs = glyphs bg.bcf = bcf bg.lut = lut @@ -764,18 +770,6 @@ def get_elevation_glyphs(surface_name, source, return bg -@dataclass -class ElevationBandedGlyphs: - glyphs: vtkGlyph3D - bcf: vtkBandedPolyDataContourFilter - lut: vtkLookupTable - lutr = vtkLookupTable - lut1: vtkLookupTable - lut1r: vtkLookupTable - scalar_range: tuple - labels: list - - class ScalarBarProperties: """ The properties needed for scalar bars. -- GitLab From a73918a924989f6bebb38df568c081f19b3a6ad7 Mon Sep 17 00:00:00 2001 From: amaclean Date: Thu, 11 Jun 2026 15:05:59 +1000 Subject: [PATCH 2/5] Updated the C++ and PythonicAPI versions --- .../CurvatureBandsWithGlyphs.cxx | 2131 +++++++++-------- .../ElevationBandsWithGlyphs.cxx | 463 ++-- .../PolyData/CurvaturesNormalsElevations.py | 1549 ++++++------ .../Visualization/CurvatureBandsWithGlyphs.py | 362 +-- .../Visualization/ElevationBandsWithGlyphs.py | 28 +- 5 files changed, 2341 insertions(+), 2192 deletions(-) diff --git a/src/Cxx/Visualization/CurvatureBandsWithGlyphs.cxx b/src/Cxx/Visualization/CurvatureBandsWithGlyphs.cxx index 6c5317a3441..cb6682e852d 100644 --- a/src/Cxx/Visualization/CurvatureBandsWithGlyphs.cxx +++ b/src/Cxx/Visualization/CurvatureBandsWithGlyphs.cxx @@ -65,6 +65,13 @@ namespace fs = std::filesystem; namespace { +/** + * Generate elevations over the surface. + * @param src - the vtkPolyData source. + * @return elev - the elevations. + */ +vtkSmartPointer GenerateElevations(vtkPolyData* src); + vtkSmartPointer GetHills(); vtkSmartPointer GetParametricHills(); vtkSmartPointer GetParametricTorus(); @@ -75,101 +82,6 @@ vtkSmartPointer GetSource(std::string const& source); vtkNew ReverseLUT(vtkLookupTable* lut); -struct BandedGlyphs -{ - std::string glyphType; - vtkSmartPointer glyphs; - vtkSmartPointer bcf; - vtkSmartPointer lut; - vtkSmartPointer lutr; - vtkSmartPointer lut1; - vtkSmartPointer lut1r; - std::array scalarRange; - std::vector labels; -}; - -/** Get curvature glyphs and the corresponding banded polydata filter for the - * surface. - * - * @param surfaceName - The surface name. - * @param source - A vtkPolyData object corresponding to the source. - * @param curvature - The curvature type: Gauss_Curvature or Mean_Curvature. - * @param precision - the precision level. - * @param frequencyTable - If true, display a frequency table corresponding to - * the bands. - * @param nearestInteger - If true, use the nearest integer when generating the - * bands. - * @return A struct holding glyphs, bcf, lut, lutr, lut1, lut1r scalar_range, - * labels. - */ -BandedGlyphs -GetCurvatureGlyphs(std::string const& surfaceName, vtkPolyData* source, - std::string const& curvature, unsigned int precision, - bool frequencyTable = false, bool nearestInteger = false); - -/** - * Get elevation glyphs and the corresponding banded polydata filter for the - * surface. - * - * @param surfaceName - The surface name. - * @param source - A vtkPolyData object corresponding to the source. - * @param needsAdjusting - Surfaces whose curvatures need to be adjusted - * along the edges of the surface or constrained. - * @param frequencyTable - If true, display a frequency table corresponding to - * the bands. - * @param nearestInteger - If true, use the nearest integer when generating the - * bands. - * @return A struct holding glyphs, bcf, lut, lutr, lut1, lut1r scalar_range, - * labels. - */ -BandedGlyphs GetElevationGlyphs(std::string const& surfaceName, - vtkPolyData* source, unsigned int precision, - bool frequencyTable = false, - bool nearestInteger = false); - -/** - * Generate elevations over the surface. - * @param src - the vtkPolyData source. - * @return elev - the elevations. - */ -vtkSmartPointer GenerateElevations(vtkPolyData* src); - -/** - * Adjust curvatures along the edges of the surface. - * - * 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 curvatureName - The name of the curvature, "Gauss_Curvature" or - * "Mean_Curvature". - * @param epsilon - Curvature values less than this will be set to zero. - * @return - */ -void AdjustEdgeCurvatures(vtkPolyData* source, std::string const& curvatureName, - double const& epsilon = 1.0e-08); - -/** - * Constrain 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 curvatureName - The name of the curvature, "Gauss_Curvature" or - * "Mean_Curvature". - * @param lowerBound - The lower bound. - * @param upperBound - The upper bound. - * @return - */ -void ConstrainCurvatures(vtkPolyData* source, std::string const& curvatureName, - double const& lowerBound = 0.0, - double const& upperBound = 0.0); - /** * Glyph the normals on the surface. * @@ -239,6 +151,42 @@ GetIntegralBands(std::array const& scalarRange); std::map GetFrequencies(std::map>& bands, vtkPolyData* src); +/** + * Adjust curvatures along the edges of the surface. + * + * 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 curvatureName - The name of the curvature, "Gauss_Curvature" or + * "Mean_Curvature". + * @param epsilon - Curvature values less than this will be set to zero. + * @return + */ +void AdjustEdgeCurvatures(vtkPolyData* source, std::string const& curvatureName, + double const& epsilon = 1.0e-08); + +/** + * Constrain 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 curvatureName - The name of the curvature, "Gauss_Curvature" or + * "Mean_Curvature". + * @param lowerBound - The lower bound. + * @param upperBound - The upper bound. + * @return + */ +void ConstrainCurvatures(vtkPolyData* source, std::string const& curvatureName, + double const& lowerBound = 0.0, + double const& upperBound = 0.0); + /** * The bands and frequencies are adjusted so that the first and last * frequencies in the range are non-zero. @@ -246,18 +194,62 @@ std::map GetFrequencies(std::map>& bands, * @param bands: The bands. * @param freq: The frequencies. */ -void AdjustFrequencyRanges(std::map>& bands, - std::map& freq); +void AdjustRanges(std::map>& bands, + std::map& freq); + +struct BandedGlyphs +{ + std::string glyphType; + vtkSmartPointer glyphs; + vtkSmartPointer bcf; + int lutMaxBands; + vtkSmartPointer lut; + vtkSmartPointer lutr; + int lut1MaxBands; + vtkSmartPointer lut1; + vtkSmartPointer lut1r; + std::array scalarRange; + std::vector labels; +}; + +/** Get curvature glyphs and the corresponding banded polydata filter for the + * surface. + * + * @param surfaceName - The surface name. + * @param source - A vtkPolyData object corresponding to the source. + * @param curvature - The curvature type: Gauss_Curvature or Mean_Curvature. + * @param precision - the precision level. + * @param frequencyTable - If true, display a frequency table corresponding to + * the bands. + * @param nearestInteger - If true, use the nearest integer when generating the + * bands. + * @return A struct holding glyphs, bcf, lut, lutr, lut1, lut1r scalar_range, + * labels. + */ +BandedGlyphs +GetCurvatureGlyphs(std::string const& surfaceName, vtkPolyData* source, + std::string const& curvature, unsigned int precision, + bool frequencyTable = false, bool nearestInteger = false); /** - * Print each band and the number of scalars in each band. + * Get elevation glyphs and the corresponding banded polydata filter for the + * surface. * - * @param bands: The bands. - * @param freq: The frequencies. - * @param precision: The precision for the ranges in each band. + * @param surfaceName - The surface name. + * @param source - A vtkPolyData object corresponding to the source. + * @param needsAdjusting - Surfaces whose curvatures need to be adjusted + * along the edges of the surface or constrained. + * @param frequencyTable - If true, display a frequency table corresponding to + * the bands. + * @param nearestInteger - If true, use the nearest integer when generating the + * bands. + * @return A struct holding glyphs, bcf, lut, lutr, lut1, lut1r scalar_range, + * labels. */ -void PrintBandsFrequencies(std::map> const& bands, - std::map& freq, int const& precision = 2); +BandedGlyphs GetElevationGlyphs(std::string const& surfaceName, + vtkPolyData* source, unsigned int precision, + bool frequencyTable = false, + bool nearestInteger = false); typedef std::map> TTextPosition; typedef std::map TTextPositions; @@ -351,6 +343,16 @@ GetTextPositions(std::vector const& names, */ std::string Title(std::string s, const std::locale& loc = std::locale()); +/** + * Print each band and the number of scalars in each band. + * + * @param bands: The bands. + * @param freq: The frequencies. + * @param precision: The precision for the ranges in each band. + */ +void PrintBandsFrequencies(std::map> const& bands, + std::map& freq, int const& precision = 2); + /** * Adjust the camera parameters. * @@ -602,57 +604,37 @@ int main(int argc, char* argv[]) textWidget->On(); // Add scalar bars. - auto curvSBTtle = curvature; - std::replace(curvSBTtle.begin(), curvSBTtle.end(), '_', ' '); - auto curvSBP = ScalarBarProperties(); - curvSBP.titleText = curvSBTtle + "\n"; + curvSBP.titleText = curvBG.glyphType + "\n"; + std::replace(curvSBP.titleText.begin(), curvSBP.titleText.end(), '_', '\n'); curvSBP.numberOfLabels = static_cast(curvBG.labels.size()); // lut puts the lowest value at the top of the scalar bar. // lutr puts the highest value at the top of the scalar bar. curvSBP.lut = curvBG.lut; curvSBP.orientation = false; - auto maxBands = 9; - if (surfaceName == "hills") - { - curvSBP.positionH = std::move(PositionSBWH(9, maxBands)); - } - else if (surfaceName == "parametric hills") - { - curvSBP.positionH = std::move(PositionSBWH(7, maxBands)); - } - else if (surfaceName == "plane" || surfaceName == "sphere") + if (curvSBP.numberOfLabels == 1) { - curvSBP.positionH = std::move(PositionSBWH(1, maxBands)); + curvSBP.positionH = + PositionSBWH(curvSBP.numberOfLabels, curvBG.lutMaxBands + 3); } else { - curvSBP.positionH = std::move(PositionSBWH(5, maxBands)); + curvSBP.positionH = + PositionSBWH(curvSBP.numberOfLabels, curvBG.lutMaxBands); } auto cbw = MakeScalarBarWidget(curvSBP, titleTextProp, labelTextProp, ren, iren); auto elevSBP = ScalarBarProperties(); - elevSBP.titleText = "Elevation\n"; + elevSBP.titleText = elevBG.glyphType + "\n"; + std::replace(elevSBP.titleText.begin(), elevSBP.titleText.end(), '_', '\n'); elevSBP.numberOfLabels = static_cast(elevBG.labels.size()); // lut puts the lowest value at the top of the vertical scalar bar. // lutr puts the highest value at the top of the vertical scalar bar. elevSBP.lut = elevBG.lutr; elevSBP.orientation = true; - maxBands = 8; - if (surfaceName == "hills" || surfaceName == "parametric hills") - { - elevSBP.positionV = std::move(PositionSBWV(8, maxBands)); - } - else if (surfaceName == "plane") - { - elevSBP.positionV = std::move(PositionSBWV(1, maxBands)); - } - else - { - elevSBP.positionV = std::move(PositionSBWV(5, maxBands)); - } + elevSBP.positionV = PositionSBWV(elevSBP.numberOfLabels, elevBG.lut1MaxBands); auto sbw = MakeScalarBarWidget(elevSBP, titleTextProp, labelTextProp, ren, iren); @@ -692,1122 +674,1100 @@ int main(int argc, char* argv[]) } namespace { -void AdjustEdgeCurvatures(vtkPolyData* source, std::string const& curvatureName, - double const& epsilon) +vtkSmartPointer GenerateElevations(vtkPolyData* src) { - auto PointNeighbourhood = - [&source](vtkIdType const& pId) -> std::set { - // Extract the topological neighbors for point pId. In two steps: - // 1) source->GetPointCells(pId, cellIds) - // 2) source->GetCellPoints(cellId, cellPointIds) for all cellId in - // cellIds - vtkNew cellIds; - source->GetPointCells(pId, cellIds); - std::set neighbours; - for (vtkIdType i = 0; i < cellIds->GetNumberOfIds(); ++i) - { - auto cellId = cellIds->GetId(i); - vtkNew cellPointIds; - source->GetCellPoints(cellId, cellPointIds); - for (vtkIdType j = 0; j < cellPointIds->GetNumberOfIds(); ++j) - { - neighbours.insert(cellPointIds->GetId(j)); - } - } - return neighbours; - }; + 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[2] = 0.0; + bounds[3] = 1.0e-8; + } + vtkNew elevFilter; + elevFilter->SetInputData(src); + elevFilter->SetLowPoint(0, bounds[2], 0); + elevFilter->SetHighPoint(0, bounds[3], 0); + elevFilter->SetScalarRange(bounds[2], bounds[3]); + elevFilter->Update(); - auto ComputeDistance = [&source](vtkIdType const& ptIdA, - vtkIdType const& ptIdB) { - std::array ptA{0.0, 0.0, 0.0}; - std::array ptB{0.0, 0.0, 0.0}; - std::array ptC{0.0, 0.0, 0.0}; - source->GetPoint(ptIdA, ptA.data()); - source->GetPoint(ptIdB, ptB.data()); - std::transform(std::begin(ptA), std::end(ptA), std::begin(ptB), - std::begin(ptC), std::minus()); - // Calculate the norm. - auto result = std::sqrt(std::inner_product(std::begin(ptC), std::end(ptC), - std::begin(ptC), 0.0)); - return result; - }; - - source->GetPointData()->SetActiveScalars(curvatureName.c_str()); - // Curvature as a vector. - auto array = source->GetPointData()->GetAbstractArray(curvatureName.c_str()); - std::vector curvatures; - for (vtkIdType i = 0; i < source->GetNumberOfPoints(); ++i) - { - curvatures.push_back(array->GetVariantValue(i).ToDouble()); - } - - // Get the boundary point IDs. - std::string name = "Ids"; - vtkNew idFilter; - idFilter->SetInputData(source); - idFilter->SetPointIds(true); - idFilter->SetCellIds(false); - idFilter->SetPointIdsArrayName(name.c_str()); - idFilter->SetCellIdsArrayName(name.c_str()); - idFilter->Update(); + return elevFilter->GetPolyDataOutput(); +} - vtkNew edges; +vtkSmartPointer GetHills() +{ + // Create four hills on a plane. + // This will have regions of negative, zero and positive Gsaussian curvatures. - edges->SetInputConnection(idFilter->GetOutputPort()); - edges->BoundaryEdgesOn(); - edges->ManifoldEdgesOff(); - edges->NonManifoldEdgesOff(); - edges->FeatureEdgesOff(); - edges->Update(); + 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); - auto edgeAarray = - edges->GetOutput()->GetPointData()->GetAbstractArray(name.c_str()); - std::vector boundaryIds; - for (vtkIdType i = 0; i < edges->GetOutput()->GetNumberOfPoints(); ++i) - { - boundaryIds.push_back(edgeAarray->GetVariantValue(i).ToInt()); - } - // Remove duplicate Ids. - std::set pIdsSet(boundaryIds.begin(), boundaryIds.end()); - for (auto const pId : boundaryIds) + // Make a grid. + vtkNew points; + for (auto i = 0; i < xRes; ++i) { - auto pIdsNeighbors = PointNeighbourhood(pId); - std::set pIdsNeighborsInterior; - std::set_difference( - pIdsNeighbors.begin(), pIdsNeighbors.end(), pIdsSet.begin(), - pIdsSet.end(), - std::inserter(pIdsNeighborsInterior, pIdsNeighborsInterior.begin())); - // Compute distances and extract curvature values. - std::vector curvs; - std::vector dists; - for (auto const pIdN : pIdsNeighborsInterior) - { - curvs.push_back(curvatures[pIdN]); - dists.push_back(ComputeDistance(pIdN, pId)); - } - std::vector nonZeroDistIds; - for (size_t i = 0; i < dists.size(); ++i) - { - if (dists[i] > 0) - { - nonZeroDistIds.push_back(i); - } - } - std::vector curvsNonZero; - std::vector distsNonZero; - for (auto const i : nonZeroDistIds) - { - curvsNonZero.push_back(curvs[i]); - distsNonZero.push_back(dists[i]); - } - // Iterate over the edge points and compute the curvature as the weighted - // average of the neighbours. - auto countInvalid = 0; - auto newCurv = 0.0; - if (curvsNonZero.size() > 0) - { - std::vector weights; - double sum = 0.0; - for (auto const d : distsNonZero) - { - sum += 1.0 / d; - weights.push_back(1.0 / d); - } - for (size_t i = 0; i < weights.size(); ++i) - { - weights[i] = weights[i] / sum; - } - newCurv = std::inner_product(curvsNonZero.begin(), curvsNonZero.end(), - weights.begin(), 0.0); - } - else + auto x = xMin + i * dx; + for (auto j = 0; j < yRes; ++j) { - // Corner case. - // countInvalid += 1; - // Assuming the curvature of the point is planar. - newCurv = 0.0; + auto y = yMin + j * dy; + points->InsertNextPoint(x, y, 0); } - // Set the new curvature value. - curvatures[pId] = newCurv; } - // Set small values to zero. - if (epsilon != 0.0) - { - auto eps = std::abs(epsilon); - for (size_t i = 0; i < curvatures.size(); ++i) - { - if (std::abs(curvatures[i]) < eps) - { - curvatures[i] = 0.0; - } - } - } + // Add the grid points to a polydata object. + vtkNew plane; + plane->SetPoints(points); - if (static_cast(source->GetNumberOfPoints()) != curvatures.size()) - { - std::string s = curvatureName; - s += ":\nCannot add the adjusted curvatures to the source.\n"; - s += " The number of points in source does not equal the\n"; - s += " number of point ids in the adjusted curvature array."; - std::cerr << s << std::endl; - return; - } - vtkNew adjustedCurvatures; - adjustedCurvatures->SetName(curvatureName.c_str()); - for (auto curvature : curvatures) - { - adjustedCurvatures->InsertNextTuple1(curvature); - } - source->GetPointData()->AddArray(adjustedCurvatures); - source->GetPointData()->SetActiveScalars(curvatureName.c_str()); -} + // Triangulate the grid. + vtkNew delaunay; + delaunay->SetInputData(plane); + delaunay->Update(); -void ConstrainCurvatures(vtkPolyData* source, std::string const& curvatureName, - double const& lowerBound, double const& upperBound) -{ - std::array bounds{0.0, 0.0}; - if (lowerBound < upperBound) - { - bounds[0] = lowerBound; - bounds[1] = upperBound; - } - else - { - bounds[0] = upperBound; - bounds[1] = lowerBound; - } + auto polydata = delaunay->GetOutput(); - source->GetPointData()->SetActiveScalars(curvatureName.c_str()); - // Curvature as a vector. - auto array = source->GetPointData()->GetAbstractArray(curvatureName.c_str()); - std::vector curvatures; - for (vtkIdType i = 0; i < source->GetNumberOfPoints(); ++i) - { - curvatures.push_back(array->GetVariantValue(i).ToDouble()); - } - // Set upper and lower bounds. - for (size_t i = 0; i < curvatures.size(); ++i) + vtkNew 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> 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}}; + // Three hills. + // std::vector> hd{{-2.5, -2.5, 2.5, 6.5, 3.5}, + // {5.0, -2.5, 1.5, 1.5, 2.5}, + // {-5.0, 5, 2.5, 3.0, 3}}; + // Two hills. + // std::vector> hd{{5.0, -2.5, 1.5, 1.5, 2.5}, + // {-5.0, 5, 2.5, 3.0, 3}}; + + std::array xx{0.0, 0.0}; + for (auto i = 0; i < points->GetNumberOfPoints(); ++i) { - if (curvatures[i] < bounds[0]) - { - curvatures[i] = bounds[0]; - } - else + auto x = polydata->GetPoint(i); + for (size_t j = 0; j < hd.size(); ++j) { - if (curvatures[i] > bounds[1]) - { - curvatures[i] = bounds[1]; - } + 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); } - } - vtkNew adjustedCurvatures; - for (auto curvature : curvatures) - { - adjustedCurvatures->InsertNextTuple1(curvature); - } - adjustedCurvatures->SetName(curvatureName.c_str()); - source->GetPointData()->RemoveArray(curvatureName.c_str()); - source->GetPointData()->AddArray(adjustedCurvatures); - source->GetPointData()->SetActiveScalars(curvatureName.c_str()); -} - -BandedGlyphs GetCurvatureGlyphs(std::string const& surfaceName, - vtkPolyData* source, - std::string const& curvature, - unsigned int precision, bool frequencyTable, - bool nearestInteger) -{ - // The length of the normal arrow glyphs. - auto scaleFactor = 1.0; - if (surfaceName == "hills") - { - scaleFactor = 0.5; - } - if (surfaceName == "sphere") - { - scaleFactor = 0.2; + polydata->GetPoints()->SetPoint(i, x); + elevation->SetValue(i, x[2]); } - source->GetPointData()->SetActiveScalars(curvature.c_str()); - auto sr = source->GetPointData()->GetScalars(curvature.c_str())->GetRange(); - std::array scalarRange{sr[0], sr[1]}; + vtkNew textures; + textures->SetNumberOfComponents(2); + textures->SetNumberOfTuples(2 * polydata->GetNumberOfPoints()); + textures->SetName("Textures"); - vtkNew colorSeries; - if (surfaceName == "hills") - { - colorSeries->SetColorScheme(colorSeries->BREWER_DIVERGING_SPECTRAL_7); - } - else if (surfaceName == "parametric hills") - { - colorSeries->SetColorScheme(colorSeries->BREWER_DIVERGING_SPECTRAL_9); - } - else if (surfaceName == "plane" || surfaceName == "sphere") - { - colorSeries->SetColorScheme(colorSeries->CITRUS); - } - else + for (auto i = 0; i < xRes; ++i) { - colorSeries->SetColorScheme(colorSeries->BREWER_DIVERGING_SPECTRAL_5); + float tc[2]; + tc[0] = i / (xRes - 1.0); + for (auto j = 0; j < yRes; ++j) + { + // tc[1] = 1.0 - j / (yRes - 1.0); + tc[1] = j / (yRes - 1.0); + textures->SetTuple(static_cast(i) * yRes + j, tc); + } } - vtkNew lut; - colorSeries->BuildLookupTable(lut, colorSeries->CATEGORICAL); - lut->SetNanColor(0, 0, 0, 1); - lut->SetTableRange(scalarRange.data()); + polydata->GetPointData()->SetScalars(elevation); + polydata->GetPointData()->GetScalars()->SetName("Elevation"); + polydata->GetPointData()->SetTCoords(textures); - vtkNew lut1; - colorSeries->BuildLookupTable(lut1, colorSeries->ORDINAL); - lut1->SetNanColor(0, 0, 0, 1); - lut1->SetTableRange(scalarRange.data()); + vtkNew normals; + normals->SetInputData(polydata); + normals->SetInputData(polydata); + normals->SetFeatureAngle(30); + normals->SplittingOff(); - auto numberOfBands = lut->GetNumberOfTableValues(); - lut1->SetNumberOfTableValues(numberOfBands); + vtkNew tr1; + tr1->RotateX(-90); - auto bands = GetBands(scalarRange, numberOfBands, precision, nearestInteger); + vtkNew tf1; + tf1->SetInputConnection(normals->GetOutputPort()); + tf1->SetTransform(tr1); + tf1->Update(); - if (curvature == "Gaussian_Curvature") - { - if (surfaceName == "random hills") - { - // These are my custom bands. - // Generated by first running: - // bands = GetBands(curvScalarRanges, curvNumBands, precision, false); - // then: - // std::vector freq = Frequencies(bands, src); - // PrintBandsFrequencies(bands, freq); - // Finally using the output to create this table: - // std::vector> 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},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) - // and some spherical regions. - std::vector> myBands = { - {-0.630, -0.190}, {-0.190, -0.043}, {-0.043, 0.0452}, - {0.0452, 0.0746}, {0.0746, 0.104}, {0.104, 0.251}, - {0.251, 1.131}}; + return tf1->GetPolyDataOutput(); +} - // Comment this out if you want to see how allocating - // equally spaced bands works. - bands = GetCustomBands(scalarRange, numberOfBands, myBands); - } - else if (surfaceName == "hills") - { - std::vector> myBands = { - {-2.104, -0.15}, {-0.15, -0.1}, {-0.1, -0.05}, - {-0.05, -0.02}, {-0.02, -0.005}, {-0.005, -0.0005}, - {-0.0005, 0.0005}, {0.0005, 0.09}, {0.09, 4.972}, - }; +vtkSmartPointer GetParametricHills() +{ + vtkNew fn; + fn->AllowRandomGenerationOn(); + fn->SetRandomSeed(1); + fn->SetNumberOfHills(30); - // Comment this out if you want to see how allocating - // equally spaced bands works. - bands = GetCustomBands(scalarRange, numberOfBands, myBands); - } - } + vtkNew source; + source->SetParametricFunction(fn); + source->SetUResolution(50); + source->SetVResolution(50); + source->SetScalarModeToZ(); + source->Update(); - // Adjust the number of table values and scalar range. - scalarRange[0] = bands.begin()->second[0]; - scalarRange[1] = std::prev(bands.end())->second[2]; - lut->SetTableRange(scalarRange.data()); - lut->SetNumberOfTableValues(bands.size()); - lut1->SetTableRange(scalarRange.data()); - lut1->SetNumberOfTableValues(bands.size()); + // Name the arrays (not needed in VTK 6.2+ for vtkParametricFunctionSource). + // source->GetOutput()->GetPointData()->GetNormals()->SetName("Normals"); + // source->GetOutput()->GetPointData()->GetScalars()->SetName("Scalars"); + // Rename the scalars to "Elevation" since we are using the Z-scalars as + // elevations. + source->GetOutput()->GetPointData()->GetScalars()->SetName("Elevation"); - if (frequencyTable) - { - auto freq = GetFrequencies(bands, source); - AdjustFrequencyRanges(bands, freq); - fmt::println("{:s} {:s}", Title(surfaceName), curvature); - PrintBandsFrequencies(bands, freq); - } + vtkNew transform; + transform->Translate(0.0, 5.0, 15.0); + transform->RotateX(-90.0); + vtkNew transformFilter; + transformFilter->SetInputConnection(source->GetOutputPort()); + transformFilter->SetTransform(transform); + transformFilter->Update(); - // We will use the midpoint of the band as the label. - std::vector labels; - for (std::map>::const_iterator p = bands.begin(); - p != bands.end(); ++p) - { - labels.push_back(fmt::format("{:4.2f}", p->second[1])); - } + return transformFilter->GetPolyDataOutput(); +} - // Annotate - vtkNew values; - for (size_t i = 0; i < labels.size(); ++i) - { - values->InsertNextValue(vtkVariant(labels[i])); - } - for (vtkIdType i = 0; i < values->GetNumberOfTuples(); ++i) - { - lut->SetAnnotation(i, values->GetValue(i).ToString()); - } +vtkSmartPointer GetParametricTorus() +{ + vtkNew fn; + fn->SetRingRadius(5); + fn->SetCrossSectionRadius(2); - // Create the contour bands. - vtkNew bcf; - bcf->SetInputData(source); - // Use either the minimum or maximum value for each band. - int i = 0; - for (std::map>::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(); + vtkNew source; + source->SetParametricFunction(fn); + source->SetUResolution(50); + source->SetVResolution(50); + source->SetScalarModeToZ(); + source->Update(); - // Generate the glyphs on the original surface. - auto glyphs = GetGlyphs(source, scaleFactor, false); + // Name the arrays (not needed in VTK 6.2+ for vtkParametricFunctionSource). + // source->GetOutput()->GetPointData()->GetNormals()->SetName("Normals"); + // source->GetOutput()->GetPointData()->GetScalars()->SetName("Scalars"); + // Rename the scalars to "Elevation" since we are using the Z-scalars as + // elevations. + source->GetOutput()->GetPointData()->GetScalars()->SetName("Elevation"); - BandedGlyphs bg; - bg.glyphType = curvature; - bg.glyphs = glyphs; - bg.bcf = bcf; - bg.lut = lut; - bg.lutr = ReverseLUT(lut); - bg.lut1 = lut1; - bg.lut1r = ReverseLUT(lut1); - bg.scalarRange = scalarRange; - bg.labels = labels; + vtkNew transform; + transform->RotateX(-90.0); + vtkNew transformFilter; + transformFilter->SetInputConnection(source->GetOutputPort()); + transformFilter->SetTransform(transform); + transformFilter->Update(); - return bg; + return transformFilter->GetPolyDataOutput(); } -BandedGlyphs GetElevationGlyphs(std::string const& surfaceName, - vtkPolyData* source, unsigned int precision, - bool frequencyTable, bool nearestInteger) +vtkSmartPointer GetPlane() { - // The length of the normal arrow glyphs. - auto scaleFactor = 1.0; - if (surfaceName == "hills") - { - scaleFactor = 0.5; - } - if (surfaceName == "sphere") - { - scaleFactor = 0.2; - } + vtkNew source; + source->SetOrigin(-10.0, -10.0, 0.0); + source->SetPoint2(-10.0, 10.0, 0.0); + source->SetPoint1(10.0, -10.0, 0.0); + source->SetXResolution(10); + source->SetYResolution(10); + source->Update(); - source->GetPointData()->SetActiveScalars("Elevation"); - auto sr = source->GetPointData()->GetScalars("Elevation")->GetRange(); - std::array scalarRange{sr[0], sr[1]}; + vtkNew transform; + transform->Translate(0.0, 0.0, 0.0); + transform->RotateX(-90.0); + vtkNew transformFilter; + transformFilter->SetInputConnection(source->GetOutputPort()); + transformFilter->SetTransform(transform); + transformFilter->Update(); - vtkNew colorSeries; - if (surfaceName == "hills" || surfaceName == "parametric hills") - { - colorSeries->SetColorScheme( - colorSeries->BREWER_DIVERGING_BROWN_BLUE_GREEN_8); - } - else - { - colorSeries->SetColorScheme( - colorSeries->BREWER_DIVERGING_BROWN_BLUE_GREEN_5); - } + // 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. + vtkNew tri; + tri->SetInputConnection(transformFilter->GetOutputPort()); - vtkNew lut; - colorSeries->BuildLookupTable(lut, colorSeries->CATEGORICAL); - lut->SetNanColor(0, 0, 0, 1); - lut->SetTableRange(scalarRange.data()); + // Pass it though a CleanPolyDataFilter and merge any points which + // are coincident, or very close + vtkNew cleaner; + cleaner->SetInputConnection(tri->GetOutputPort()); + cleaner->SetTolerance(0.005); + cleaner->Update(); + + return cleaner->GetOutput(); +} + +vtkSmartPointer GetSphere() +{ + vtkNew source; + source->SetCenter(0.0, 0.0, 0.0); + source->SetRadius(1.0); + source->SetThetaResolution(32); + source->SetPhiResolution(32); + source->Update(); + + return source->GetOutput(); +} + +vtkSmartPointer GetTorus() +{ + vtkNew source; + source->SetCenter(0.0, 0.0, 0.0); + source->SetPhiResolution(64); + source->SetThetaResolution(64); + source->SetThetaRoundness(1); + source->SetThickness(0.5); + source->SetSize(10); + source->SetToroidal(1); + source->Update(); - vtkNew lut1; - colorSeries->BuildLookupTable(lut1, colorSeries->ORDINAL); - lut1->SetNanColor(0, 0, 0, 1); - lut1->SetTableRange(scalarRange.data()); + // The quadric is made of strips, so pass it through a triangle filter as + // the curvature filter only operates on polys + vtkNew tri; + tri->SetInputConnection(source->GetOutputPort()); - auto numberOfBands = lut->GetNumberOfTableValues(); - lut1->SetNumberOfTableValues(numberOfBands); + // 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 + vtkNew cleaner; + cleaner->SetInputConnection(tri->GetOutputPort()); + cleaner->SetTolerance(0.005); + cleaner->Update(); - auto bands = GetBands(scalarRange, numberOfBands, precision, nearestInteger); + return cleaner->GetOutput(); +} - if (surfaceName == "parametric hills") +vtkSmartPointer GetSource(std::string const& source) +{ + std::string surface = source; + std::transform(surface.begin(), surface.end(), surface.begin(), + [](unsigned char c) { return std::tolower(c); }); + if (surface == "hills") { - // These are my custom bands. - // Generated by first running: - // bands = GetBands(scalarRange, numberOfBands, precision, false); - // then: - // std::vector freq = Frequencies(bands, source); - // PrintBandsFrequencies(bands, freq); - // Finally using the output to create this table: - std::vector> myBands = { - {0, 1.0}, {1.0, 2.0}, {2.0, 3.0}, {3.0, 4.0}, - {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 = GetCustomBands(scalarRange, numberOfBands, myBands); + return GetHills(); } - - // Adjust the number of table values and scalar range. - scalarRange[0] = bands.begin()->second[0]; - scalarRange[1] = std::prev(bands.end())->second[2]; - lut->SetTableRange(scalarRange.data()); - lut->SetNumberOfTableValues(bands.size()); - lut1->SetTableRange(scalarRange.data()); - lut1->SetNumberOfTableValues(bands.size()); - - if (frequencyTable) + else if (surface == "parametric hills") { - auto freq = GetFrequencies(bands, source); - AdjustFrequencyRanges(bands, freq); - fmt::println("{:s} Elevation", Title(surfaceName)); - PrintBandsFrequencies(bands, freq); + return GetParametricHills(); } - - // We will use the midpoint of the band as the label. - std::vector labels; - for (std::map>::const_iterator p = bands.begin(); - p != bands.end(); ++p) + else if (surface == "parametric torus") { - labels.push_back(fmt::format("{:4.2f}", p->second[1])); + return GetParametricTorus(); } - - // Annotate - vtkNew values; - for (size_t i = 0; i < labels.size(); ++i) + else if (surface == "plane") { - values->InsertNextValue(vtkVariant(labels[i])); + return GenerateElevations(GetPlane()); } - for (vtkIdType i = 0; i < values->GetNumberOfTuples(); ++i) + else if (surface == "sphere") { - lut->SetAnnotation(i, values->GetValue(i).ToString()); + return GenerateElevations(GetSphere()); } - - // Create the contour bands. - vtkNew bcf; - bcf->SetInputData(source); - // Use either the minimum or maximum value for each band. - int i = 0; - for (std::map>::const_iterator p = bands.begin(); - p != bands.end(); ++p) + else if (surface == "torus") { - bcf->SetValue(i, p->second[2]); - ++i; + return GenerateElevations(GetTorus()); } - // We will use an indexed lookup table. - bcf->SetScalarModeToIndex(); - bcf->GenerateContourEdgesOn(); - - // Generate the glyphs on the original surface. - auto glyphs = GetGlyphs(source, scaleFactor, false); + std::cout << "The surface is not available." << std::endl; + std::cout << "Using parametric hills instead." << std::endl; + return GetParametricHills(); +} - BandedGlyphs bg; - bg.glyphType = "Elevation"; - bg.glyphs = glyphs; - bg.bcf = bcf; - bg.lut = lut; - bg.lutr = ReverseLUT(lut); - bg.lut1 = lut1; - bg.lut1r = ReverseLUT(lut1); - bg.scalarRange = scalarRange; - bg.labels = labels; +vtkNew ReverseLUT(vtkLookupTable* lut) +{ + // First do a deep copy just to get the whole structure + // and then reverse the colors and annotations. + vtkNew lutr; + lutr->DeepCopy(lut); + vtkIdType t = lut->GetNumberOfTableValues() - 1; + for (vtkIdType i = t; i >= 0; --i) + { + std::array rgb{0.0, 0.0, 0.0}; + std::array rgba{0.0, 0.0, 0.0, 1.0}; + lut->GetColor(i, rgb.data()); + std::copy(std::begin(rgb), std::end(rgb), std::begin(rgba)); + rgba[3] = lut->GetOpacity(i); + lutr->SetTableValue(t - i, rgba.data()); + } + t = lut->GetNumberOfAnnotatedValues() - 1; + for (vtkIdType i = t; i >= 0; --i) + { + lutr->SetAnnotation(t - i, lut->GetAnnotation(i)); + } - return bg; + return lutr; } -vtkSmartPointer GenerateElevations(vtkPolyData* src) +vtkNew GetGlyphs(vtkPolyData* src, double const& scaleFactor, + bool const& reverseNormals) { - 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) + // 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 reverse; + vtkNew maskPts; + maskPts->SetOnRatio(5); + maskPts->RandomModeOn(); + if (reverseNormals) { - bounds[2] = 0.0; - bounds[3] = 1.0e-8; + reverse->SetInputData(src); + reverse->ReverseCellsOn(); + reverse->ReverseNormalsOn(); + maskPts->SetInputConnection(reverse->GetOutputPort()); + } + else + { + maskPts->SetInputData(src); } - vtkNew 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(); + // Source for the glyph filter + vtkNew arrow; + arrow->SetTipResolution(16); + arrow->SetTipLength(0.3); + arrow->SetTipRadius(0.1); + + vtkNew glyph; + glyph->SetSourceConnection(arrow->GetOutputPort()); + glyph->SetInputConnection(maskPts->GetOutputPort()); + glyph->SetVectorModeToUseNormal(); + glyph->SetScaleFactor(scaleFactor); + glyph->SetColorModeToColorByVector(); + glyph->SetScaleModeToScaleByVector(); + glyph->OrientOn(); + glyph->Update(); + + return glyph; } -vtkSmartPointer GetHills() +std::map> +GetBands(std::array const& scalarRange, int const& numberOfBands, + int const& precision, bool const& nearestInteger) { - // Create four hills on a plane. - // This will have regions of negative, zero and positive Gsaussian curvatures. + auto prec = abs(precision); + prec = (prec > 14) ? 14 : prec; - 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); + auto RoundOff = [&prec](const double& x) { + auto pow_10 = std::pow(10.0, prec); + return std::round(x * pow_10) / pow_10; + }; - // Make a grid. - vtkNew points; - for (auto i = 0; i < xRes; ++i) + std::map> bands; + if ((scalarRange[1] < scalarRange[0]) || (numberOfBands <= 0)) { - auto x = xMin + i * dx; - for (auto j = 0; j < yRes; ++j) + return bands; + } + double x[2]; + for (int i = 0; i < 2; ++i) + { + x[i] = scalarRange[i]; + } + if (nearestInteger) + { + x[0] = std::floor(x[0]); + x[1] = std::ceil(x[1]); + } + double dx = (x[1] - x[0]) / static_cast(numberOfBands); + std::vector 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) { - auto y = yMin + j * dy; - points->InsertNextPoint(x, y, 0); + for (std::vector::iterator p = b.begin(); p != b.end(); ++p) + { + *p = RoundOff(*p); + } + b[0] = x[0]; + } + bands[i] = b; + for (std::vector::iterator p = b.begin(); p != b.end(); ++p) + { + *p = RoundOff(*p + dx); } } + return bands; +} - // Add the grid points to a polydata object. - vtkNew plane; - plane->SetPoints(points); - - // Triangulate the grid. - vtkNew delaunay; - delaunay->SetInputData(plane); - delaunay->Update(); - - auto polydata = delaunay->GetOutput(); - - vtkNew elevation; - elevation->SetNumberOfTuples(points->GetNumberOfPoints()); +std::map> +GetCustomBands(std::array const& scalarRange, + int const& numberOfBands, + std::vector> const& myBands) +{ + std::map> bands; + if ((scalarRange[1] < scalarRange[0]) || (numberOfBands <= 0)) + { + return bands; + } - // We define the parameters for the hills here. - // [[0: x0, 1: y0, 2: x variance, 3: y variance, 4: amplitude]...] - std::vector> 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}}; - // Three hills. - // std::vector> hd{{-2.5, -2.5, 2.5, 6.5, 3.5}, - // {5.0, -2.5, 1.5, 1.5, 2.5}, - // {-5.0, 5, 2.5, 3.0, 3}}; - // Two hills. - // std::vector> hd{{5.0, -2.5, 1.5, 1.5, 2.5}, - // {-5.0, 5, 2.5, 3.0, 3}}; + std::vector> x; + std::copy(myBands.begin(), myBands.end(), std::back_inserter(x)); - std::array xx{0.0, 0.0}; - for (auto i = 0; i < points->GetNumberOfPoints(); ++i) + // Determine the index of the range minimum and range maximum. + int idxMin = 0; + for (auto idx = 0; idx < static_cast(myBands.size()); ++idx) { - auto x = polydata->GetPoint(i); - for (size_t j = 0; j < hd.size(); ++j) + if (scalarRange[0] < myBands[idx][1] && scalarRange[0] >= myBands[idx][0]) { - 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); + idxMin = idx; + break; } - polydata->GetPoints()->SetPoint(i, x); - elevation->SetValue(i, x[2]); } - - vtkNew textures; - textures->SetNumberOfComponents(2); - textures->SetNumberOfTuples(2 * polydata->GetNumberOfPoints()); - textures->SetName("Textures"); - - for (auto i = 0; i < xRes; ++i) + int idxMax = static_cast(myBands.size()) - 1; + for (int idx = static_cast(myBands.size()) - 1; idx >= 0; --idx) { - float tc[2]; - tc[0] = i / (xRes - 1.0); - for (auto j = 0; j < yRes; ++j) + if (scalarRange[1] < myBands[idx][1] && scalarRange[1] >= myBands[idx][0]) { - // tc[1] = 1.0 - j / (yRes - 1.0); - tc[1] = j / (yRes - 1.0); - textures->SetTuple(static_cast(i) * yRes + j, tc); + idxMax = static_cast(idx); + break; } } - polydata->GetPointData()->SetScalars(elevation); - polydata->GetPointData()->GetScalars()->SetName("Elevation"); - polydata->GetPointData()->SetTCoords(textures); - - vtkNew normals; - normals->SetInputData(polydata); - normals->SetInputData(polydata); - normals->SetFeatureAngle(30); - normals->SplittingOff(); - - vtkNew tr1; - tr1->RotateX(-90); - - vtkNew tf1; - tf1->SetInputConnection(normals->GetOutputPort()); - tf1->SetTransform(tr1); - tf1->Update(); - - return tf1->GetPolyDataOutput(); + // Set the minimum to match the range minimum. + x[idxMin][0] = scalarRange[0]; + x[idxMax][1] = scalarRange[1]; + for (int i = idxMin; i < idxMax + 1; ++i) + { + std::vector 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; } -vtkSmartPointer GetParametricHills() +std::map> +GetIntegralBands(std::array const& scalarRange) { - vtkNew fn; - fn->AllowRandomGenerationOn(); - fn->SetRandomSeed(1); - fn->SetNumberOfHills(30); - - vtkNew source; - source->SetParametricFunction(fn); - source->SetUResolution(50); - source->SetVResolution(50); - source->SetScalarModeToZ(); - source->Update(); - - // Name the arrays (not needed in VTK 6.2+ for vtkParametricFunctionSource). - // source->GetOutput()->GetPointData()->GetNormals()->SetName("Normals"); - // source->GetOutput()->GetPointData()->GetScalars()->SetName("Scalars"); - // Rename the scalars to "Elevation" since we are using the Z-scalars as - // elevations. - source->GetOutput()->GetPointData()->GetScalars()->SetName("Elevation"); - - vtkNew transform; - transform->Translate(0.0, 5.0, 15.0); - transform->RotateX(-90.0); - vtkNew transformFilter; - transformFilter->SetInputConnection(source->GetOutputPort()); - transformFilter->SetTransform(transform); - transformFilter->Update(); - - return transformFilter->GetPolyDataOutput(); + if (scalarRange[1] < scalarRange[0]) + { + std::map> bands; + return bands; + } + std::array x{0, 0}; + for (int i = 0; i < 2; ++i) + { + x[i] = scalarRange[i]; + } + x[0] = std::floor(x[0]); + x[1] = std::ceil(x[1]); + int numberOfBands = static_cast(std::abs(x[1]) - std::abs(x[0])); + return GetBands(x, numberOfBands, false); } -vtkSmartPointer GetParametricTorus() +std::map GetFrequencies(std::map>& bands, + vtkPolyData* src) { - vtkNew fn; - fn->SetRingRadius(5); - fn->SetCrossSectionRadius(2); - - vtkNew source; - source->SetParametricFunction(fn); - source->SetUResolution(50); - source->SetVResolution(50); - source->SetScalarModeToZ(); - source->Update(); - - // Name the arrays (not needed in VTK 6.2+ for vtkParametricFunctionSource). - // source->GetOutput()->GetPointData()->GetNormals()->SetName("Normals"); - // source->GetOutput()->GetPointData()->GetScalars()->SetName("Scalars"); - // Rename the scalars to "Elevation" since we are using the Z-scalars as - // elevations. - source->GetOutput()->GetPointData()->GetScalars()->SetName("Elevation"); - - vtkNew transform; - transform->RotateX(-90.0); - vtkNew transformFilter; - transformFilter->SetInputConnection(source->GetOutputPort()); - transformFilter->SetTransform(transform); - transformFilter->Update(); - - return transformFilter->GetPolyDataOutput(); + std::map freq; + for (auto i = 0; i < static_cast(bands.size()); ++i) + { + freq[i] = 0; + } + vtkIdType tuples = src->GetPointData()->GetScalars()->GetNumberOfTuples(); + for (int i = 0; i < tuples; ++i) + { + const double* x = src->GetPointData()->GetScalars()->GetTuple(i); + for (auto j = 0; j < static_cast(bands.size()); ++j) + { + if (*x <= bands[j][2]) + { + freq[j] = freq[j] + 1; + break; + } + } + } + return freq; } -vtkSmartPointer GetPlane() +void AdjustEdgeCurvatures(vtkPolyData* source, std::string const& curvatureName, + double const& epsilon) { - vtkNew source; - source->SetOrigin(-10.0, -10.0, 0.0); - source->SetPoint2(-10.0, 10.0, 0.0); - source->SetPoint1(10.0, -10.0, 0.0); - source->SetXResolution(10); - source->SetYResolution(10); - source->Update(); - - vtkNew transform; - transform->Translate(0.0, 0.0, 0.0); - transform->RotateX(-90.0); - vtkNew transformFilter; - transformFilter->SetInputConnection(source->GetOutputPort()); - transformFilter->SetTransform(transform); - transformFilter->Update(); - - // 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. - vtkNew tri; - tri->SetInputConnection(transformFilter->GetOutputPort()); - - // Pass it though a CleanPolyDataFilter and merge any points which - // are coincident, or very close - vtkNew cleaner; - cleaner->SetInputConnection(tri->GetOutputPort()); - cleaner->SetTolerance(0.005); - cleaner->Update(); + auto PointNeighbourhood = + [&source](vtkIdType const& pId) -> std::set { + // Extract the topological neighbors for point pId. In two steps: + // 1) source->GetPointCells(pId, cellIds) + // 2) source->GetCellPoints(cellId, cellPointIds) for all cellId in + // cellIds + vtkNew cellIds; + source->GetPointCells(pId, cellIds); + std::set neighbours; + for (vtkIdType i = 0; i < cellIds->GetNumberOfIds(); ++i) + { + auto cellId = cellIds->GetId(i); + vtkNew cellPointIds; + source->GetCellPoints(cellId, cellPointIds); + for (vtkIdType j = 0; j < cellPointIds->GetNumberOfIds(); ++j) + { + neighbours.insert(cellPointIds->GetId(j)); + } + } + return neighbours; + }; - return cleaner->GetOutput(); -} + auto ComputeDistance = [&source](vtkIdType const& ptIdA, + vtkIdType const& ptIdB) { + std::array ptA{0.0, 0.0, 0.0}; + std::array ptB{0.0, 0.0, 0.0}; + std::array ptC{0.0, 0.0, 0.0}; + source->GetPoint(ptIdA, ptA.data()); + source->GetPoint(ptIdB, ptB.data()); + std::transform(std::begin(ptA), std::end(ptA), std::begin(ptB), + std::begin(ptC), std::minus()); + // Calculate the norm. + auto result = std::sqrt(std::inner_product(std::begin(ptC), std::end(ptC), + std::begin(ptC), 0.0)); + return result; + }; -vtkSmartPointer GetSphere() -{ - vtkNew source; - source->SetCenter(0.0, 0.0, 0.0); - source->SetRadius(1.0); - source->SetThetaResolution(32); - source->SetPhiResolution(32); - source->Update(); + source->GetPointData()->SetActiveScalars(curvatureName.c_str()); + // Curvature as a vector. + auto array = source->GetPointData()->GetAbstractArray(curvatureName.c_str()); + std::vector curvatures; + for (vtkIdType i = 0; i < source->GetNumberOfPoints(); ++i) + { + curvatures.push_back(array->GetVariantValue(i).ToDouble()); + } - return source->GetOutput(); -} + // Get the boundary point IDs. + std::string name = "Ids"; + vtkNew idFilter; + idFilter->SetInputData(source); + idFilter->SetPointIds(true); + idFilter->SetCellIds(false); + idFilter->SetPointIdsArrayName(name.c_str()); + idFilter->SetCellIdsArrayName(name.c_str()); + idFilter->Update(); -vtkSmartPointer GetTorus() -{ - vtkNew source; - source->SetCenter(0.0, 0.0, 0.0); - source->SetPhiResolution(64); - source->SetThetaResolution(64); - source->SetThetaRoundness(1); - source->SetThickness(0.5); - source->SetSize(10); - source->SetToroidal(1); - source->Update(); + vtkNew edges; - // The quadric is made of strips, so pass it through a triangle filter as - // the curvature filter only operates on polys - vtkNew tri; - tri->SetInputConnection(source->GetOutputPort()); + edges->SetInputConnection(idFilter->GetOutputPort()); + edges->BoundaryEdgesOn(); + edges->ManifoldEdgesOff(); + edges->NonManifoldEdgesOff(); + edges->FeatureEdgesOff(); + edges->Update(); - // 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 - vtkNew cleaner; - cleaner->SetInputConnection(tri->GetOutputPort()); - cleaner->SetTolerance(0.005); - cleaner->Update(); + auto edgeAarray = + edges->GetOutput()->GetPointData()->GetAbstractArray(name.c_str()); + std::vector boundaryIds; + for (vtkIdType i = 0; i < edges->GetOutput()->GetNumberOfPoints(); ++i) + { + boundaryIds.push_back(edgeAarray->GetVariantValue(i).ToInt()); + } + // Remove duplicate Ids. + std::set pIdsSet(boundaryIds.begin(), boundaryIds.end()); + for (auto const pId : boundaryIds) + { + auto pIdsNeighbors = PointNeighbourhood(pId); + std::set pIdsNeighborsInterior; + std::set_difference( + pIdsNeighbors.begin(), pIdsNeighbors.end(), pIdsSet.begin(), + pIdsSet.end(), + std::inserter(pIdsNeighborsInterior, pIdsNeighborsInterior.begin())); + // Compute distances and extract curvature values. + std::vector curvs; + std::vector dists; + for (auto const pIdN : pIdsNeighborsInterior) + { + curvs.push_back(curvatures[pIdN]); + dists.push_back(ComputeDistance(pIdN, pId)); + } + std::vector nonZeroDistIds; + for (size_t i = 0; i < dists.size(); ++i) + { + if (dists[i] > 0) + { + nonZeroDistIds.push_back(i); + } + } + std::vector curvsNonZero; + std::vector distsNonZero; + for (auto const i : nonZeroDistIds) + { + curvsNonZero.push_back(curvs[i]); + distsNonZero.push_back(dists[i]); + } + // Iterate over the edge points and compute the curvature as the weighted + // average of the neighbours. + auto countInvalid = 0; + auto newCurv = 0.0; + if (curvsNonZero.size() > 0) + { + std::vector weights; + double sum = 0.0; + for (auto const d : distsNonZero) + { + sum += 1.0 / d; + weights.push_back(1.0 / d); + } + for (size_t i = 0; i < weights.size(); ++i) + { + weights[i] = weights[i] / sum; + } + newCurv = std::inner_product(curvsNonZero.begin(), curvsNonZero.end(), + weights.begin(), 0.0); + } + else + { + // Corner case. + // countInvalid += 1; + // Assuming the curvature of the point is planar. + newCurv = 0.0; + } + // Set the new curvature value. + curvatures[pId] = newCurv; + } - return cleaner->GetOutput(); -} + // Set small values to zero. + if (epsilon != 0.0) + { + auto eps = std::abs(epsilon); + for (size_t i = 0; i < curvatures.size(); ++i) + { + if (std::abs(curvatures[i]) < eps) + { + curvatures[i] = 0.0; + } + } + } -vtkSmartPointer GetSource(std::string const& source) -{ - std::string surface = source; - std::transform(surface.begin(), surface.end(), surface.begin(), - [](unsigned char c) { return std::tolower(c); }); - if (surface == "hills") + if (static_cast(source->GetNumberOfPoints()) != curvatures.size()) { - return GetHills(); + std::string s = curvatureName; + s += ":\nCannot add the adjusted curvatures to the source.\n"; + s += " The number of points in source does not equal the\n"; + s += " number of point ids in the adjusted curvature array."; + std::cerr << s << std::endl; + return; } - else if (surface == "parametric hills") + vtkNew adjustedCurvatures; + adjustedCurvatures->SetName(curvatureName.c_str()); + for (auto curvature : curvatures) { - return GetParametricHills(); + adjustedCurvatures->InsertNextTuple1(curvature); } - else if (surface == "parametric torus") + source->GetPointData()->AddArray(adjustedCurvatures); + source->GetPointData()->SetActiveScalars(curvatureName.c_str()); +} + +void ConstrainCurvatures(vtkPolyData* source, std::string const& curvatureName, + double const& lowerBound, double const& upperBound) +{ + std::array bounds{0.0, 0.0}; + if (lowerBound < upperBound) { - return GetParametricTorus(); + bounds[0] = lowerBound; + bounds[1] = upperBound; } - else if (surface == "plane") + else { - return GenerateElevations(GetPlane()); + bounds[0] = upperBound; + bounds[1] = lowerBound; } - else if (surface == "sphere") + + source->GetPointData()->SetActiveScalars(curvatureName.c_str()); + // Curvature as a vector. + auto array = source->GetPointData()->GetAbstractArray(curvatureName.c_str()); + std::vector curvatures; + for (vtkIdType i = 0; i < source->GetNumberOfPoints(); ++i) { - return GenerateElevations(GetSphere()); + curvatures.push_back(array->GetVariantValue(i).ToDouble()); } - else if (surface == "torus") + // Set upper and lower bounds. + for (size_t i = 0; i < curvatures.size(); ++i) { - return GenerateElevations(GetTorus()); + if (curvatures[i] < bounds[0]) + { + curvatures[i] = bounds[0]; + } + else + { + if (curvatures[i] > bounds[1]) + { + curvatures[i] = bounds[1]; + } + } } - std::cout << "The surface is not available." << std::endl; - std::cout << "Using parametric hills instead." << std::endl; - return GetParametricHills(); + vtkNew adjustedCurvatures; + for (auto curvature : curvatures) + { + adjustedCurvatures->InsertNextTuple1(curvature); + } + adjustedCurvatures->SetName(curvatureName.c_str()); + source->GetPointData()->RemoveArray(curvatureName.c_str()); + source->GetPointData()->AddArray(adjustedCurvatures); + source->GetPointData()->SetActiveScalars(curvatureName.c_str()); } -vtkNew ReverseLUT(vtkLookupTable* lut) +void AdjustRanges(std::map>& bands, + std::map& freq) { - // First do a deep copy just to get the whole structure - // and then reverse the colors and annotations. - vtkNew lutr; - lutr->DeepCopy(lut); - vtkIdType t = lut->GetNumberOfTableValues() - 1; - for (vtkIdType i = t; i >= 0; --i) + // Get the indices of the first and last non-zero elements. + auto first = 0; + for (auto i = 0; i < static_cast(freq.size()); ++i) { - std::array rgb{0.0, 0.0, 0.0}; - std::array rgba{0.0, 0.0, 0.0, 1.0}; - lut->GetColor(i, rgb.data()); - std::copy(std::begin(rgb), std::end(rgb), std::begin(rgba)); - rgba[3] = lut->GetOpacity(i); - lutr->SetTableValue(t - i, rgba.data()); + if (freq[i] != 0) + { + first = i; + break; + } + } + std::vector keys; + for (std::map::iterator it = freq.begin(); it != freq.end(); ++it) + { + keys.push_back(it->first); + } + std::reverse(keys.begin(), keys.end()); + auto last = keys[0]; + for (size_t i = 0; i < keys.size(); ++i) + { + if (freq[keys[i]] != 0) + { + last = keys[i]; + break; + } + } + // Now adjust the ranges. + std::map::iterator freqItr; + freqItr = freq.find(first); + freq.erase(freq.begin(), freqItr); + freqItr = ++freq.find(last); + freq.erase(freqItr, freq.end()); + std::map>::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 adjFreq; + int idx = 0; + for (auto p : freq) + { + adjFreq[idx] = p.second; + ++idx; } - t = lut->GetNumberOfAnnotatedValues() - 1; - for (vtkIdType i = t; i >= 0; --i) + std::map> adjBands; + idx = 0; + for (auto const& p : bands) { - lutr->SetAnnotation(t - i, lut->GetAnnotation(i)); + adjBands[idx] = p.second; + ++idx; } - - return lutr; + bands = adjBands; + freq = adjFreq; } -vtkNew GetGlyphs(vtkPolyData* src, double const& scaleFactor, - bool const& reverseNormals) +BandedGlyphs GetCurvatureGlyphs(std::string const& surfaceName, + vtkPolyData* source, + std::string const& curvature, + unsigned int precision, bool frequencyTable, + bool nearestInteger) { - // 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 reverse; - vtkNew maskPts; - maskPts->SetOnRatio(5); - maskPts->RandomModeOn(); - if (reverseNormals) + // The length of the normal arrow glyphs. + auto scaleFactor = 1.0; + if (surfaceName == "hills") { - reverse->SetInputData(src); - reverse->ReverseCellsOn(); - reverse->ReverseNormalsOn(); - maskPts->SetInputConnection(reverse->GetOutputPort()); + scaleFactor = 0.5; } - else + if (surfaceName == "sphere") { - maskPts->SetInputData(src); + scaleFactor = 0.2; } - // Source for the glyph filter - vtkNew arrow; - arrow->SetTipResolution(16); - arrow->SetTipLength(0.3); - arrow->SetTipRadius(0.1); - - vtkNew 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> -GetBands(std::array const& scalarRange, int const& numberOfBands, - int const& precision, bool const& nearestInteger) -{ - auto prec = abs(precision); - prec = (prec > 14) ? 14 : prec; + source->GetPointData()->SetActiveScalars(curvature.c_str()); + auto sr = source->GetPointData()->GetScalars(curvature.c_str())->GetRange(); + std::array scalarRange{sr[0], sr[1]}; - auto RoundOff = [&prec](const double& x) { - auto pow_10 = std::pow(10.0, prec); - return std::round(x * pow_10) / pow_10; - }; + int maxBands = 0; - std::map> bands; - if ((scalarRange[1] < scalarRange[0]) || (numberOfBands <= 0)) + vtkNew colorSeries; + if (surfaceName == "hills") { - return bands; + colorSeries->SetColorScheme(colorSeries->BREWER_DIVERGING_SPECTRAL_7); + maxBands = 7; } - double x[2]; - for (int i = 0; i < 2; ++i) + else if (surfaceName == "parametric hills") { - x[i] = scalarRange[i]; + colorSeries->SetColorScheme(colorSeries->BREWER_DIVERGING_SPECTRAL_9); + maxBands = 9; } - if (nearestInteger) + else if (surfaceName == "plane" || surfaceName == "sphere") { - x[0] = std::floor(x[0]); - x[1] = std::ceil(x[1]); + colorSeries->SetColorScheme( + colorSeries->BREWER_SEQUENTIAL_YELLOW_ORANGE_BROWN_3); + maxBands = 3; } - double dx = (x[1] - x[0]) / static_cast(numberOfBands); - std::vector 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) + else { - if (i == 0) - { - for (std::vector::iterator p = b.begin(); p != b.end(); ++p) - { - *p = RoundOff(*p); - } - b[0] = x[0]; - } - bands[i] = b; - for (std::vector::iterator p = b.begin(); p != b.end(); ++p) - { - *p = RoundOff(*p + dx); - } + colorSeries->SetColorScheme(colorSeries->BREWER_DIVERGING_SPECTRAL_6); + maxBands = 6; } - return bands; -} -std::map> -GetCustomBands(std::array const& scalarRange, - int const& numberOfBands, - std::vector> const& myBands) -{ - std::map> bands; - if ((scalarRange[1] < scalarRange[0]) || (numberOfBands <= 0)) + vtkNew lut; + colorSeries->BuildLookupTable(lut, colorSeries->CATEGORICAL); + lut->SetNanColor(0, 0, 0, 1); + lut->SetTableRange(scalarRange.data()); + + if (surfaceName == "plane" || surfaceName == "sphere") { - return bands; + lut->SetNumberOfTableValues(1); } + auto numberOfBands = lut->GetNumberOfTableValues(); - std::vector> x; - std::copy(myBands.begin(), myBands.end(), std::back_inserter(x)); + vtkNew lut1; + colorSeries->BuildLookupTable(lut1, colorSeries->ORDINAL); + lut1->SetNanColor(0, 0, 0, 1); + lut1->SetTableRange(scalarRange.data()); + lut1->SetNumberOfTableValues(numberOfBands); - // Determine the index of the range minimum and range maximum. - int idxMin = 0; - for (auto idx = 0; idx < static_cast(myBands.size()); ++idx) + auto bands = GetBands(scalarRange, numberOfBands, precision, nearestInteger); + + if (curvature == "Gaussian_Curvature") { - if (scalarRange[0] < myBands[idx][1] && scalarRange[0] >= myBands[idx][0]) + if (surfaceName == "random hills") { - idxMin = idx; - break; + // These are my custom bands. + // Generated by first running: + // bands = GetBands(curvScalarRanges, curvNumBands, precision, false); + // then: + // std::vector freq = Frequencies(bands, src); + // PrintBandsFrequencies(bands, freq); + // Finally using the output to create this table: + // std::vector> 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},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) + // and some spherical regions. + std::vector> myBands = { + {-0.630, -0.190}, {-0.190, -0.043}, {-0.043, 0.0452}, + {0.0452, 0.0746}, {0.0746, 0.104}, {0.104, 0.251}, + {0.251, 1.131}}; + + // Comment this out if you want to see how allocating + // equally spaced bands works. + bands = GetCustomBands(scalarRange, numberOfBands, myBands); } - } - int idxMax = static_cast(myBands.size()) - 1; - for (int idx = static_cast(myBands.size()) - 1; idx >= 0; --idx) - { - if (scalarRange[1] < myBands[idx][1] && scalarRange[1] >= myBands[idx][0]) + else if (surfaceName == "hills") { - idxMax = static_cast(idx); - break; + std::vector> myBands = { + {-2.104, -0.15}, {-0.15, -0.1}, {-0.1, -0.05}, + {-0.05, -0.02}, {-0.02, -0.005}, {-0.005, -0.0005}, + {-0.0005, 0.0005}, {0.0005, 0.09}, {0.09, 4.972}, + }; + + // Comment this out if you want to see how allocating + // equally spaced bands works. + bands = GetCustomBands(scalarRange, numberOfBands, myBands); } } - // Set the minimum to match the range minimum. - x[idxMin][0] = scalarRange[0]; - x[idxMax][1] = scalarRange[1]; - for (int i = idxMin; i < idxMax + 1; ++i) + // Adjust the number of table values and scalar range. + scalarRange[0] = bands.begin()->second[0]; + scalarRange[1] = std::prev(bands.end())->second[2]; + lut->SetTableRange(scalarRange.data()); + lut->SetNumberOfTableValues(bands.size()); + lut1->SetTableRange(scalarRange.data()); + lut1->SetNumberOfTableValues(bands.size()); + + if (frequencyTable) { - std::vector 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; + auto freq = GetFrequencies(bands, source); + AdjustRanges(bands, freq); + fmt::println("{:s} {:s}", Title(surfaceName), curvature); + PrintBandsFrequencies(bands, freq); } - return bands; -} -std::map> -GetIntegralBands(std::array const& scalarRange) -{ - if (scalarRange[1] < scalarRange[0]) + // We will use the midpoint of the band as the label. + std::vector labels; + for (std::map>::const_iterator p = bands.begin(); + p != bands.end(); ++p) { - std::map> bands; - return bands; + labels.push_back(fmt::format("{:4.2f}", p->second[1])); } - std::array x{0, 0}; - for (int i = 0; i < 2; ++i) + + // Annotate + vtkNew values; + for (size_t i = 0; i < labels.size(); ++i) { - x[i] = scalarRange[i]; + values->InsertNextValue(vtkVariant(labels[i])); } - x[0] = std::floor(x[0]); - x[1] = std::ceil(x[1]); - int numberOfBands = static_cast(std::abs(x[1]) - std::abs(x[0])); - return GetBands(x, numberOfBands, false); + for (vtkIdType i = 0; i < values->GetNumberOfTuples(); ++i) + { + lut->SetAnnotation(i, values->GetValue(i).ToString()); + } + + // Create the contour bands. + vtkNew bcf; + bcf->SetInputData(source); + // Use either the minimum or maximum value for each band. + int i = 0; + for (std::map>::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 glyphs = GetGlyphs(source, scaleFactor, false); + + BandedGlyphs bg; + bg.glyphType = curvature; + bg.glyphs = glyphs; + bg.lutMaxBands = maxBands; + bg.bcf = bcf; + bg.lut = lut; + bg.lutr = ReverseLUT(lut); + bg.lut1MaxBands = maxBands; + bg.lut1 = lut1; + bg.lut1r = ReverseLUT(lut1); + bg.scalarRange = scalarRange; + bg.labels = labels; + + return bg; } -std::map GetFrequencies(std::map>& bands, - vtkPolyData* src) +BandedGlyphs GetElevationGlyphs(std::string const& surfaceName, + vtkPolyData* source, unsigned int precision, + bool frequencyTable, bool nearestInteger) { - std::map freq; - for (auto i = 0; i < static_cast(bands.size()); ++i) + // The length of the normal arrow glyphs. + auto scaleFactor = 1.0; + if (surfaceName == "hills") { - freq[i] = 0; + scaleFactor = 0.5; } - vtkIdType tuples = src->GetPointData()->GetScalars()->GetNumberOfTuples(); - for (int i = 0; i < tuples; ++i) + if (surfaceName == "sphere") { - const double* x = src->GetPointData()->GetScalars()->GetTuple(i); - for (auto j = 0; j < static_cast(bands.size()); ++j) - { - if (*x <= bands[j][2]) - { - freq[j] = freq[j] + 1; - break; - } - } + scaleFactor = 0.2; } - return freq; -} -void AdjustFrequencyRanges(std::map>& bands, - std::map& freq) -{ - // Get the indices of the first and last non-zero elements. - auto first = 0; - for (auto i = 0; i < static_cast(freq.size()); ++i) + source->GetPointData()->SetActiveScalars("Elevation"); + auto sr = source->GetPointData()->GetScalars("Elevation")->GetRange(); + std::array scalarRange{sr[0], sr[1]}; + + int maxBands = 0; + + vtkNew colorSeries; + if (surfaceName == "hills" || surfaceName == "parametric hills") { - if (freq[i] != 0) - { - first = i; - break; - } + colorSeries->SetColorScheme( + colorSeries->BREWER_DIVERGING_BROWN_BLUE_GREEN_8); + maxBands = 8; } - std::vector keys; - for (std::map::iterator it = freq.begin(); it != freq.end(); ++it) + else { - keys.push_back(it->first); + colorSeries->SetColorScheme( + colorSeries->BREWER_DIVERGING_BROWN_BLUE_GREEN_6); + maxBands = 6; } - std::reverse(keys.begin(), keys.end()); - auto last = keys[0]; - for (size_t i = 0; i < keys.size(); ++i) + + vtkNew lut; + colorSeries->BuildLookupTable(lut, colorSeries->CATEGORICAL); + lut->SetNanColor(0, 0, 0, 1); + lut->SetTableRange(scalarRange.data()); + + if (surfaceName == "plane") { - if (freq[keys[i]] != 0) - { - last = keys[i]; - break; - } + lut->SetNumberOfTableValues(1); } - // Now adjust the ranges. - std::map::iterator freqItr; - freqItr = freq.find(first); - freq.erase(freq.begin(), freqItr); - freqItr = ++freq.find(last); - freq.erase(freqItr, freq.end()); - std::map>::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 adjFreq; - int idx = 0; - for (auto p : freq) + auto numberOfBands = lut->GetNumberOfTableValues(); + + vtkNew lut1; + colorSeries->BuildLookupTable(lut1, colorSeries->ORDINAL); + lut1->SetNanColor(0, 0, 0, 1); + lut1->SetTableRange(scalarRange.data()); + lut1->SetNumberOfTableValues(numberOfBands); + + auto bands = GetBands(scalarRange, numberOfBands, precision, nearestInteger); + + if (surfaceName == "parametric hills") { - adjFreq[idx] = p.second; - ++idx; + // These are my custom bands. + // Generated by first running: + // bands = GetBands(scalarRange, numberOfBands, precision, false); + // then: + // std::vector freq = Frequencies(bands, source); + // PrintBandsFrequencies(bands, freq); + // Finally using the output to create this table: + std::vector> myBands = { + {0, 1.0}, {1.0, 2.0}, {2.0, 3.0}, {3.0, 4.0}, + {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 = GetCustomBands(scalarRange, numberOfBands, myBands); } - std::map> adjBands; - idx = 0; - for (auto const& p : bands) + + // Adjust the number of table values and scalar range. + scalarRange[0] = bands.begin()->second[0]; + scalarRange[1] = std::prev(bands.end())->second[2]; + lut->SetTableRange(scalarRange.data()); + lut->SetNumberOfTableValues(bands.size()); + lut1->SetTableRange(scalarRange.data()); + lut1->SetNumberOfTableValues(bands.size()); + + if (frequencyTable) { - adjBands[idx] = p.second; - ++idx; + auto freq = GetFrequencies(bands, source); + AdjustRanges(bands, freq); + fmt::println("{:s} Elevation", Title(surfaceName)); + PrintBandsFrequencies(bands, freq); } - bands = adjBands; - freq = adjFreq; -} -void PrintBandsFrequencies(std::map> const& bands, - std::map& freq, int const& precision) -{ - auto prec = abs(precision); - prec = (prec > 14) ? 14 : prec; + // We will use the midpoint of the band as the label. + std::vector labels; + for (std::map>::const_iterator p = bands.begin(); + p != bands.end(); ++p) + { + labels.push_back(fmt::format("{:4.2f}", p->second[1])); + } - if (bands.size() != freq.size()) + // Annotate + vtkNew values; + for (size_t i = 0; i < labels.size(); ++i) { - std::cout << "Bands and frequencies must be the same size." << std::endl; - return; + values->InsertNextValue(vtkVariant(labels[i])); } - std::ostringstream os; - os << "Bands & Frequencies:\n"; - size_t idx = 0; - auto total = 0; - auto width = prec + 6; - std::string res{""}; + for (vtkIdType i = 0; i < values->GetNumberOfTuples(); ++i) + { + lut->SetAnnotation(i, values->GetValue(i).ToString()); + } + + // Create the contour bands. + vtkNew bcf; + bcf->SetInputData(source); + // Use either the minimum or maximum value for each band. + int i = 0; for (std::map>::const_iterator p = bands.begin(); p != bands.end(); ++p) { - total += freq[p->first]; - for (std::vector::const_iterator q = p->second.begin(); - q != p->second.end(); ++q) - { - if (q == p->second.begin()) - { - res += fmt::format("{:4d} [", idx); - } - if (q == std::prev(p->second.end())) - { - res += fmt::format("{:{}.{}f}]: {:{}d}\n", *q, width, precision, - freq[p->first], width); - } - else - { - res += fmt::format("{:{}.{}f}, ", *q, width, precision); - } - } - ++idx; + bcf->SetValue(i, p->second[2]); + ++i; } - auto width1 = 3 * width + 13; - fmt::println("{:s}\n{:<{}s}{:>{}d}", res, "Total", width1, total, width); - std::cout << std::endl; + // We will use an indexed lookup table. + bcf->SetScalarModeToIndex(); + bcf->GenerateContourEdgesOn(); + + // Generate the glyphs on the original surface. + auto glyphs = GetGlyphs(source, scaleFactor, false); + + BandedGlyphs bg; + bg.glyphType = "Elevation"; + bg.glyphs = glyphs; + bg.lutMaxBands = maxBands; + bg.bcf = bcf; + bg.lut = lut; + bg.lutr = ReverseLUT(lut); + bg.lut1MaxBands = maxBands; + bg.lut1 = lut1; + bg.lut1r = ReverseLUT(lut1); + bg.scalarRange = scalarRange; + bg.labels = labels; + + return bg; } vtkNew @@ -2009,6 +1969,51 @@ std::string Title(std::string s, const std::locale& loc) return s; } +void PrintBandsFrequencies(std::map> const& bands, + std::map& 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"; + size_t idx = 0; + auto total = 0; + auto width = prec + 6; + std::string res{""}; + for (std::map>::const_iterator p = bands.begin(); + p != bands.end(); ++p) + { + total += freq[p->first]; + for (std::vector::const_iterator q = p->second.begin(); + q != p->second.end(); ++q) + { + if (q == p->second.begin()) + { + res += fmt::format("{:4d} [", idx); + } + if (q == std::prev(p->second.end())) + { + res += fmt::format("{:{}.{}f}]: {:{}d}\n", *q, width, precision, + freq[p->first], width); + } + else + { + res += fmt::format("{:{}.{}f}, ", *q, width, precision); + } + } + ++idx; + } + auto width1 = 3 * width + 13; + fmt::println("{:s}\n{:<{}s}{:>{}d}", res, "Total", width1, total, width); + std::cout << std::endl; +} + void AdjustCameraParameters(std::string const& surfaceName, vtkRenderer* ren) { if (surfaceName == "hills") diff --git a/src/Cxx/Visualization/ElevationBandsWithGlyphs.cxx b/src/Cxx/Visualization/ElevationBandsWithGlyphs.cxx index e04e3d4615e..af41d8c560c 100644 --- a/src/Cxx/Visualization/ElevationBandsWithGlyphs.cxx +++ b/src/Cxx/Visualization/ElevationBandsWithGlyphs.cxx @@ -56,6 +56,14 @@ namespace fs = std::filesystem; namespace { +/** + * Generate elevations over the surface. + * + * @param src - The vtkPolyData source. + * @return elev - The elevations. + */ +vtkSmartPointer GenerateElevations(vtkPolyData* src); + vtkSmartPointer GetHills(); vtkSmartPointer GetParametricHills(); vtkSmartPointer GetParametricTorus(); @@ -66,47 +74,6 @@ vtkSmartPointer GetSource(std::string const& source); vtkSmartPointer ReverseLUT(vtkLookupTable* lut); -struct BandedGlyphs -{ - std::string glyphType; - vtkSmartPointer glyphs; - vtkSmartPointer bcf; - vtkSmartPointer lut; - vtkSmartPointer lutr; - vtkSmartPointer lut1; - vtkSmartPointer lut1r; - std::array scalarRange; - std::vector labels; -}; - -/** - * Get elevation glyphs and the corresponding banded polydata filter for the - * surface. - * - * @param surfaceName - The surface name. - * @param source - A vtkPolyData object corresponding to the source. - * @param needsAdjusting - Surfaces whose curvatures need to be adjusted - * along the edges of the surface or constrained. - * @param frequencyTable - If true, display a frequency table corresponding to - * the bands. - * @param nearestInteger - If true, use the nearest integer when generating the - * bands. - * @return A struct holding glyphs, bcf, lut, lutr, lut1, lut1r scalar_range, - * labels. - */ -BandedGlyphs GetElevationGlyphs(std::string const& surfaceName, - vtkPolyData* source, unsigned int precision, - bool frequencyTable = false, - bool nearestInteger = false); - -/** - * Generate elevations over the surface. - * - * @param src - The vtkPolyData source. - * @return elev - The elevations. - */ -vtkSmartPointer GenerateElevations(vtkPolyData* src); - /** * Glyph the normals on the surface. * @@ -183,18 +150,43 @@ std::map GetFrequencies(std::map>& bands, * @param bands: The bands. * @param freq: The frequencies. */ -void AdjustFrequencyRanges(std::map>& bands, - std::map& freq); +void AdjustRanges(std::map>& bands, + std::map& freq); + +struct BandedGlyphs +{ + std::string glyphType; + vtkSmartPointer glyphs; + vtkSmartPointer bcf; + int lutMaxBands; + vtkSmartPointer lut; + vtkSmartPointer lutr; + int lut1MaxBands; + vtkSmartPointer lut1; + vtkSmartPointer lut1r; + std::array scalarRange; + std::vector labels; +}; /** - * Print each band and the number of scalars in each band. + * Get elevation glyphs and the corresponding banded polydata filter for the + * surface. * - * @param bands: The bands. - * @param freq: The frequencies. - * @param precision: The precision for the ranges in each band. + * @param surfaceName - The surface name. + * @param source - A vtkPolyData object corresponding to the source. + * @param needsAdjusting - Surfaces whose curvatures need to be adjusted + * along the edges of the surface or constrained. + * @param frequencyTable - If true, display a frequency table corresponding to + * the bands. + * @param nearestInteger - If true, use the nearest integer when generating the + * bands. + * @return A struct holding glyphs, bcf, lut, lutr, lut1, lut1r scalar_range, + * labels. */ -void PrintBandsFrequencies(std::map> const& bands, - std::map& freq, int const& precision = 2); +BandedGlyphs GetElevationGlyphs(std::string const& surfaceName, + vtkPolyData* source, unsigned int precision, + bool frequencyTable = false, + bool nearestInteger = false); typedef std::map> TTextPosition; typedef std::map TTextPositions; @@ -288,6 +280,16 @@ GetTextPositions(std::vector const& names, */ std::string Title(std::string s, const std::locale& loc = std::locale()); +/** + * Print each band and the number of scalars in each band. + * + * @param bands: The bands. + * @param freq: The frequencies. + * @param precision: The precision for the ranges in each band. + */ +void PrintBandsFrequencies(std::map> const& bands, + std::map& freq, int const& precision = 2); + /** * Adjust the camera parameters. * @@ -477,25 +479,14 @@ int main(int argc, char* argv[]) textWidget->EnabledOn(); auto elevSBP = ScalarBarProperties(); - elevSBP.titleText = "Elevation\n"; + elevSBP.titleText = elevBG.glyphType + "\n"; + std::replace(elevSBP.titleText.begin(), elevSBP.titleText.end(), '_', '\n'); elevSBP.numberOfLabels = static_cast(elevBG.labels.size()); // lut puts the lowest value at the top of the vertical scalar bar. // lutr puts the highest value at the top of the vertical scalar bar. elevSBP.lut = elevBG.lutr; elevSBP.orientation = true; - auto maxBands = 8; - if (surfaceName == "hills" || surfaceName == "parametric hills") - { - elevSBP.positionV = std::move(PositionSBWV(8, maxBands)); - } - else if (surfaceName == "plane") - { - elevSBP.positionV = std::move(PositionSBWV(1, maxBands)); - } - else - { - elevSBP.positionV = std::move(PositionSBWV(5, maxBands)); - } + elevSBP.positionV = PositionSBWV(elevSBP.numberOfLabels, elevBG.lut1MaxBands); auto sbw = MakeScalarBarWidget(elevSBP, titleTextProperty, labelTextProperty, ren, iren); @@ -535,138 +526,6 @@ int main(int argc, char* argv[]) } namespace { -BandedGlyphs GetElevationGlyphs(std::string const& surfaceName, - vtkPolyData* source, unsigned int precision, - bool frequencyTable, bool nearestInteger) -{ - // The length of the normal arrow glyphs. - auto scaleFactor = 1.0; - if (surfaceName == "hills") - { - scaleFactor = 0.5; - } - if (surfaceName == "sphere") - { - scaleFactor = 0.2; - } - - source->GetPointData()->SetActiveScalars("Elevation"); - auto sr = source->GetPointData()->GetScalars("Elevation")->GetRange(); - std::array scalarRange{sr[0], sr[1]}; - - vtkNew colorSeries; - if (surfaceName == "hills" || surfaceName == "parametric hills") - { - colorSeries->SetColorScheme( - colorSeries->BREWER_DIVERGING_BROWN_BLUE_GREEN_8); - } - else - { - colorSeries->SetColorScheme( - colorSeries->BREWER_DIVERGING_BROWN_BLUE_GREEN_5); - } - - vtkNew lut; - colorSeries->BuildLookupTable(lut, colorSeries->CATEGORICAL); - lut->SetNanColor(0, 0, 0, 1); - lut->SetTableRange(scalarRange.data()); - - vtkNew lut1; - colorSeries->BuildLookupTable(lut1, colorSeries->ORDINAL); - lut1->SetNanColor(0, 0, 0, 1); - lut1->SetTableRange(scalarRange.data()); - - auto numberOfBands = lut->GetNumberOfTableValues(); - lut1->SetNumberOfTableValues(numberOfBands); - - auto bands = GetBands(scalarRange, numberOfBands, precision, nearestInteger); - - if (surfaceName == "parametric hills") - { - // These are my custom bands. - // Generated by first running: - // bands = GetBands(scalarRange, numberOfBands, precision, false); - // then: - // std::vector freq = Frequencies(bands, source); - // PrintBandsFrequencies(bands, freq); - // Finally using the output to create this table: - std::vector> myBands = { - {0, 1.0}, {1.0, 2.0}, {2.0, 3.0}, {3.0, 4.0}, - {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 = GetCustomBands(scalarRange, numberOfBands, myBands); - // bands = GetBands(scalarRange, numberOfBands, precision, false); - // Adjust the number of table values - } - - // Adjust the number of table values and scalar range. - scalarRange[0] = bands.begin()->second[0]; - scalarRange[1] = std::prev(bands.end())->second[2]; - lut->SetTableRange(scalarRange.data()); - lut->SetNumberOfTableValues(bands.size()); - lut1->SetTableRange(scalarRange.data()); - lut1->SetNumberOfTableValues(bands.size()); - - if (frequencyTable) - { - auto freq = GetFrequencies(bands, source); - AdjustFrequencyRanges(bands, freq); - fmt::println("{:s} Elevation", Title(surfaceName)); - PrintBandsFrequencies(bands, freq); - } - - // We will use the midpoint of the band as the label. - std::vector labels; - for (std::map>::const_iterator p = bands.begin(); - p != bands.end(); ++p) - { - labels.push_back(fmt::format("{:4.2f}", p->second[1])); - } - - // Annotate - vtkNew values; - for (size_t i = 0; i < labels.size(); ++i) - { - values->InsertNextValue(vtkVariant(labels[i])); - } - for (vtkIdType i = 0; i < values->GetNumberOfTuples(); ++i) - { - lut->SetAnnotation(i, values->GetValue(i).ToString()); - } - - // Create the contour bands. - vtkNew bcf; - bcf->SetInputData(source); - // Use either the minimum or maximum value for each band. - int i = 0; - for (std::map>::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 glyphs = GetGlyphs(source, scaleFactor, false); - - BandedGlyphs bg; - bg.glyphType = "Elevation"; - bg.glyphs = glyphs; - bg.bcf = bcf; - bg.lut = lut; - bg.lutr = ReverseLUT(lut); - bg.lut1 = lut1; - bg.lut1r = ReverseLUT(lut1); - bg.scalarRange = scalarRange; - bg.labels = labels; - - return bg; -} - vtkSmartPointer GenerateElevations(vtkPolyData* src) { double bounds[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0}; @@ -1171,8 +1030,8 @@ std::map GetFrequencies(std::map>& bands, return freq; } -void AdjustFrequencyRanges(std::map>& bands, - std::map& freq) +void AdjustRanges(std::map>& bands, + std::map& freq) { // Get the indices of the first and last non-zero elements. auto first = 0; @@ -1229,49 +1088,144 @@ void AdjustFrequencyRanges(std::map>& bands, freq = adjFreq; } -void PrintBandsFrequencies(std::map> const& bands, - std::map& freq, int const& precision) +BandedGlyphs GetElevationGlyphs(std::string const& surfaceName, + vtkPolyData* source, unsigned int precision, + bool frequencyTable, bool nearestInteger) { - auto prec = abs(precision); - prec = (prec > 14) ? 14 : prec; + // The length of the normal arrow glyphs. + auto scaleFactor = 1.0; + if (surfaceName == "hills") + { + scaleFactor = 0.5; + } + if (surfaceName == "sphere") + { + scaleFactor = 0.2; + } - if (bands.size() != freq.size()) + source->GetPointData()->SetActiveScalars("Elevation"); + auto sr = source->GetPointData()->GetScalars("Elevation")->GetRange(); + std::array scalarRange{sr[0], sr[1]}; + + int maxBands = 0; + + vtkNew colorSeries; + if (surfaceName == "hills" || surfaceName == "parametric hills") { - std::cout << "Bands and frequencies must be the same size." << std::endl; - return; + colorSeries->SetColorScheme( + colorSeries->BREWER_DIVERGING_BROWN_BLUE_GREEN_8); + maxBands = 8; } - std::ostringstream os; - os << "Bands & Frequencies:\n"; - size_t idx = 0; - auto total = 0; - auto width = prec + 6; - std::string res{""}; + else + { + colorSeries->SetColorScheme( + colorSeries->BREWER_DIVERGING_BROWN_BLUE_GREEN_6); + maxBands = 6; + } + + vtkNew lut; + colorSeries->BuildLookupTable(lut, colorSeries->CATEGORICAL); + lut->SetNanColor(0, 0, 0, 1); + lut->SetTableRange(scalarRange.data()); + + if (surfaceName == "plane") + { + lut->SetNumberOfTableValues(1); + } + auto numberOfBands = lut->GetNumberOfTableValues(); + + vtkNew lut1; + colorSeries->BuildLookupTable(lut1, colorSeries->ORDINAL); + lut1->SetNanColor(0, 0, 0, 1); + lut1->SetTableRange(scalarRange.data()); + lut1->SetNumberOfTableValues(numberOfBands); + + auto bands = GetBands(scalarRange, numberOfBands, precision, nearestInteger); + + if (surfaceName == "parametric hills") + { + // These are my custom bands. + // Generated by first running: + // bands = GetBands(scalarRange, numberOfBands, precision, false); + // then: + // std::vector freq = Frequencies(bands, source); + // PrintBandsFrequencies(bands, freq); + // Finally using the output to create this table: + std::vector> myBands = { + {0, 1.0}, {1.0, 2.0}, {2.0, 3.0}, {3.0, 4.0}, + {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 = GetCustomBands(scalarRange, numberOfBands, myBands); + } + + // Adjust the number of table values and scalar range. + scalarRange[0] = bands.begin()->second[0]; + scalarRange[1] = std::prev(bands.end())->second[2]; + lut->SetTableRange(scalarRange.data()); + lut->SetNumberOfTableValues(bands.size()); + lut1->SetTableRange(scalarRange.data()); + lut1->SetNumberOfTableValues(bands.size()); + + if (frequencyTable) + { + auto freq = GetFrequencies(bands, source); + AdjustRanges(bands, freq); + fmt::println("{:s} Elevation", Title(surfaceName)); + PrintBandsFrequencies(bands, freq); + } + + // We will use the midpoint of the band as the label. + std::vector labels; for (std::map>::const_iterator p = bands.begin(); p != bands.end(); ++p) { - total += freq[p->first]; - for (std::vector::const_iterator q = p->second.begin(); - q != p->second.end(); ++q) - { - if (q == p->second.begin()) - { - res += fmt::format("{:4d} [", idx); - } - if (q == std::prev(p->second.end())) - { - res += fmt::format("{:{}.{}f}]: {:{}d}\n", *q, width, precision, - freq[p->first], width); - } - else - { - res += fmt::format("{:{}.{}f}, ", *q, width, precision); - } - } - ++idx; + labels.push_back(fmt::format("{:4.2f}", p->second[1])); } - auto width1 = 3 * width + 13; - fmt::println("{:s}\n{:<{}s}{:>{}d}", res, "Total", width1, total, width); - std::cout << std::endl; + + // Annotate + vtkNew values; + for (size_t i = 0; i < labels.size(); ++i) + { + values->InsertNextValue(vtkVariant(labels[i])); + } + for (vtkIdType i = 0; i < values->GetNumberOfTuples(); ++i) + { + lut->SetAnnotation(i, values->GetValue(i).ToString()); + } + + // Create the contour bands. + vtkNew bcf; + bcf->SetInputData(source); + // Use either the minimum or maximum value for each band. + int i = 0; + for (std::map>::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 glyphs = GetGlyphs(source, scaleFactor, false); + + BandedGlyphs bg; + bg.glyphType = "Elevation"; + bg.glyphs = glyphs; + bg.lutMaxBands = maxBands; + bg.bcf = bcf; + bg.lut = lut; + bg.lutr = ReverseLUT(lut); + bg.lut1MaxBands = maxBands; + bg.lut1 = lut1; + bg.lut1r = ReverseLUT(lut1); + bg.scalarRange = scalarRange; + bg.labels = labels; + + return bg; } vtkNew @@ -1473,6 +1427,51 @@ std::string Title(std::string s, const std::locale& loc) return s; } +void PrintBandsFrequencies(std::map> const& bands, + std::map& 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"; + size_t idx = 0; + auto total = 0; + auto width = prec + 6; + std::string res{""}; + for (std::map>::const_iterator p = bands.begin(); + p != bands.end(); ++p) + { + total += freq[p->first]; + for (std::vector::const_iterator q = p->second.begin(); + q != p->second.end(); ++q) + { + if (q == p->second.begin()) + { + res += fmt::format("{:4d} [", idx); + } + if (q == std::prev(p->second.end())) + { + res += fmt::format("{:{}.{}f}]: {:{}d}\n", *q, width, precision, + freq[p->first], width); + } + else + { + res += fmt::format("{:{}.{}f}, ", *q, width, precision); + } + } + ++idx; + } + auto width1 = 3 * width + 13; + fmt::println("{:s}\n{:<{}s}{:>{}d}", res, "Total", width1, total, width); + std::cout << std::endl; +} + void AdjustCameraParameters(std::string const& surfaceName, vtkRenderer* ren) { if (surfaceName == "hills") diff --git a/src/PythonicAPI/PolyData/CurvaturesNormalsElevations.py b/src/PythonicAPI/PolyData/CurvaturesNormalsElevations.py index 40bc7d2b8e2..2a5f64a3d68 100644 --- a/src/PythonicAPI/PolyData/CurvaturesNormalsElevations.py +++ b/src/PythonicAPI/PolyData/CurvaturesNormalsElevations.py @@ -23,6 +23,7 @@ from vtkmodules.vtkCommonComputationalGeometry import ( ) from vtkmodules.vtkCommonCore import ( VTK_DOUBLE, + vtkDoubleArray, vtkFloatArray, vtkIdList, vtkLookupTable, @@ -30,20 +31,18 @@ from vtkmodules.vtkCommonCore import ( vtkVariant, vtkVariantArray ) -from vtkmodules.vtkCommonDataModel import ( - vtkCellArray, - vtkTriangle -) from vtkmodules.vtkCommonDataModel import vtkPolyData from vtkmodules.vtkCommonTransforms import vtkTransform from vtkmodules.vtkFiltersCore import ( vtkCleanPolyData, + vtkDelaunay2D, vtkElevationFilter, vtkFeatureEdges, vtkGlyph3D, vtkGenerateIds, vtkMaskPoints, vtkPolyDataNormals, + vtkPolyDataTangents, vtkReverseSense, vtkTriangleFilter ) @@ -71,7 +70,6 @@ from vtkmodules.vtkInteractionWidgets import ( from vtkmodules.vtkRenderingAnnotation import vtkAxesActor, vtkScalarBarActor from vtkmodules.vtkRenderingCore import ( vtkActor, - vtkColorTransferFunction, vtkPolyDataMapper, vtkRenderWindow, vtkRenderWindowInteractor, @@ -85,15 +83,15 @@ def get_program_parameters(): import argparse description = 'Demonstrates Gaussian and Mean curvatures on a surface, along with normals colored by elevation.' epilogue = ''' - For example: "Random Hills" -f + For example: "parametric hills" -f Will display the curvatures along with normals on the surface colored by elevation. ''' parser = argparse.ArgumentParser(description=description, epilog=epilogue, formatter_class=argparse.RawDescriptionHelpFormatter) - parser.add_argument('surface_name', nargs='?', default='random hills', + parser.add_argument('surface_name', nargs='?', default='parametric hills', help='The name of the surface - enclose the name in quotes if it has spaces.') parser.add_argument('-f', '--frequency_table', action='store_true', help='Display the frequency table.') - parser.add_argument('-omw', action='store_false', + parser.add_argument('-o', '--omw', action='store_true', help='Use an OrientationMarkerWidget instead of a CameraOrientationWidget.') args = parser.parse_args() @@ -101,15 +99,15 @@ def get_program_parameters(): def main(argv): - surface_name, frequency_table, use_camera_omw = get_program_parameters() + surface_name, frequency_table, use_omw = get_program_parameters() - available_surfaces = ['hills', 'parametric torus', 'plane', 'random hills', 'sphere', 'torus'] + available_surfaces = ['hills', 'parametric hills', 'parametric torus', 'plane', 'sphere', 'torus'] # Surfaces whose curvatures need to be adjusted along the edges of the surface or constrained. - needs_adjusting = ['hills', 'parametric torus', 'plane', 'random hills'] + needs_adjusting = ['hills', 'parametric hills', 'parametric torus', 'plane'] surface_name = ' '.join(surface_name.lower().replace('_', ' ').split()) - if surface_name == 'randomhills': - surface_name = 'random hills' + if surface_name in ['parametrichills', 'random hills', 'randomhills']: + surface_name = 'parametric hills' if surface_name == 'parametrictorus': surface_name = 'parametric torus' if surface_name.lower() not in available_surfaces: @@ -123,24 +121,55 @@ def main(argv): if i < len(asl) - 1: s += ',' print(f' {s}') - print('If a name has spaces in it, delineate the name with quotes e.g. "random hills"') + print('If a name has spaces in it, delineate the name with quotes e.g. "parametric hills"') return - Surface = namedtuple('Surface', 'name source') - surface = Surface(surface_name, get_source(surface_name, available_surfaces)) + source = get_source(surface_name) + if not source: + print('The surface is not available.') + return - # -------------------------------------------------------------------------------------- - # Get the filters, scalar range of curvatures and elevation along with the lookup tables. - # -------------------------------------------------------------------------------------- # Use an ordered dictionary as we want the keys in a specific order. curvatures = OrderedDict() - curvatures['Gauss_Curvature'] = generate_gaussian_curvatures(surface, needs_adjusting, - frequency_table=frequency_table) - curvatures['Mean_Curvature'] = generate_mean_curvatures(surface, needs_adjusting, frequency_table=frequency_table) + # Gaussian Curvature + cc_gauss = vtkCurvatures(input_data=source) + cc_gauss.SetCurvatureTypeToGaussian() + cc_gauss.update() + curvature = cc_gauss.output.point_data.scalars.name + if surface_name in needs_adjusting: + adjust_edge_curvatures(cc_gauss.output, curvature) + if surface_name == 'plane': + constrain_curvatures(cc_gauss.output, curvature, 0.0, 0.0) + if surface_name == 'sphere': + # The sphere radius is 1. + # Gaussian curvature is 1/r^2 + constrain_curvatures(cc_gauss.output, curvature, 1.0, 1.0) + curvatures[curvature] = get_curvature_glyphs(surface_name, cc_gauss, curvature, + precision=10, frequency_table=frequency_table, + nearest_integer=False) + + # Mean Curvature + cc_mean = vtkCurvatures(input_data=source) + cc_mean.SetCurvatureTypeToMean() + cc_mean.update() + curvature = cc_mean.output.point_data.scalars.name + if surface_name in needs_adjusting: + adjust_edge_curvatures(cc_mean.output, curvature) + if surface_name == 'plane': + constrain_curvatures(cc_mean.output, curvature, 0.0, 0.0) + if surface_name == 'sphere': + # The sphere radius is 1. + # Mean curvature is 1/r + constrain_curvatures(cc_mean.output, curvature, 1.0, 1.0) + curvatures[curvature] = get_curvature_glyphs(surface_name, cc_mean, curvature, + precision=10, frequency_table=frequency_table, + nearest_integer=False) + + elev_bg = get_elevation_glyphs(surface_name, source, + precision=10, frequency_table=frequency_table, + nearest_integer=False) colors = vtkNamedColors() - colors.SetColor('ParaViewBlueGrayBkg', 84, 89, 109, 255) - colors.SetColor('ParaViewWarmGrayBkg', 98, 93, 90, 255) # Define viewport ranges [x_min, y_min, x_max, y_max] viewports = dict() @@ -162,309 +191,154 @@ def main(argv): renderers = list() contour_widgets = dict() elevation_widgets = dict() - # Set up the scalar bar properties. - scalar_bar_properties = ScalarBarProperties() # Position the source name according to its length. text_positions = get_text_positions(available_surfaces, justification=TextProperty.Justification.VTK_TEXT_LEFT, vertical_justification=TextProperty.VerticalJustification.VTK_TEXT_TOP, - width=0.45) - - text_property = vtkTextProperty(color=colors.GetColor3d('AliceBlue'), bold=True, italic=True, shadow=True, - font_size=16, - justification=TextProperty.Justification.VTK_TEXT_LEFT) + width=0.225) + + title_text_property = vtkTextProperty(color=colors.GetColor3d('AliceBlue'), bold=True, italic=True, shadow=True, + font_size=16, + justification=TextProperty.Justification.VTK_TEXT_LEFT) + label_text_property = vtkTextProperty(color=colors.GetColor3d('AliceBlue'), bold=False, italic=False, shadow=True, + font_size=12, + justification=TextProperty.Justification.VTK_TEXT_LEFT) text_actor = vtkTextActor(input=surface_name.title(), text_scale_mode=vtkTextActor.TEXT_SCALE_MODE_NONE, - text_property=text_property) + text_property=title_text_property) # Create the text representation. Used for positioning the text actor. text_representation = vtkTextRepresentation(enforce_normalized_viewport_bounds=True) - text_representation.position_coordinate.value = text_positions[surface.name]['p'] - text_representation.position2_coordinate.value = text_positions[surface.name]['p2'] + text_representation.position_coordinate.value = text_positions[surface_name]['p'] + text_representation.position2_coordinate.value = text_positions[surface_name]['p2'] text_widget = vtkTextWidget(representation=text_representation, text_actor=text_actor, interactor=iren, selectable=False) - first = True + # Scalar bar properties for the curvatures and elevation. + curv_sbps = dict() for k, v in curvatures.items(): - src_mapper = vtkPolyDataMapper(scalar_range=v['curv_scalar_range'], - lookup_table=v['curv_lut'], - scalar_mode=Mapper.ScalarMode().VTK_SCALAR_MODE_USE_CELL_DATA) + curv_sbp = ScalarBarProperties() + curv_sbp.title_text = v.glyph_type.replace('_', '\n') + '\n' + curv_sbp.number_of_labels = len(v.labels) + curv_sbp.lut = v.lut + curv_sbp.orientation = False + if curv_sbp.number_of_labels == 1: + curv_sbp.position_h = position_sbw_h(curv_sbp.number_of_labels, v.lut_max_bands + 3) + else: + curv_sbp.position_h = position_sbw_h(curv_sbp.number_of_labels, v.lut_max_bands) + curv_sbps[k] = curv_sbp + + elev_sbp = ScalarBarProperties() + elev_sbp.title_text = elev_bg.glyph_type.replace('_', '\n') + '\n' + elev_sbp.number_of_labels = len(elev_bg.labels) + elev_sbp.lut = elev_bg.lutr + elev_sbp.orientation = True + if elev_sbp.number_of_labels == 1: + elev_sbp.position_v = position_sbw_v(elev_sbp.number_of_labels, elev_bg.lut_max_bands) + else: + elev_sbp.position_v = position_sbw_v(elev_sbp.number_of_labels, elev_bg.lut_max_bands) + # Put it all together, + # Gaussian curvature in the left viewport, mean curvature in the right. + for k, v in curvatures.items(): + src_mapper = vtkPolyDataMapper(input_connection=v.bcf.output_port, + lookup_table=v.lut, + scalar_range=v.scalar_range, + scalar_mode=Mapper.ScalarMode.VTK_SCALAR_MODE_USE_CELL_DATA) src_actor = vtkActor(mapper=src_mapper) - v['bcf'] >> src_mapper # Create contour edges - edge_mapper = vtkPolyDataMapper( - resolve_coincident_topology=Mapper.ResolveCoincidentTopology.VTK_RESOLVE_POLYGON_OFFSET) + edge_mapper = vtkPolyDataMapper(input_data=v.bcf.contour_edges_output, + resolve_coincident_topology=Mapper.ResolveCoincidentTopology.VTK_RESOLVE_POLYGON_OFFSET) edge_actor = vtkActor(mapper=edge_mapper) edge_actor.property.color = colors.GetColor3d('Black') - v['bcf'].GetContourEdgesOutput() >> edge_mapper - glyph_mapper = vtkPolyDataMapper(scalar_range=v['elev_scalar_range'], - lookup_table=v['elev_lut'], - scalar_mode=Mapper.ScalarMode.VTK_SCALAR_MODE_USE_POINT_FIELD_DATA, + glyph_mapper = vtkPolyDataMapper(input_connection=elev_bg.glyphs.output_port, + lookup_table=elev_bg.lut1, + scalar_range=elev_bg.scalar_range, + color_mode=Mapper.ColorMode.VTK_COLOR_MODE_MAP_SCALARS, scalar_visibility=True, - color_mode=Mapper.ColorMode.VTK_COLOR_MODE_MAP_SCALARS) - glyph_mapper.SelectColorArray('Elevation') + scalar_mode=Mapper.ScalarMode.VTK_SCALAR_MODE_USE_POINT_FIELD_DATA) + glyph_mapper.SelectColorArray('Elevation') glyph_actor = vtkActor(mapper=glyph_mapper) - v['glyph'] >> glyph_mapper - renderer = vtkRenderer(background=colors.GetColor3d('ParaViewBlueGrayBkg')) + ren = vtkRenderer(background=colors.GetColor3d('ParaViewBlueGrayBkg')) - if surface_name in ['plane', 'sphere']: - scalar_bar_properties.position_h = {'p': (0.3, 0.1), 'p2': (0.35, 0.1)} - scalar_bar_properties.lut = curvatures[k]['curv_lut'] - scalar_bar_properties.orientation = False - scalar_bar_properties.title_text = k.replace('_', ' ') + '\n' - scalar_bar_properties.number_of_labels = curvatures[k]['elev_lut'].GetNumberOfTableValues() - contour_widgets[k] = make_scalar_bar_widget(scalar_bar_properties, text_property, text_property, renderer, iren) - if surface_name in ['plane', 'sphere']: - scalar_bar_properties.position_h = copy.deepcopy(scalar_bar_properties.default_h) + contour_widgets[k] = make_scalar_bar_widget(curv_sbps[k], title_text_property, + label_text_property, ren, iren) + contour_widgets[k].Off() - if first: - text_widget.default_renderer = renderer - first = False - renderer.SetViewport(*viewports[k]) - renderer.AddActor(src_actor) - renderer.AddActor(edge_actor) - renderer.AddActor(glyph_actor) - - # Now for the elevation, it is the same for both surface actors. - # So we will just make it for the maan curvature. - if k == 'Mean_Curvature': - if surface_name == 'plane': - scalar_bar_properties.position_v = {'p': (0.85, 0.5), 'p2': (0.1, 0.1)} - scalar_bar_properties.lut = curvatures[k]['elev_lut'] - scalar_bar_properties.orientation = True - scalar_bar_properties.title_text = 'Elevation\n' - scalar_bar_properties.number_of_labels = curvatures[k]['elev_labels'] - elevation_widgets[k] = make_scalar_bar_widget(scalar_bar_properties, text_property, text_property, renderer, - iren) - if surface_name == 'plane': - scalar_bar_properties.position_v = copy.deepcopy(scalar_bar_properties.default_v) - - contour_widgets[k].default_renderer = renderer - if k == 'Mean_Curvature': - elevation_widgets[k].default_renderer = renderer - - renderers.append(renderer) - - for renderer in renderers: - ren_win.AddRenderer(renderer) + ren.SetViewport(*viewports[k]) + ren.AddActor(src_actor) + ren.AddActor(edge_actor) + ren.AddActor(glyph_actor) + + elevation_widgets[k] = make_scalar_bar_widget(elev_sbp, title_text_property, + label_text_property, ren, iren) + elevation_widgets[k].Off() + renderers.append(ren) + ren_win.AddRenderer(ren) + + ren.ResetCamera() + + adjust_camera_parameters(surface_name, ren) + + # Share the camera. + renderers[1].active_camera = renderers[0].active_camera + # Set the default renderer for the text widget. + text_widget.default_renderer = renderers[0] + + # Switch on the contour and elevation widgets. for k in curvatures.keys(): if k == 'Gauss_Curvature': contour_widgets[k].On() + elevation_widgets[k].On() else: contour_widgets[k].On() elevation_widgets[k].On() text_widget.On() - if use_camera_omw: - cow = vtkCameraOrientationWidget(parent_renderer=renderers[0]) - # Enable the widget. - cow.On() - else: + # Important: The interactor must be set prior to enabling the widgets. + if use_omw: rgb = [0.0] * 4 colors.GetColor("Carrot", rgb) rgb = tuple(rgb[:3]) widget = vtkOrientationMarkerWidget(orientation_marker=vtkAxesActor(), - interactor=iren, default_renderer=renderers[1], - outline_color=rgb, viewport=(0.7, 0.8, 0.9, 1.0), zoom=1.5, enabled=True, + interactor=iren, default_renderer=ren, + outline_color=rgb, viewport=(0.8, 0.8, 1.0, 1.0), zoom=1.0, enabled=True, interactive=True) + else: + cow = vtkCameraOrientationWidget(parent_renderer=ren) + # Enable the widget. + cow.On() - camera = None - for i in range(0, len(renderers)): - if i == 0: - camera = renderers[i].active_camera - camera.Elevation(60) - # This moves the window center slightly to ensure that - # the whole surface is not obscured by the scalar bars. - camera.window_center = (0.0, -0.15) - else: - renderers[i].active_camera = camera - renderers[i].ResetCamera() - - if surface_name == 'plane': - renderers[0].active_camera.Zoom(0.8) ren_win.Render() - iren.Start() -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. - - :param source: 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): - """ - Extract the topological neighbors for point. - - :param pt_id: The point id. - :return: The neighbour 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: First point. - :param pt_id_b: Second point. - :return: The distance. - """ - 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.point_data.SetActiveScalars(curvature_name) - np_source = dsa.WrapDataObject(source) - curvatures = np_source.PointData[curvature_name] - - # Get the boundary point IDs. - array_name = 'ids' - id_filter = vtkGenerateIds(point_ids=True, cell_ids=False, - point_ids_array_name=array_name, - cell_ids_array_name=array_name) - - edges = vtkFeatureEdges(boundary_edges=True, manifold_edges=False, - non_manifold_edges=False, feature_edges=False) - - (source >> id_filter >> edges).update() - - edge_array = edges.output.GetPointData().GetArray(array_name) - boundary_ids = [] - for i in range(edges.output.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) - curv = numpy_support.numpy_to_vtk(num_array=curvatures.ravel(), - deep=True, - array_type=VTK_DOUBLE) - curv.name = curvature_name - source.point_data.RemoveArray(curvature_name) - source.point_data.AddArray(curv) - source.point_data.active_scalars = 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.point_data.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) - curv = numpy_support.numpy_to_vtk(num_array=curvatures.ravel(), - deep=True, - array_type=VTK_DOUBLE) - curv.name = curvature_name - source.point_data.RemoveArray(curvature_name) - source.point_data.AddArray(curv) - source.point_data.active_scalars = curvature_name - - -def get_source(source, available_surfaces): +def generate_elevations(src): """ - - :param source: The name of the source. - :param available_surfaces: The surfaces - :return: + Generate elevations over the surface. + :param: src - the vtkPolyData source. + :return: - vtkPolyData source with elevations. """ - surface = source.lower() - if surface not in available_surfaces: - return None - elif surface == 'hills': - return get_hills() - elif surface == 'parametric torus': - return get_parametric_torus() - elif surface == 'plane': - return get_plane() - elif surface == 'random hills': - return get_parametric_hills() - elif surface == 'sphere': - return get_sphere() - elif surface == 'torus': - return get_torus() - return None + bounds = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0] + src.GetBounds(bounds) + if abs(bounds[2]) < 1.0e-8 and abs(bounds[3]) < 1.0e-8: + bounds[3] = bounds[2] + 1 + elev_filter = vtkElevationFilter(input_data=src, + low_point=(0, bounds[2], 0), + high_point=(0, bounds[3], 0), + scalar_range=(bounds[2], bounds[3])) + elev_filter.update() + return elev_filter.GetPolyDataOutput() def get_hills(): - """ - Create four hills on a plane. - This will have regions of negative, zero and positive Gaussian curvatures. - - :return: - """ + # Create four hills on a plane. + # This will have regions of negative, zero and positive Gaussian curvatures. x_res = 50 y_res = 50 @@ -476,75 +350,38 @@ def get_hills(): dy = (y_max - y_min) / (x_res - 1) # Make a grid. - # We define the parameters for the hills here. - # There are four hills. + 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(points=points) + + # Triangulate the grid. + delaunay = vtkDelaunay2D(input_data=plane) + + polydata = delaunay.update().output + + elevation = vtkDoubleArray(number_of_tuples=points.number_of_points) + + # 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]] - # Three hills. - # hd = [[-2.5, -2.5, 2.5, 6.5, 3.5], [-5.0, 5, 2.5, 3.0, 3], [5.0, -2.5, 1.5, 1.5, 2.5]] - # Two hills. - # hd = [[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.number_of_points): + 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]) - # Generate the points. - y = y_min - xyz = list() - for r in range(x_res): - x = x_min - for c in range(y_res): - z = 0 - for j in range(0, len(hd)): - xx[0] = (x - hd[j][0] / hd[j][2]) ** 2.0 - xx[1] = (y - hd[j][1] / hd[j][3]) ** 2.0 - z += hd[j][4] * math.exp(-(xx[0] + xx[1]) / 2.0) - xyz.append([x, y, z]) - x += dx - y += dy - np_xyz = np.array(xyz) - points = vtkPoints(data=numpy_support.numpy_to_vtk(np_xyz)) - # print(points) - - # Triangulate each quad. - triangles = vtkCellArray() - t_num = 0 - for r in range(0, x_res - 1): - for c in range(0, y_res - 1): - indices = list() - - # We select the index of each point so that - # the ordering is counterclockwise. - - rc = r * y_res - # Triangle 1. - indices.append(rc + c) - indices.append(indices[0] + 1) - indices.append(rc + y_res + c) - # print(f't_num: {t_num} :- {indices[0]}, {indices[1]}, {indices[2]}') - triangle = vtkTriangle() - for i in range(0, len(indices)): - triangle.GetPointIds().SetId(i, indices[i]) - triangles.InsertNextCell(triangle) - indices.clear() - t_num += 1 - - # Triangle 2. - indices.append(rc + c + 1) - indices.append(rc + y_res + c + 1) - indices.append(rc + y_res + c) - # print(f't_num: {t_num} :- {indices[0]}, {indices[1]}, {indices[2]}') - triangle = vtkTriangle() - for i in range(0, len(indices)): - triangle.GetPointIds().SetId(i, indices[i]) - triangles.InsertNextCell(triangle) - indices.clear() - t_num += 1 - - poly_data = vtkPolyData(points=points, polys=triangles) - - textures = vtkFloatArray(name='Textures', number_of_components=2, number_of_tuples=2 * poly_data.number_of_points) - poly_data.GetPointData().SetTCoords(textures) + textures = vtkFloatArray(name='Textures', number_of_components=2, number_of_tuples=2 * polydata.number_of_points) for i in range(0, x_res): tc = [i / (x_res - 1.0), 0.0] @@ -553,68 +390,99 @@ def get_hills(): tc[1] = j / (y_res - 1.0) textures.SetTuple(i * y_res + j, tc) - bounds = poly_data.bounds - elevation_filter = vtkElevationFilter(low_point=(0.0, 0.0, bounds[4]), high_point=(0.0, 0.0, bounds[5])) + polydata.GetPointData().SetScalars(elevation) + polydata.GetPointData().GetScalars().name = 'Elevation' + polydata.GetPointData().TCoords = textures normals = vtkPolyDataNormals(feature_angle=30, splitting=False) tr = vtkTransform() - tr.RotateX(-90.0) + tr.RotateX(-90) + tf = vtkTransformFilter(transform=tr) - return poly_data >> elevation_filter >> normals >> tf + return (polydata >> normals >> tf).update().output def get_parametric_hills(): - fn = vtkParametricRandomHills(random_seed=1, number_of_hills=30) + """ + Make a parametric hills surface as the source. + :return: vtkPolyData with normal and scalar data. + """ + random_seed = 1 + number_of_hills = 30 + # If you want a plane + # hill_amplitude=0 + fn = vtkParametricRandomHills(random_seed=random_seed, number_of_hills=number_of_hills) fn.AllowRandomGenerationOn() - source = vtkParametricFunctionSource(parametric_function=fn, u_resolution=50, v_resolution=50, - scalar_mode=vtkParametricFunctionSource.SCALAR_Z) + u_resolution = 50 + v_resolution = 50 + source = vtkParametricFunctionSource(parametric_function=fn, + u_resolution=u_resolution, v_resolution=v_resolution, + generate_texture_coordinates=True) source.SetScalarModeToZ() - src = source.update().output - + source.update() # Rename the scalars to 'Elevation' since we are using the Z-scalars as elevations. - src.point_data.scalars.SetName('Elevation') + source.output.GetPointData().GetScalars().SetName('Elevation') + + # Build the tangents. + tangents = vtkPolyDataTangents() - transform = vtkTransform() - transform.Translate(0.0, 5.0, 15.0) - transform.RotateX(-90.0) - transform_filter = vtkTransformFilter(transform=transform) + tr = vtkTransform() + tr.Translate(0.0, 5.0, 15.0) + tr.RotateX(-90.0) + + tf = vtkTransformFilter(transform=tr) - return src >> transform_filter + return (source >> tangents >> tf).update().output def get_parametric_torus(): - fn = vtkParametricTorus(ring_radius=5, cross_section_radius=2) + """ + Make a parametric torus as the source. + :return: vtkPolyData with normal and scalar data. + """ - source = vtkParametricFunctionSource(parametric_function=fn, u_resolution=50, v_resolution=50, - scalar_mode=vtkParametricFunctionSource.SCALAR_Z) - src = source.update().output + fn = vtkParametricTorus(ring_radius=5, cross_section_radius=2) + source = vtkParametricFunctionSource(parametric_function=fn, + u_resolution=50, v_resolution=50, + generate_texture_coordinates=True) + source.SetScalarModeToZ() + source.update() # Rename the scalars to 'Elevation' since we are using the Z-scalars as elevations. - src.point_data.scalars.SetName('Elevation') + source.output.GetPointData().GetScalars().SetName('Elevation') + + # Build the tangents. + tangents = vtkPolyDataTangents() + + tr = vtkTransform() + # transform.Translate(0.0, 0.0, 0.0) + tr.RotateX(-90.0) - transform = vtkTransform() - transform.Translate(0.0, 0.0, 0.0) - transform.RotateX(-90.0) - transform_filter = vtkTransformFilter(transform=transform) + t = vtkTransformFilter(transform=tr) - return src >> transform_filter + return (source >> tangents >> t).update().output def get_plane(): - source = vtkPlaneSource(origin=(-10.0, -10.0, 0.0), point1=(10.0, -10.0, 0.0), point2=(-10.0, 10.0, 0.0), - x_resolution=5, y_resolution=5) - src = source.update().output + """ + Make a plane as the source. + :return: vtkPolyData with normal and scalar data. + """ + + source = vtkPlaneSource(origin=(-10.0, -10.0, 0.0), + point2=(-10.0, 10.0, 0.0), point1=(10.0, -10.0, 0.0), + x_resolution=10, y_resolution=10) - transform = vtkTransform() - transform.Translate(0.0, 0.0, 0.0) - transform.RotateX(-90.0) + tr = vtkTransform() + tr.Translate(0.0, 0.0, 0.0) + tr.RotateX(-90.0) - transform_filter = vtkTransformFilter(transform=transform) + tf = vtkTransformFilter(transform=tr) - # We have an m x n array of quadrilaterals arranged as a regular tiling in a + # 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 = vtkTriangleFilter() @@ -623,25 +491,25 @@ def get_plane(): # are coincident, or very close cleaner = vtkCleanPolyData(tolerance=0.005) - elev_filter = vtkElevationFilter(low_point=(0, 0, 0), high_point=(0, 0, 1), scalar_range=(0, 0)) - - return src >> transform_filter >> tri >> cleaner >> elev_filter + return (source >> tf >> tri >> cleaner).update().output def get_sphere(): - source = vtkSphereSource(center=(0.0, 0.0, 0.0), radius=10.0, theta_resolution=32, phi_resolution=32) - src = source.update().output + source = vtkSphereSource(center=(0.0, 0.0, 0.0), radius=1.0, + theta_resolution=32, phi_resolution=32) - elev_filter = vtkElevationFilter(low_point=(0, src.bounds[2], 0), high_point=(0, src.bounds[3], 0), - scalar_range=(src.bounds[2], src.bounds[3])) - - return src >> elev_filter + return source.update().output def get_torus(): - source = vtkSuperquadricSource(center=(0.0, 0.0, 0.0), scale=(1.0, 1.0, 1.0), phi_resolution=64, - theta_resolution=64, theta_roundness=1, thickness=0.5, size=10, toroidal=True) - src = source.update().output + """ + Make a torus as the source. + :return: vtkPolyData with normal and scalar data. + """ + source = vtkSuperquadricSource(center=(0.0, 0.0, 0.0), scale=(1.0, 1.0, 1.0), + phi_resolution=64, + theta_resolution=64, theta_roundness=1, + thickness=0.5, size=10, toroidal=True) # The quadric is made of strips, so pass it through a triangle filter as # the curvature filter only operates on polys @@ -652,67 +520,26 @@ def get_torus(): # are coincident, or very close cleaner = vtkCleanPolyData(tolerance=0.005) - elev_filter = vtkElevationFilter(low_point=(0, src.bounds[2], 0), high_point=(0, src.bounds[3], 0), - scalar_range=(src.bounds[2], src.bounds[3])) - - return src >> tri >> cleaner >> elev_filter + return (source >> tri >> cleaner).update().output -def get_categorical_lut(): - """ - Make a lookup table using vtkColorSeries. - :return: An indexed (categorical) lookup table. - """ - color_series = vtkColorSeries(color_scheme=vtkColorSeries.BREWER_QUALITATIVE_SET3) - # Make the lookup table. - lut = vtkLookupTable() - color_series.BuildLookupTable(lut, color_series.CATEGORICAL) - lut.nan_color = (0, 1, 0, 1) - return lut - - -def get_ordinal_lut(): - """ - Make a lookup table using vtkColorSeries. - :return: An ordinal (not indexed) lookup table. - """ - color_series = vtkColorSeries(color_scheme=vtkColorSeries.BREWER_DIVERGING_BROWN_BLUE_GREEN_11) - # Make the lookup table. - lut = vtkLookupTable() - color_series.BuildLookupTable(lut, color_series.ORDINAL) - lut.nan_color = (0, 1, 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: The lookup table. - """ - ctf = vtkColorTransferFunction(color_space=ColorTransferFunction.ColorSpace.VTK_CTF_DIVERGING) - # 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 get_source(source): + surface = source.lower() + if surface == 'hills': + return get_hills() + elif surface == 'parametric hills': + return get_parametric_hills() + elif surface == 'parametric torus': + return get_parametric_torus() + elif surface == 'plane': + return generate_elevations(get_plane()) + elif surface == 'sphere': + return generate_elevations(get_sphere()) + elif surface == 'torus': + return generate_elevations(get_torus()) + print('The surface is not available.') + print('Using parametric hills instead.') + return get_parametric_hills() def reverse_lut(lut): @@ -723,7 +550,7 @@ def reverse_lut(lut): """ lutr = vtkLookupTable() lutr.DeepCopy(lut) - t = lut.GetNumberOfTableValues() - 1 + t = lut.number_of_table_values - 1 rev_range = reversed(list(range(t + 1))) for i in rev_range: rgba = [0.0] * 3 @@ -734,68 +561,52 @@ def reverse_lut(lut): t = lut.number_of_annotated_values - 1 rev_range = reversed(list(range(t + 1))) for i in rev_range: - lutr.annotation = (t - i, lut.GetAnnotation(i)) + lutr.SetAnnotation(t - i, lut.GetAnnotation(i)) return lutr -def get_glyphs(surface, arrow_scale=None, scale_factor=None, reverse_normals=False): +def get_glyphs(src, scale_factor=1.0, reverse_normals=False): """ - Glyph the surface. + Glyph the normals on the surface. - :param surface: The surface to glyph. - :param arrow_scale: Scaling for the arrows, default is [1, 1, 1]. - :param scale_factor: The scaling factor for the arrow, default is 1.0. - :param reverse_normals: If True the normals on the surface are reversed. - :return: The glyph filter. - """ - name = surface.name - source = surface.source - - if arrow_scale is None: - arrow_scale = [1, 1, 1] - # The length of the arrow glyph. - if scale_factor is None: - scale_factor = 1.0 - - # Choose a random subset of points. - if name == 'plane': - mask_pts = vtkMaskPoints(on_ratio=1, random_mode=True) - else: - mask_pts = vtkMaskPoints(on_ratio=5, random_mode=True) + You may need to adjust the parameters for mask_pts, arrow and glyph for a + nice appearance. - # Sometimes the contouring algorithm can create a volume whose gradient - # vector and ordering of the polygon (using the right hand rule) are - # inconsistent. vtkReverseSense cures this problem. + :param: src - the surface to glyph. + :param: reverse_normals - if True the normals on the surface are reversed. + :return: The glyph object. + + """ if 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 = vtkReverseSense(reverse_cells=True, reverse_normals=True) - source >> reverse >> mask_pts + # Choose a random subset of points. + mask_pts = vtkMaskPoints(on_ratio=5, random_mode=True) + src >> reverse >> mask_pts else: - source >> mask_pts + # Choose a random subset of points. + mask_pts = vtkMaskPoints(on_ratio=5, random_mode=True) + src >> mask_pts - # Source for the glyph filter. - arrow = vtkArrowSource(shaft_resolution=16, shaft_radius=0.03, tip_resolution=16, tip_length=0.3, tip_radius=0.1) - # Scale the arrow. - tr = vtkTransform() - tr.Scale(arrow_scale) - tf = vtkTransformFilter(transform=tr) - - p = (arrow >> tf).update().output + # Source for the glyph filter + arrow = vtkArrowSource(tip_resolution=16, tip_length=0.3, tip_radius=0.1) - glyph = vtkGlyph3D(source_data=p, scale_factor=scale_factor, - vector_mode=Glyph3D.VectorMode.VTK_USE_NORMAL, - color_mode=Glyph3D.ColorMode.VTK_COLOR_BY_VECTOR, - scale_mode=Glyph3D.ScaleMode.VTK_SCALE_BY_VECTOR - ) - glyph.OrientOn() + # glyph = vtkGlyph3D() + glyphs = vtkGlyph3D(source_connection=arrow.output_port, + input_connection=mask_pts.output_port, + scaling=True, scale_mode=Glyph3D.ScaleMode.VTK_SCALE_BY_VECTOR, + scale_factor=scale_factor, orient=True, clamping=False, + vector_mode=Glyph3D.VectorMode.VTK_USE_NORMAL, + color_mode=Glyph3D.ColorMode.VTK_COLOR_BY_VECTOR) + return glyphs - return mask_pts >> glyph - -def get_bands(d_r, number_of_bands, precision=2, nearest_integer=False): +def get_bands(scalar_range, 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. + Divide a range into bands + :param: scalar_range - [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. @@ -806,9 +617,9 @@ def get_bands(d_r, number_of_bands, precision=2, nearest_integer=False): prec = 14 bands = dict() - if (d_r[1] < d_r[0]) or (number_of_bands <= 0): + if (scalar_range[1] < scalar_range[0]) or (number_of_bands <= 0): return bands - x = list(d_r) + x = list(scalar_range) if nearest_integer: x[0] = math.floor(x[0]) x[1] = math.ceil(x[1]) @@ -825,39 +636,37 @@ 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): +def get_custom_bands(scalar_range, number_of_bands, my_bands): """ Divide a range into custom bands. - You need to specify each band as a 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]...] + You need to specify each band as a 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: scalar_range - [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): + if (scalar_range[1] < scalar_range[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]: + if my_bands[idx][1] > scalar_range[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]: + if my_bands[idx][1] > scalar_range[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[idx_min][0] = scalar_range[0] + x[idx_max][1] = scalar_range[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]] @@ -886,16 +695,155 @@ def get_frequencies(bands, src): return freq +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. + + :param source: 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): + """ + Extract the topological neighbors for point. + + :param pt_id: The point id. + :return: The neighbour 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: First point. + :param pt_id_b: Second point. + :return: The distance. + """ + 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.point_data.SetActiveScalars(curvature_name) + np_source = dsa.WrapDataObject(source) + curvatures = np_source.PointData[curvature_name] + + # Get the boundary point IDs. + array_name = 'ids' + id_filter = vtkGenerateIds(point_ids=True, cell_ids=False, + point_ids_array_name=array_name, + cell_ids_array_name=array_name) + + edges = vtkFeatureEdges(boundary_edges=True, manifold_edges=False, + non_manifold_edges=False, feature_edges=False) + + (source >> id_filter >> edges).update() + + edge_array = edges.output.GetPointData().GetArray(array_name) + boundary_ids = [] + for i in range(edges.output.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) + curv = numpy_support.numpy_to_vtk(num_array=curvatures.ravel(), + deep=True, + array_type=VTK_DOUBLE) + curv.name = curvature_name + source.point_data.RemoveArray(curvature_name) + source.point_data.AddArray(curv) + source.point_data.active_scalars = 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.point_data.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) + curv = numpy_support.numpy_to_vtk(num_array=curvatures.ravel(), + deep=True, + array_type=VTK_DOUBLE) + curv.name = curvature_name + source.point_data.RemoveArray(curvature_name) + source.point_data.AddArray(curv) + source.point_data.active_scalars = curvature_name + + 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 bands: The dictionary containing the bands. :param freq: The frequency dictionary. :return: Adjusted bands and frequencies. """ - # Get the indices of the first and last non-zero frequency elements. + # Get the indices of the first and last non-zero elements. first = 0 for k, v in freq.items(): if v != 0: @@ -927,37 +875,120 @@ def adjust_ranges(bands, freq): return adj_bands, adj_freq -def print_bands_frequencies(curvature, bands, freq, precision=2): - prec = abs(precision) - if prec > 14: - prec = 14 +def get_curvature_glyphs(surface_name, source, curvature, + precision, frequency_table=False, nearest_integer=False): + """ + Get curvature glyphs and the corresponding banded polydata filter for the surface. + :param: surface_name - the name of the surface. + :param: src - the polydata surface to glyph. + :param: curvature - The curvature type: Gauss_Curvature or Mean_Curvature + :param: precision - the precision level. + :param: frequency_table - True if you want a frequency table printed. + :param: nearest_integer - If True use the nearest integer when generating the bands. + :return: A dataclass holding glyphs, bcf, lut, lutr, lut1, lut1r, scalar_range, labels - if len(bands) != len(freq): - print('Bands and Frequencies must be the same size.') - return - s = f'Bands & Frequencies:\n{" ".join(curvature.lower().replace("_", " ").split()).title()}\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:{width}.{prec}f}]: {freq[k]:8d}\n' - else: - s += f'{q:{width}.{prec}f}, ' - width = 3 * width + 13 - s += f'{"Total":{width}s}{total:8d}\n' - print(s) + """ + # The length of the normal arrow glyphs. + scale_factor = 1.0 + if surface_name == 'hills': + scale_factor = 0.5 + elif surface_name == 'sphere': + scale_factor = 0.2 + + source.output.GetPointData().SetActiveScalars(curvature) + scalar_range = source.output.point_data.GetScalars(curvature).range + + max_lut_bands = 0 + + color_series = vtkColorSeries() + color_series = vtkColorSeries(color_series.BREWER_QUALITATIVE_SET3) + + if surface_name == 'hills': + color_series.color_scheme = color_series.BREWER_DIVERGING_SPECTRAL_9 + max_lut_bands = 9 + elif surface_name == 'parametric hills': + color_series.color_scheme = color_series.BREWER_DIVERGING_SPECTRAL_7 + max_lut_bands = 7 + elif surface_name in ['plane', 'sphere']: + color_series.color_scheme = color_series.BREWER_SEQUENTIAL_YELLOW_ORANGE_BROWN_3 + max_lut_bands = 3 + else: + color_series.color_scheme = color_series.BREWER_DIVERGING_SPECTRAL_6 + max_lut_bands = 6 + lut = vtkLookupTable() + color_series.BuildLookupTable(lut, color_series.CATEGORICAL) + lut.SetNanColor(0, 0, 0, 1) + lut.SetTableRange(scalar_range) + + if surface_name in ['plane', 'sphere']: + lut.number_of_table_values = 1 + + lut1 = vtkLookupTable() + color_series.BuildLookupTable(lut1, color_series.ORDINAL) + lut1.SetNanColor(0, 0, 0, 1) + lut1.SetTableRange(scalar_range) + + number_of_bands = lut.number_of_table_values + lut1.number_of_table_values = number_of_bands + + bands = get_bands(scalar_range, number_of_bands, precision, nearest_integer) + + if curvature == 'Gauss_Curvature': + if surface_name == 'parametric hills': + # These are my custom bands. + # Generated by first running: + # bands = get_bands(scalar_range_curvatures, number_of_bands, False) + # then: + # freq = frequencies(bands, src) + # print_bands_frequencies(bands, freq) + # Finally using the output to create this table: + # my_bands = [ + # [-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.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) + # and some spherical regions. + my_bands = [ + [-0.630, -0.190], [-0.190, -0.043], [-0.043, 0.0452], [0.0452, 0.0746], + [0.0746, 0.104], [0.104, 0.251], [0.251, 1.131]] + # Comment this out if you want to see how allocating + # equally spaced bands works. + bands = get_custom_bands(scalar_range, number_of_bands, my_bands) + elif surface_name == 'hills': + my_bands = [ + [-2.104, -0.15], [-0.15, -0.1], [-0.1, -0.05], + [-0.05, -0.02], [-0.02, -0.005], [-0.005, -0.0005], + [-0.0005, 0.0005], [0.0005, 0.09], [0.09, 4.972]] + # Comment this out if you want to see how allocating + # equally spaced bands works. + bands = get_custom_bands(scalar_range, number_of_bands, my_bands) + + # Adjust the number of table values and scalar range. + scalar_range = (bands[0][0], bands[len(bands) - 1][2]) + lut.TableRange = scalar_range + lut.number_of_table_values = len(bands) + lut1.TableRange = scalar_range + lut1.number_of_table_values = len(bands) + + if frequency_table: + print(f'{surface_name.title()} {curvature}') + # The number of scalars in each band. + freq = get_frequencies(bands, source.output) + bands, freq = adjust_ranges(bands, freq) + print_bands_frequencies(bands, freq) + + scalar_range = (bands[0][0], bands[len(bands) - 1][2]) + lut.TableRange = scalar_range + lut.number_of_table_values = len(bands) + lut1.TableRange = scalar_range + lut1.number_of_table_values = len(bands) -def annotate_lut(lut, bands): # We will use the midpoint of the band as the label. labels = [] - keys = sorted(bands.keys()) - for k in keys: - labels.append(f'{bands[k][1]:4.2f}') + for k in bands: + labels.append('{:4.2f}'.format(bands[k][1])) # Annotate values = vtkVariantArray() @@ -966,204 +997,160 @@ def annotate_lut(lut, bands): for i in range(values.GetNumberOfTuples()): lut.SetAnnotation(i, values.GetValue(i).ToString()) - -def generate_gaussian_curvatures(surface, needs_adjusting, frequency_table=False): + # Create the contour bands. + # We will use an indexed lookup table. + bcf = vtkBandedPolyDataContourFilter(input_data=source.output, + scalar_mode=BandedPolyDataContourFilter.ScalarMode.VTK_SCALAR_MODE_INDEX, + generate_contour_edges=True) + # Use either the minimum or maximum value for each band. + for i in range(len(bands)): + bcf.SetValue(i, bands[i][2]) + + glyphs = get_glyphs(source, scale_factor, reverse_normals=False) + + BandedGlyphs = namedtuple('BandedGlyphs', + ['glyph_type', 'glyphs', 'bcf', + 'lut_max_bands', 'lut', 'lutr', + 'lut1_max_bands', 'lut1', 'lut1r', + 'scalar_range', 'labels']) + bg = BandedGlyphs + bg.glyph_type = curvature + bg.glyphs = glyphs + bg.bcf = bcf + bg.lut_max_bands = max_lut_bands + bg.lut = lut + bg.lutr = reverse_lut(lut) + bg.lut1_max_bands = max_lut_bands + bg.lut1 = lut1 + bg.lut1r = reverse_lut(lut1) + bg.scalar_range = scalar_range + bg.labels = labels + + return bg + + +def get_elevation_glyphs(surface_name, source, + precision, frequency_table=False, nearest_integer=False): """ - Generate the filters for calculating Gaussian curvatures on the surface. + Get elevation glyphs and the corresponding banded polydata filter for the surface. + :param: surface_name - the name of the surface. + :param: src - the polydata surface to glyph. + :param: precision - the precision level. + :param: frequency_table - If true, display a frequency table corresponding to the bands. + :param: nearest_integer - If true, use the nearest integer when generating the bands. + :return: A dataclass holding glyphs, bcf, lut, lutr, lut1, lut1r scalar_range, labels - :param surface: The surface. - :param needs_adjusting: Surfaces whose curvatures need to be adjusted along the edges of the surface or constrained. - :param frequency_table: True if a frequency table is to be displayed. - :return: Return the filters, scalar ranges of curvatures and elevation along with the lookup tables. """ - name = surface.name - source = surface.source - curvature = 'Gauss_Curvature' - - curvatures = vtkCurvatures(curvature_type=Curvatures.CurvatureType().VTK_CURVATURE_GAUSS) - p = (source >> curvatures).update().output - - if name in needs_adjusting: - adjust_edge_curvatures(p, curvature) - if name == 'plane': - constrain_curvatures(p, curvature, 0.0, 0.0) - if name == 'sphere': - # Gaussian curvature is 1/r^2. - radius = 10 - gauss_curvature = 1.0 / radius ** 2 - constrain_curvatures(p, curvature, gauss_curvature, gauss_curvature) - - p.GetPointData().SetActiveScalars('Elevation') - elev_scalar_range = p.GetPointData().GetScalars('Elevation').range - # elev_lut = get_ordinal_lut() - # elev_labels = elev_lut.GetNumberOfTableValues(); - elev_lut = get_diverging_lut() - elev_labels = 5 - if name == 'plane': - elev_labels = 1 - elev_lut.SetTableRange(elev_scalar_range) - - p.GetPointData().SetActiveScalars(curvature) - curv_scalar_range = p.GetPointData().GetScalars(curvature).range - curv_lut = get_categorical_lut() - curv_lut.SetTableRange(curv_scalar_range) - curv_num_bands = curv_lut.GetNumberOfTableValues() - curv_bands = get_bands(curv_scalar_range, number_of_bands=curv_num_bands, precision=10, nearest_integer=False) - - if name == 'random hills': - # These are my custom bands. - # Generated by first running: - # bands = get_bands(curv_scalar_range, number_of_bands=curv_num_bands, - # precision=2, nearest_integer=False) - # then: - # freq = frequencies(bands, curvatures_output) - # print_bands_frequencies(curvature, bands, freq) - # Finally using the output to create this table: - # my_bands = [ - # [-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.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) - # and some spherical regions. - my_bands = [ - [-0.630, -0.190], [-0.190, -0.043], [-0.043, 0.0452], [0.0452, 0.0746], - [0.0746, 0.104], [0.104, 0.251], [0.251, 1.131]] - # Comment this out if you want to see how allocating - # equally spaced bands works. - curv_bands = get_custom_bands(curv_scalar_range, number_of_bands=curv_num_bands, my_bands=my_bands) - # Adjust the number of table values - curv_lut.SetNumberOfTableValues(len(curv_bands)) - if name == 'hills': - my_bands = [ - [-2.104, -0.15], [-0.15, -0.1], [-0.1, -0.05], - [-0.05, -0.02], [-0.02, -0.005], [-0.005, -0.0005], - [-0.0005, 0.0005], [0.0005, 0.09], [0.09, 4.972]] - # Comment this out if you want to see how allocating - # equally spaced bands works. - curv_bands = get_custom_bands(curv_scalar_range, number_of_bands=curv_num_bands, my_bands=my_bands) - # Adjust the number of table values - curv_lut.SetNumberOfTableValues(len(curv_bands)) + # The length of the normal arrow glyphs. + scale_factor = 1.0 + if surface_name == 'hills': + scale_factor = 0.5 + elif surface_name == 'sphere': + scale_factor = 0.2 - curv_freq = get_frequencies(curv_bands, p) - curv_bands, curv_freq = adjust_ranges(curv_bands, curv_freq) - if frequency_table: - # Let's do a frequency table with the number of scalars in each band. - print_bands_frequencies(curvature, curv_bands, curv_freq) + source.point_data.active_scalars = 'Elevation' + scalar_range = source.point_data.GetScalars('Elevation').range - curv_lut.SetTableRange(curv_scalar_range) - curv_lut.SetNumberOfTableValues(len(curv_bands)) + max_lut_bands = 0 - annotate_lut(curv_lut, curv_bands) + color_series = vtkColorSeries() + if surface_name in ['hills', 'parametric hills']: + color_series.color_scheme = color_series.BREWER_DIVERGING_BROWN_BLUE_GREEN_8 + max_lut_bands = 8 + else: + color_series.color_scheme = color_series.BREWER_DIVERGING_BROWN_BLUE_GREEN_6 + max_lut_bands = 6 - # Create the contour bands. - # We will use an indexed lookup table. - bcf = vtkBandedPolyDataContourFilter(input_data=p, - scalar_mode=BandedPolyDataContourFilter.ScalarMode.VTK_SCALAR_MODE_INDEX, - generate_contour_edges=True) + lut = vtkLookupTable() + color_series.BuildLookupTable(lut, color_series.CATEGORICAL) + lut.SetNanColor(0, 0, 0, 1) + lut.SetTableRange(scalar_range) - # Use either the minimum or maximum value for each band. - for k in curv_bands: - bcf.SetValue(k, curv_bands[k][2]) + if surface_name in ['plane']: + lut.number_of_table_values = 1 - # Generate the glyphs on the original surface. - arrow_scale = [1, 1, 1] - if name == 'plane': - arrow_scale = [1.25, 1.25, 1.25] - scale_factor = 1.0 - if name == 'hills': - scale_factor = 0.5 - if name == 'sphere': - scale_factor = 2.0 + lut1 = vtkLookupTable() + color_series.BuildLookupTable(lut1, color_series.ORDINAL) + lut1.SetNanColor(0, 0, 0, 1) + lut1.SetTableRange(scalar_range) - glyph = get_glyphs(surface, arrow_scale=arrow_scale, scale_factor=scale_factor, reverse_normals=False) + number_of_bands = lut.number_of_table_values + lut1.number_of_table_values = number_of_bands - return {'bcf': bcf, 'glyph': glyph, 'curv_scalar_range': curv_scalar_range, - 'elev_scalar_range': elev_scalar_range, 'curv_lut': curv_lut, - 'elev_lut': elev_lut, 'elev_labels': elev_labels} + bands = get_bands(scalar_range, number_of_bands, precision, nearest_integer) + if surface_name == 'parametric hills': + # These are my custom bands. + # Generated by first running: + # bands = get_bands(scalar_range, number_of_bands, precision, False) + # then: + # freq = get_frequencies(bands, source) + # print_bands_frequencies(bands, freq) + # Finally using the output to create this table: + my_bands = [ + [0, 1.0], [1.0, 2.0], [2.0, 3.0], + [3.0, 4.0], [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 = get_custom_bands(scalar_range, number_of_bands, my_bands) -def generate_mean_curvatures(surface, needs_adjusting, frequency_table=False): - """ - Generate the filters for calculating mean curvatures on the surface. + # Adjust the number of table values and scalar range. + scalar_range = (bands[0][0], bands[len(bands) - 1][2]) + lut.TableRange = scalar_range + lut.number_of_table_values = len(bands) + lut1.TableRange = scalar_range + lut1.number_of_table_values = len(bands) - :param surface: The surface. - :param needs_adjusting: Surfaces whose curvatures need to be adjusted along the edges of the surface or constrained. - :param frequency_table: True if a frequency table is to be displayed. - :return: Return the filters, scalar ranges of curvatures and elevation along with the lookup tables. - """ - name = surface.name - source = surface.source - curvature = 'Mean_Curvature' - - curvatures = vtkCurvatures(curvature_type=Curvatures.CurvatureType().VTK_CURVATURE_MEAN) - p = (source >> curvatures).update().output - - if name in needs_adjusting: - adjust_edge_curvatures(p, curvature) - if name == 'plane': - constrain_curvatures(p, curvature, 0.0, 0.0) - if name == 'sphere': - # Mean curvature is 1/r - radius = 10 - mean_curvature = 1.0 / radius - constrain_curvatures(p, curvature, mean_curvature, mean_curvature) - - p.GetPointData().SetActiveScalars('Elevation') - elev_scalar_range = p.GetPointData().GetScalars('Elevation').range - # elev_lut = get_ordinal_lut() - # elev_labels = elev_lut.GetNumberOfTableValues(); - elev_lut = get_diverging_lut() - elev_labels = 5 - if name == 'plane': - elev_labels = 1 - elev_lut.SetTableRange(elev_scalar_range) - - p.GetPointData().SetActiveScalars(curvature) - curv_scalar_range = p.GetPointData().GetScalars(curvature).range - curv_lut = get_categorical_lut() - curv_lut.SetTableRange(curv_scalar_range) - curv_num_bands = curv_lut.GetNumberOfTableValues() - curv_bands = get_bands(curv_scalar_range, number_of_bands=curv_num_bands, precision=10, nearest_integer=False) - - # If any bands need adjusting, we would do it here. - - curv_freq = get_frequencies(curv_bands, p) - curv_bands, curv_freq = adjust_ranges(curv_bands, curv_freq) if frequency_table: - # Let's do a frequency table with the number of scalars in each band. - print_bands_frequencies(curvature, curv_bands, curv_freq) + print(f'{surface_name.title()} Elevation') + # The number of scalars in each band. + freq = get_frequencies(bands, source) + bands, freq = adjust_ranges(bands, freq) + print_bands_frequencies(bands, freq) - curv_lut.SetTableRange(curv_scalar_range) - curv_lut.SetNumberOfTableValues(len(curv_bands)) + # We will use the midpoint of the band as the label. + labels = [] + for k in bands: + labels.append('{:4.2f}'.format(bands[k][1])) - annotate_lut(curv_lut, curv_bands) + # Annotate + values = vtkVariantArray() + for i in range(len(labels)): + values.InsertNextValue(vtkVariant(labels[i])) + for i in range(values.GetNumberOfTuples()): + lut.SetAnnotation(i, values.GetValue(i).ToString()) # Create the contour bands. # We will use an indexed lookup table. - bcf = vtkBandedPolyDataContourFilter(input_data=p, + bcf = vtkBandedPolyDataContourFilter(input_data=source, scalar_mode=BandedPolyDataContourFilter.ScalarMode.VTK_SCALAR_MODE_INDEX, generate_contour_edges=True) - # Use either the minimum or maximum value for each band. - for k in curv_bands: - bcf.SetValue(k, curv_bands[k][2]) - - # Generate the glyphs on the original surface. - arrow_scale = (1, 1, 1) - if name == 'plane': - arrow_scale = (1.25, 1.25, 1.25) - scale_factor = 1.0 - if name == 'hills': - scale_factor = 0.5 - if name == 'sphere': - scale_factor = 2.0 - - glyph = get_glyphs(surface, arrow_scale=arrow_scale, scale_factor=scale_factor, reverse_normals=False) - - return {'bcf': bcf, 'glyph': glyph, 'curv_scalar_range': curv_scalar_range, - 'elev_scalar_range': elev_scalar_range, 'curv_lut': curv_lut, - 'elev_lut': elev_lut, 'elev_labels': elev_labels} + for i in range(len(bands)): + bcf.SetValue(i, bands[i][2]) + + glyphs = get_glyphs(source, scale_factor, reverse_normals=False) + + BandedGlyphs = namedtuple('BandedGlyphs', + ['glyph_type', 'glyphs', 'bcf', + 'lut_max_bands', 'lut', 'lutr', + 'lut1_max_bands', 'lut1', 'lut1r', + 'scalar_range', 'labels']) + bg = BandedGlyphs + bg.glyph_type = 'Elevation' + bg.glyphs = glyphs + bg.bcf = bcf + bg.lut_max_bands = max_lut_bands + bg.lut = lut + bg.lutr = reverse_lut(lut) + bg.lut1_max_bands = max_lut_bands + bg.lut1 = lut1 + bg.lut1r = reverse_lut(lut1) + bg.scalar_range = scalar_range + bg.labels = labels + + return bg class ScalarBarProperties: @@ -1221,11 +1208,72 @@ def make_scalar_bar_widget(scalar_bar_properties, title_text_property, label_tex sb_rep.position_coordinate.value = scalar_bar_properties.position_h['p'] sb_rep.position2_coordinate.value = scalar_bar_properties.position_h['p2'] - widget = vtkScalarBarWidget(representation=sb_rep, scalar_bar_actor=sb_actor, default_renderer=renderer, interactor=interactor, enabled=True) + widget = vtkScalarBarWidget(representation=sb_rep, scalar_bar_actor=sb_actor, default_renderer=renderer, + interactor=interactor, enabled=True) return widget +def position_sbw_h(num_bands, max_bands): + """ + Position the vertical scalar bar widget. + :param: num_bands - the number of bands in the scalar bar. + :param: max_bands - the maximum number of bands. + :return: The scalar bar position. + """ + + max_bands = abs(max_bands) + num_bands = abs(num_bands) + if num_bands > max_bands: + num_bands = max_bands + if num_bands == 0: + num_bands = 1 + # Origin of the scalar bar. + xy0 = [0.125, 0.05] + # Width and height of the scalar bar. + dxy = [0.75, 0.1] + if num_bands >= max_bands: + return {'p': tuple(xy0), 'p2': tuple(dxy)} + + dx = dxy[0] - xy0[0] * num_bands / max_bands + dxy[0] = dxy[0] * num_bands / max_bands + if num_bands == 1: + xy0[0] = 0.5 - dx * num_bands / (max_bands * 2) + else: + xy0[0] = 0.5 - dx * (num_bands + 1) / (max_bands * 2) + return {'p': tuple(xy0), 'p2': tuple(dxy)} + + +def position_sbw_v(num_bands, max_bands): + """ + Position the vertical scalar bar widget. + :param: num_bands - the number of bands in the scalar bar. + :param: max_bands - the maximum number of bands. + :return: The scalar bar position. + """ + + max_bands = abs(max_bands) + num_bands = abs(num_bands) + if num_bands > max_bands: + num_bands = max_bands + if num_bands == 0: + num_bands = 1 + # Origin of the scalar bar. + xy0 = [0.9, 0.25] + # Width and height of the scalar bar. + dxy = [0.08, 0.5] + if num_bands >= max_bands: + return {'p': tuple(xy0), 'p2': tuple(dxy)} + + dy = dxy[1] - xy0[1] * num_bands / max_bands + dxy[1] = dxy[1] * num_bands / max_bands + if num_bands == 1: + xy0[1] = 0.5 - dy * num_bands / (max_bands * 2) + else: + xy0[1] = 0.5 - dy * (num_bands + 1) / (max_bands * 2) + return {'p': tuple(xy0), 'p2': tuple(dxy)} + + def get_text_positions(names, justification=0, vertical_justification=0, width=0.96, height=0.1): """ Get viewport positioning information for a list of names. @@ -1288,6 +1336,83 @@ def get_text_positions(names, justification=0, vertical_justification=0, width=0 return text_positions +def print_bands_frequencies(bands, freq, precision=2): + """ + Print each band and the number of scalars in each band. + + :param bands: The bands. + :param freq: The frequencies. + :param precision: The precision for the ranges in each band. + """ + prec = abs(precision) + if prec > 14: + prec = 14 + + if len(bands) != len(freq): + print('Bands and Frequencies must be the same size.') + return + 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:{width}.{prec}f}]: {freq[k]:8d}\n' + else: + s += f'{q:{width}.{prec}f}, ' + width = 3 * width + 13 + s += f'{"Total":{width}s}{total:8d}\n' + print(s) + + +def adjust_camera_parameters(surface_name, ren): + """ + Adjust the camera parameters. + + :param surface_name: The name of the surface. + :param ren: The surface renderer. + + """ + if surface_name == 'hills': + camera = ren.active_camera + camera.position = (16.3424, 19.8311, 0.46492) + camera.focal_point = (0.209609, 0.432443, -1.18699) + camera.view_up = (-0.755535, 0.640179, -0.13906) + camera.distance = 25.2845 + camera.clipping_range = (13.1133, 37.6179) + elif surface_name == 'parametric hills': + camera = ren.GetActiveCamera() + camera.position = (10.9299, 59.1505, 24.9823) + camera.focal_point = (2.21692, 7.97545, 7.75135) + camera.view_up = (-0.230136, 0.345504, -0.909761) + camera.distance = 54.6966 + camera.clipping_range = (36.3006, 77.9852) + elif surface_name == 'parametric torus': + camera = ren.active_camera + camera.position = (-1.38419, 24.2883, 34.9246) + camera.focal_point = (-2.07248e-07, 3.63658e-06, 0.016056) + camera.view_up = (0.010284, 0.821007, -0.570825) + camera.distance = 42.5493 + camera.clipping_range = (25.2917, 64.5115) + elif surface_name == 'plane': + camera = ren.active_camera + camera.position = (-0.516003, 22.5763, 51.9171) + camera.focal_point = (-5.77108e-08, 0.500002, 5.80651e-06) + camera.view_up = (-0.000956134, 0.920254, -0.391321) + camera.distance = 56.4182 + camera.clipping_range = (36.7854, 81.268) + elif surface_name == 'torus': + camera = ren.active_camera + camera.position = (-2.02659, 35.5605, 51.1256) + camera.focal_point = (0, 0, 0.0160508) + camera.view_up = (0.010284, 0.821007, -0.570825) + camera.distance = 62.2964 + camera.clipping_range = (38.14, 92.8545) + + @dataclass(frozen=True) class BandedPolyDataContourFilter: @dataclass(frozen=True) diff --git a/src/PythonicAPI/Visualization/CurvatureBandsWithGlyphs.py b/src/PythonicAPI/Visualization/CurvatureBandsWithGlyphs.py index d1eed81cdd6..8952f5d6e6c 100644 --- a/src/PythonicAPI/Visualization/CurvatureBandsWithGlyphs.py +++ b/src/PythonicAPI/Visualization/CurvatureBandsWithGlyphs.py @@ -242,36 +242,26 @@ def main(argv): selectable=False, enabled=True) curv_sbp = ScalarBarProperties() - curv_sbp.title_text = curvature.replace('_', ' ') + '\n' + curv_sbp.title_text = curv_bg.glyph_type.replace('_', '\n') + '\n' curv_sbp.number_of_labels = len(curv_bg.labels) - # lut puts the lowest value at the top of the scalar bar. - # lutr puts the highest value at the top of the scalar bar. curv_sbp.lut = curv_bg.lut curv_sbp.orientation = False - max_bands = 9 - if surface_name == 'hills': - curv_sbp.position_h = position_sbw_h(9, max_bands) - elif surface_name == 'parametric hills': - curv_sbp.position_h = position_sbw_h(7, max_bands) - elif surface_name in ['plane', 'sphere']: - curv_sbp.position_h = position_sbw_h(1, max_bands) + if curv_sbp.number_of_labels == 1: + curv_sbp.position_h = position_sbw_h(curv_sbp.number_of_labels, curv_bg.lut_max_bands + 3) else: - curv_sbp.position_h = position_sbw_h(5, max_bands) + curv_sbp.position_h = position_sbw_h(curv_sbp.number_of_labels, curv_bg.lut_max_bands) curv_sb_widget = make_scalar_bar_widget(curv_sbp, title_text_property, label_text_property, ren, iren) elev_sbp = ScalarBarProperties() - elev_sbp.title_text = 'Elevation\n' + elev_sbp.title_text = elev_bg.glyph_type.replace('_', '\n') + '\n' elev_sbp.number_of_labels = len(elev_bg.labels) elev_sbp.lut = elev_bg.lutr elev_sbp.orientation = True - max_bands = 8 - if surface_name in ['hills', 'parametric hills']: - elev_sbp.position_v = position_sbw_v(8, max_bands) - elif surface_name == 'plane': - elev_sbp.position_v = position_sbw_v(1, max_bands) + if elev_sbp.number_of_labels == 1: + elev_sbp.position_v = position_sbw_v(elev_sbp.number_of_labels, elev_bg.lut_max_bands) else: - elev_sbp.position_v = position_sbw_v(5, max_bands) + elev_sbp.position_v = position_sbw_v(elev_sbp.number_of_labels, elev_bg.lut_max_bands) elev_sb_widget = make_scalar_bar_widget(elev_sbp, title_text_property, label_text_property, ren, iren) @@ -297,154 +287,6 @@ def main(argv): iren.Start() -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: The vtkPolyData object with the adjusted edge curvatures. - """ - - def point_neighbourhood(pt_id): - """ - Find the ids of the neighbors 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 = vtkGenerateIds(input_data=source, point_ids=True, cell_ids=False, - point_ids_array_name=array_name, cell_ids_array_name=array_name) - - edges = vtkFeatureEdges(boundary_edges=True, manifold_edges=False, - non_manifold_edges=False, feature_edges=False) - - (source >> id_filter >> edges).update() - - edge_array = edges.output.GetPointData().GetArray(array_name) - boundary_ids = [] - for i in range(edges.output.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: The vtkPolyData object with the constrained curvatures. - """ - - 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 generate_elevations(src): """ Generate elevations over the surface. @@ -822,6 +664,154 @@ def get_frequencies(bands, src): return freq +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: The vtkPolyData object with the adjusted edge curvatures. + """ + + def point_neighbourhood(pt_id): + """ + Find the ids of the neighbors 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 = vtkGenerateIds(input_data=source, point_ids=True, cell_ids=False, + point_ids_array_name=array_name, cell_ids_array_name=array_name) + + edges = vtkFeatureEdges(boundary_edges=True, manifold_edges=False, + non_manifold_edges=False, feature_edges=False) + + (source >> id_filter >> edges).update() + + edge_array = edges.output.GetPointData().GetArray(array_name) + boundary_ids = [] + for i in range(edges.output.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: The vtkPolyData object with the constrained curvatures. + """ + + 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 adjust_ranges(bands, freq): """ The bands and frequencies are adjusted so that the first and last @@ -885,21 +875,32 @@ def get_curvature_glyphs(surface_name, source, curvature, source.output.GetPointData().SetActiveScalars(curvature) scalar_range = source.output.point_data.GetScalars(curvature).range + max_lut_bands = 0 + color_series = vtkColorSeries() - if surface_name == 'parametric hills': - color_series.color_scheme = color_series.BREWER_DIVERGING_SPECTRAL_7 - elif surface_name == 'hills': + color_series = vtkColorSeries(color_series.BREWER_QUALITATIVE_SET3) + + if surface_name == 'hills': color_series.color_scheme = color_series.BREWER_DIVERGING_SPECTRAL_9 + max_lut_bands = 9 + elif surface_name == 'parametric hills': + color_series.color_scheme = color_series.BREWER_DIVERGING_SPECTRAL_7 + max_lut_bands = 7 elif surface_name in ['plane', 'sphere']: - color_series.color_scheme = color_series.CITRUS + color_series.color_scheme = color_series.BREWER_SEQUENTIAL_YELLOW_ORANGE_BROWN_3 + max_lut_bands = 3 else: - color_series.color_scheme = color_series.BREWER_DIVERGING_SPECTRAL_5 + color_series.color_scheme = color_series.BREWER_DIVERGING_SPECTRAL_6 + max_lut_bands = 6 lut = vtkLookupTable() color_series.BuildLookupTable(lut, color_series.CATEGORICAL) lut.SetNanColor(0, 0, 0, 1) lut.SetTableRange(scalar_range) + if surface_name in ['plane', 'sphere']: + lut.number_of_table_values = 1 + lut1 = vtkLookupTable() color_series.BuildLookupTable(lut1, color_series.ORDINAL) lut1.SetNanColor(0, 0, 0, 1) @@ -986,14 +987,17 @@ def get_curvature_glyphs(surface_name, source, curvature, BandedGlyphs = namedtuple('BandedGlyphs', ['glyph_type', 'glyphs', 'bcf', - 'lut', 'lutr', 'lut1', 'lut1r', + 'lut_max_bands', 'lut', 'lutr', + 'lut1_max_bands', 'lut1', 'lut1r', 'scalar_range', 'labels']) bg = BandedGlyphs bg.glyph_type = curvature bg.glyphs = glyphs bg.bcf = bcf + bg.lut_max_bands = max_lut_bands bg.lut = lut bg.lutr = reverse_lut(lut) + bg.lut1_max_bands = max_lut_bands bg.lut1 = lut1 bg.lut1r = reverse_lut(lut1) bg.scalar_range = scalar_range @@ -1024,17 +1028,24 @@ def get_elevation_glyphs(surface_name, source, source.point_data.active_scalars = 'Elevation' scalar_range = source.point_data.GetScalars('Elevation').range + max_lut_bands = 0 + color_series = vtkColorSeries() if surface_name in ['hills', 'parametric hills']: color_series.color_scheme = color_series.BREWER_DIVERGING_BROWN_BLUE_GREEN_8 + max_lut_bands = 8 else: - color_series.color_scheme = color_series.BREWER_DIVERGING_BROWN_BLUE_GREEN_5 + color_series.color_scheme = color_series.BREWER_DIVERGING_BROWN_BLUE_GREEN_6 + max_lut_bands = 6 lut = vtkLookupTable() color_series.BuildLookupTable(lut, color_series.CATEGORICAL) lut.SetNanColor(0, 0, 0, 1) lut.SetTableRange(scalar_range) + if surface_name in ['plane']: + lut.number_of_table_values = 1 + lut1 = vtkLookupTable() color_series.BuildLookupTable(lut1, color_series.ORDINAL) lut1.SetNanColor(0, 0, 0, 1) @@ -1100,14 +1111,17 @@ def get_elevation_glyphs(surface_name, source, BandedGlyphs = namedtuple('BandedGlyphs', ['glyph_type', 'glyphs', 'bcf', - 'lut', 'lutr', 'lut1', 'lut1r', + 'lut_max_bands', 'lut', 'lutr', + 'lut1_max_bands', 'lut1', 'lut1r', 'scalar_range', 'labels']) bg = BandedGlyphs bg.glyph_type = 'Elevation' bg.glyphs = glyphs bg.bcf = bcf + bg.lut_max_bands = max_lut_bands bg.lut = lut bg.lutr = reverse_lut(lut) + bg.lut1_max_bands = max_lut_bands bg.lut1 = lut1 bg.lut1r = reverse_lut(lut1) bg.scalar_range = scalar_range diff --git a/src/PythonicAPI/Visualization/ElevationBandsWithGlyphs.py b/src/PythonicAPI/Visualization/ElevationBandsWithGlyphs.py index dd822a8efab..2e1ed0141bf 100644 --- a/src/PythonicAPI/Visualization/ElevationBandsWithGlyphs.py +++ b/src/PythonicAPI/Visualization/ElevationBandsWithGlyphs.py @@ -195,20 +195,16 @@ def main(argv): selectable=False, enabled=True) elev_sbp = ScalarBarProperties() - elev_sbp.title_text = 'Elevation\n' + elev_sbp.title_text = elev_bg.glyph_type.replace('_', '\n') + '\n' elev_sbp.number_of_labels = len(elev_bg.labels) # lut puts the lowest value at the top of the vertical scalar bar. # lutr puts the highest value at the top of the vertical scalar bar. elev_sbp.lut = elev_bg.lutr elev_sbp.orientation = True - max_bands = 8 - if surface_name in ['hills', 'parametric hills']: - elev_sbp.position_v = position_sbw_v(8, max_bands) - elif surface_name == 'plane': - elev_sbp.position_v = position_sbw_v(1, max_bands) + if elev_sbp.number_of_labels == 1: + elev_sbp.position_v = position_sbw_v(elev_sbp.number_of_labels, elev_bg.lut_max_bands) else: - elev_sbp.position_v = position_sbw_v(5, max_bands) - + elev_sbp.position_v = position_sbw_v(elev_sbp.number_of_labels, elev_bg.lut_max_bands) scalar_bar_widget = make_scalar_bar_widget(elev_sbp, title_text_property, label_text_property, ren, iren) @@ -673,22 +669,29 @@ def get_elevation_glyphs(surface_name, source, if surface_name == 'hills': scale_factor = 0.5 elif surface_name == 'sphere': - scale_factor = 0.25 + scale_factor = 0.2 source.point_data.active_scalars = 'Elevation' scalar_range = source.point_data.GetScalars('Elevation').range + max_lut_bands = 0 + color_series = vtkColorSeries() if surface_name in ['hills', 'parametric hills']: color_series.color_scheme = color_series.BREWER_DIVERGING_BROWN_BLUE_GREEN_8 + max_lut_bands = 8 else: - color_series.color_scheme = color_series.BREWER_DIVERGING_BROWN_BLUE_GREEN_5 + color_series.color_scheme = color_series.BREWER_DIVERGING_BROWN_BLUE_GREEN_6 + max_lut_bands = 6 lut = vtkLookupTable() color_series.BuildLookupTable(lut, color_series.CATEGORICAL) lut.SetNanColor(0, 0, 0, 1) lut.SetTableRange(scalar_range) + if surface_name in ['plane']: + lut.number_of_table_values = 1 + lut1 = vtkLookupTable() color_series.BuildLookupTable(lut1, color_series.ORDINAL) lut1.SetNanColor(0, 0, 0, 1) @@ -754,14 +757,17 @@ def get_elevation_glyphs(surface_name, source, BandedGlyphs = namedtuple('BandedGlyphs', ['glyph_type', 'glyphs', 'bcf', - 'lut', 'lutr', 'lut1', 'lut1r', + 'lut_max_bands', 'lut', 'lutr', + 'lut1_max_bands', 'lut1', 'lut1r', 'scalar_range', 'labels']) bg = BandedGlyphs bg.glyph_type = 'Elevation' bg.glyphs = glyphs bg.bcf = bcf + bg.lut_max_bands = max_lut_bands bg.lut = lut bg.lutr = reverse_lut(lut) + bg.lut1_max_bands = max_lut_bands bg.lut1 = lut1 bg.lut1r = reverse_lut(lut1) bg.scalar_range = scalar_range -- GitLab From 560487a53f72f17315c699930a21cbb11fda7468 Mon Sep 17 00:00:00 2001 From: amaclean Date: Thu, 11 Jun 2026 16:02:11 +1000 Subject: [PATCH 3/5] Moving some code from PolyData to Visualization --- src/Cxx.md | 6 +++--- src/Cxx/PolyData/CMakeLists.txt | 4 ---- src/Cxx/Visualization/CMakeLists.txt | 4 ++++ src/Cxx/{PolyData => Visualization}/Curvatures.cxx | 0 .../{PolyData => Visualization}/CurvaturesAdjustEdges.cxx | 0 .../{PolyData => Visualization}/CurvaturesAdjustEdges.md | 0 src/Cxx/{PolyData => Visualization}/CurvaturesDemo.cxx | 0 src/Cxx/{PolyData => Visualization}/CurvaturesDemo.md | 0 src/Python.md | 6 +++--- src/Python/{PolyData => Visualization}/Curvatures.py | 0 .../{PolyData => Visualization}/CurvaturesAdjustEdges.md | 0 .../{PolyData => Visualization}/CurvaturesAdjustEdges.py | 0 src/Python/{PolyData => Visualization}/CurvaturesDemo.md | 0 src/Python/{PolyData => Visualization}/CurvaturesDemo.py | 0 src/PythonicAPI.md | 6 +++--- src/PythonicAPI/{PolyData => Visualization}/Curvatures.py | 0 .../{PolyData => Visualization}/CurvaturesAdjustEdges.md | 0 .../{PolyData => Visualization}/CurvaturesAdjustEdges.py | 0 .../CurvaturesNormalsElevations.md | 0 .../CurvaturesNormalsElevations.py | 0 .../Cxx/{PolyData => Visualization}/TestCurvatures.png | 0 .../TestCurvaturesAdjustEdges.png | 0 .../Cxx/{PolyData => Visualization}/TestCurvaturesDemo.png | 0 .../PolyData => Python/Visualization}/TestCurvatures.png | 0 .../TestCurvaturesAdjustEdges.png | 0 .../{PolyData => Visualization}/TestCurvaturesDemo.png | 0 .../Baseline/PythonicAPI/Visualization/TestCurvatures.png | 3 +++ .../TestCurvaturesAdjustEdges.png | 0 .../TestCurvaturesNormalsElevations.png | 0 29 files changed, 16 insertions(+), 13 deletions(-) rename src/Cxx/{PolyData => Visualization}/Curvatures.cxx (100%) rename src/Cxx/{PolyData => Visualization}/CurvaturesAdjustEdges.cxx (100%) rename src/Cxx/{PolyData => Visualization}/CurvaturesAdjustEdges.md (100%) rename src/Cxx/{PolyData => Visualization}/CurvaturesDemo.cxx (100%) rename src/Cxx/{PolyData => Visualization}/CurvaturesDemo.md (100%) rename src/Python/{PolyData => Visualization}/Curvatures.py (100%) mode change 100755 => 100644 rename src/Python/{PolyData => Visualization}/CurvaturesAdjustEdges.md (100%) rename src/Python/{PolyData => Visualization}/CurvaturesAdjustEdges.py (100%) mode change 100755 => 100644 rename src/Python/{PolyData => Visualization}/CurvaturesDemo.md (100%) rename src/Python/{PolyData => Visualization}/CurvaturesDemo.py (100%) mode change 100755 => 100644 rename src/PythonicAPI/{PolyData => Visualization}/Curvatures.py (100%) mode change 100755 => 100644 rename src/PythonicAPI/{PolyData => Visualization}/CurvaturesAdjustEdges.md (100%) rename src/PythonicAPI/{PolyData => Visualization}/CurvaturesAdjustEdges.py (100%) mode change 100755 => 100644 rename src/PythonicAPI/{PolyData => Visualization}/CurvaturesNormalsElevations.md (100%) rename src/PythonicAPI/{PolyData => Visualization}/CurvaturesNormalsElevations.py (100%) rename src/Testing/Baseline/Cxx/{PolyData => Visualization}/TestCurvatures.png (100%) rename src/Testing/Baseline/Cxx/{PolyData => Visualization}/TestCurvaturesAdjustEdges.png (100%) rename src/Testing/Baseline/Cxx/{PolyData => Visualization}/TestCurvaturesDemo.png (100%) rename src/Testing/Baseline/{PythonicAPI/PolyData => Python/Visualization}/TestCurvatures.png (100%) rename src/Testing/Baseline/Python/{PolyData => Visualization}/TestCurvaturesAdjustEdges.png (100%) rename src/Testing/Baseline/Python/{PolyData => Visualization}/TestCurvaturesDemo.png (100%) create mode 100644 src/Testing/Baseline/PythonicAPI/Visualization/TestCurvatures.png rename src/Testing/Baseline/PythonicAPI/{PolyData => Visualization}/TestCurvaturesAdjustEdges.png (100%) rename src/Testing/Baseline/PythonicAPI/{PolyData => Visualization}/TestCurvaturesNormalsElevations.png (100%) diff --git a/src/Cxx.md b/src/Cxx.md index a94b4fea5cf..fe7fd051f6b 100644 --- a/src/Cxx.md +++ b/src/Cxx.md @@ -319,9 +319,6 @@ These examples demonstrate how to create an display one of the many vtkParametri [ConvexHull](/Cxx/PolyData/ConvexHull) | Convex hull using vtkHull. [ConvexHullShrinkWrap](/Cxx/PolyData/ConvexHullShrinkWrap) | Convex hull using shrink wrapping. [CopyAllArrays](/Cxx/PolyData/CopyAllArrays) | Copy all arrays from one vtkPolyData to another. -[Curvatures](/Cxx/PolyData/Curvatures) | Compute Gaussian and Mean Curvatures. -[CurvaturesAdjustEdges](/Cxx/PolyData/CurvaturesAdjustEdges) | Get the Gaussian and Mean curvatures of a surface with adjustments for edge effects. -[CurvaturesDemo](/Cxx/PolyData/CurvaturesDemo) | Demonstrates how to get the Gaussian and Mean curvatures of a surface. [DataBounds](/Cxx/PolyData/DataBounds) | Get the minimum and maximum value in each dimension. (Axis aligned bounding box) [DataSetSurfaceFilter](/Cxx/PolyData/DataSetSurfaceFilter) | Convert vtkUnstructuredGrid to vtkPolyData. [DecimatePolyline](/Cxx/PolyData/DecimatePolyline) | Decimate polyline. @@ -1039,6 +1036,9 @@ See [this tutorial](http://www.vtk.org/Wiki/VTK/Tutorials/3DDataTypes) for a bri [Cursor3D](/Cxx/Visualization/Cursor3D) | [CursorShape](/Cxx/Visualization/CursorShape) | Change the shape of the cursor. [CurvatureBandsWithGlyphs](/Cxx/Visualization/CurvatureBandsWithGlyphs) | Demonstrates the coloring of a surface by partitioning the gaussian curvature of a surface into bands and using arrows to display the normals on the surface. +[Curvatures](/Cxx/Visualization/Curvatures) | Compute Gaussian and Mean Curvatures. +[CurvaturesAdjustEdges](/Cxx/Visualization/CurvaturesAdjustEdges) | Get the Gaussian and Mean curvatures of a surface with adjustments for edge effects. +[CurvaturesDemo](/Cxx/Visualization/CurvaturesDemo) | Demonstrates how to get the Gaussian and Mean curvatures of a surface. [CurvedReformation](/Cxx/Visualization/CurvedReformation) | Sample a volume with a curved surface. In medical imaging, this is often called curved multi planar reformation. [CutStructuredGrid](/Cxx/VisualizationAlgorithms/CutStructuredGrid) | Cut through structured grid with plane. The cut plane is shown solid shaded. A computational plane of constant k value is shown in wireframe for comparison. The colors correspond to flow density. Cutting surfaces are not necessarily planes: implicit functions such as spheres, cylinders, and quadrics can also be used. [Cutter](/Cxx/VisualizationAlgorithms/Cutter) | How to use vtkCutter by cutting through a cube. diff --git a/src/Cxx/PolyData/CMakeLists.txt b/src/Cxx/PolyData/CMakeLists.txt index 91a5376d9ea..cb114e93380 100644 --- a/src/Cxx/PolyData/CMakeLists.txt +++ b/src/Cxx/PolyData/CMakeLists.txt @@ -82,7 +82,6 @@ if (BUILD_TESTING) ColoredPoints CombineImportedActors ConvexHull - Curvatures DataSetSurfaceFilter DetermineArrayDataTypes DijkstraGraphGeodesicPath @@ -145,9 +144,6 @@ if (BUILD_TESTING) add_test(${KIT}-ConvexHull ${CMAKE_RUNTIME_OUTPUT_DIRECTORY}/${KIT}CxxTests TestConvexHull ${DATA}/cowHead.vtp) - add_test(${KIT}-Curvatures ${CMAKE_RUNTIME_OUTPUT_DIRECTORY}/${KIT}CxxTests - TestCurvatures ${DATA}/cowHead.vtp) - add_test(${KIT}-DataSetSurfaceFilter ${CMAKE_RUNTIME_OUTPUT_DIRECTORY}/${KIT}CxxTests TestDataSetSurfaceFilter ${DATA}/uGridEx.vtk) diff --git a/src/Cxx/Visualization/CMakeLists.txt b/src/Cxx/Visualization/CMakeLists.txt index bb8b4c24a6c..250d5df30db 100644 --- a/src/Cxx/Visualization/CMakeLists.txt +++ b/src/Cxx/Visualization/CMakeLists.txt @@ -111,6 +111,7 @@ if (BUILD_TESTING) CorrectlyRenderTranslucentGeometry CreateColorSeriesDemo CurvedReformation + Curvatures DepthSortPolyData EdgePoints ExtrudePolyDataAlongLine @@ -181,6 +182,9 @@ if (BUILD_TESTING) add_test(${KIT}-CreateColorSeriesDemo ${CMAKE_RUNTIME_OUTPUT_DIRECTORY}/${KIT}CxxTests TestCreateColorSeriesDemo Yellow) + add_test(${KIT}-Curvatures ${CMAKE_RUNTIME_OUTPUT_DIRECTORY}/${KIT}CxxTests + TestCurvatures ${DATA}/cowHead.vtp) + add_test(${KIT}-CurvedReformation ${CMAKE_RUNTIME_OUTPUT_DIRECTORY}/${KIT}CxxTests TestCurvedReformation ${DATA}/HeadMRVolume.mhd ${DATA}/polyline.vtk 200) diff --git a/src/Cxx/PolyData/Curvatures.cxx b/src/Cxx/Visualization/Curvatures.cxx similarity index 100% rename from src/Cxx/PolyData/Curvatures.cxx rename to src/Cxx/Visualization/Curvatures.cxx diff --git a/src/Cxx/PolyData/CurvaturesAdjustEdges.cxx b/src/Cxx/Visualization/CurvaturesAdjustEdges.cxx similarity index 100% rename from src/Cxx/PolyData/CurvaturesAdjustEdges.cxx rename to src/Cxx/Visualization/CurvaturesAdjustEdges.cxx diff --git a/src/Cxx/PolyData/CurvaturesAdjustEdges.md b/src/Cxx/Visualization/CurvaturesAdjustEdges.md similarity index 100% rename from src/Cxx/PolyData/CurvaturesAdjustEdges.md rename to src/Cxx/Visualization/CurvaturesAdjustEdges.md diff --git a/src/Cxx/PolyData/CurvaturesDemo.cxx b/src/Cxx/Visualization/CurvaturesDemo.cxx similarity index 100% rename from src/Cxx/PolyData/CurvaturesDemo.cxx rename to src/Cxx/Visualization/CurvaturesDemo.cxx diff --git a/src/Cxx/PolyData/CurvaturesDemo.md b/src/Cxx/Visualization/CurvaturesDemo.md similarity index 100% rename from src/Cxx/PolyData/CurvaturesDemo.md rename to src/Cxx/Visualization/CurvaturesDemo.md diff --git a/src/Python.md b/src/Python.md index 10174d7df75..8eb1ae74fe7 100644 --- a/src/Python.md +++ b/src/Python.md @@ -218,9 +218,6 @@ If you are new to VTK then these [tutorials](#tutorial) will help to get you sta [ColoredTriangle](/Python/PolyData/ColoredTriangle) | Creates a file TriangleColored.vtp. [CombinePolyData](/Python/Filtering/CombinePolyData) | This example combines Polydata objects, and displays the result to the screen. [ConnectivityFilter](/Python/Filtering/ConnectivityFilter) | Color any dataset type based on connectivity. -[Curvatures](/Python/PolyData/Curvatures) | Compute Gaussian, and Mean Curvatures. -[CurvaturesAdjustEdges](/Python/PolyData/CurvaturesAdjustEdges) | Get the Gaussian and Mean curvatures of a surface with adjustments for edge effects. -[CurvaturesDemo](/Python/PolyData/CurvaturesDemo) | Demonstrates how to get the Gaussian and Mean curvatures of a surface. [ExtractPolyLinesFromPolyData](/Python/PolyData/ExtractPolyLinesFromPolyData) | Extract polylines from polydata. [ExtractSelection](/Python/PolyData/ExtractSelection) |Extract selected points. [ExtractSelectionUsingCells](/Python/PolyData/ExtractSelectionUsingCells) | Extract cell, select cell. @@ -549,6 +546,9 @@ See [this tutorial](http://www.vtk.org/Wiki/VTK/Tutorials/3DDataTypes) for a bri [CreateColorSeriesDemo](/Python/Visualization/CreateColorSeriesDemo) | Create a custom vtkColorSeries. [CubeAxesActor](/Python/Visualization/CubeAxesActor) | Display three orthogonal axes with labels. [CurvatureBandsWithGlyphs](/Python/Visualization/CurvatureBandsWithGlyphs) | Demonstrates the coloring of a surface by partitioning the gaussian curvature of a surface into bands and using arrows to display the normals on the surface. +[Curvatures](/Python/Visualization/Curvatures) | Compute Gaussian, and Mean Curvatures. +[CurvaturesAdjustEdges](/Python/Visualization/CurvaturesAdjustEdges) | Get the Gaussian and Mean curvatures of a surface with adjustments for edge effects. +[CurvaturesDemo](/Python/Visualization/CurvaturesDemo) | Demonstrates how to get the Gaussian and Mean curvatures of a surface. [CutStructuredGrid](/Python/VisualizationAlgorithms/CutStructuredGrid) | Cut through structured grid with plane. The cut plane is shown solid shaded. A computational plane of constant k value is shown in wireframe for comparison. The colors correspond to flow density. Cutting surfaces are not necessarily planes: implicit functions such as spheres, cylinders, and quadrics can also be used. [Cutter](/Python/VisualizationAlgorithms/Cutter) | How to use vtkCutter by cutting through a cube. [DataSetSurface](/Python/VisualizationAlgorithms/DataSetSurface) | Cutting a hexahedron with a plane. The red line on the surface shows the cut. diff --git a/src/Python/PolyData/Curvatures.py b/src/Python/Visualization/Curvatures.py old mode 100755 new mode 100644 similarity index 100% rename from src/Python/PolyData/Curvatures.py rename to src/Python/Visualization/Curvatures.py diff --git a/src/Python/PolyData/CurvaturesAdjustEdges.md b/src/Python/Visualization/CurvaturesAdjustEdges.md similarity index 100% rename from src/Python/PolyData/CurvaturesAdjustEdges.md rename to src/Python/Visualization/CurvaturesAdjustEdges.md diff --git a/src/Python/PolyData/CurvaturesAdjustEdges.py b/src/Python/Visualization/CurvaturesAdjustEdges.py old mode 100755 new mode 100644 similarity index 100% rename from src/Python/PolyData/CurvaturesAdjustEdges.py rename to src/Python/Visualization/CurvaturesAdjustEdges.py diff --git a/src/Python/PolyData/CurvaturesDemo.md b/src/Python/Visualization/CurvaturesDemo.md similarity index 100% rename from src/Python/PolyData/CurvaturesDemo.md rename to src/Python/Visualization/CurvaturesDemo.md diff --git a/src/Python/PolyData/CurvaturesDemo.py b/src/Python/Visualization/CurvaturesDemo.py old mode 100755 new mode 100644 similarity index 100% rename from src/Python/PolyData/CurvaturesDemo.py rename to src/Python/Visualization/CurvaturesDemo.py diff --git a/src/PythonicAPI.md b/src/PythonicAPI.md index 5081f1fb643..59856f07380 100644 --- a/src/PythonicAPI.md +++ b/src/PythonicAPI.md @@ -206,9 +206,6 @@ This Python script, [SelectExamples](../PythonicAPI/Utilities/SelectExamples), w [CellsInsideObject](/PythonicAPI/PolyData/CellsInsideObject) | Extract cells inside a closed surface. [CenterOfMass](/PythonicAPI/PolyData/CenterOfMass) | Compute the center of mass of the points. [ConnectivityFilter](/PythonicAPI/Filtering/ConnectivityFilter) | Color any dataset type based on connectivity. -[Curvatures](/PythonicAPI/PolyData/Curvatures) | Compute either Mean or Gaussian Curvatures with adjustments for edge effects. -[CurvaturesAdjustEdges](/PythonicAPI/PolyData/CurvaturesAdjustEdges) | Get the Gaussian and Mean curvatures of a surface with adjustments for edge effects. -[CurvaturesNormalsElevations](/PythonicAPI/PolyData/CurvaturesNormalsElevations) | Gaussian and Mean curvatures of a surface with arrows colored by elevation to display the normals. The surfaces are also adjusted for edge effects. [DecimatePolyline](/PythonicAPI/PolyData/DecimatePolyline) | Decimate polyline. [DistancePolyDataFilter](/PythonicAPI/PolyData/DistancePolyDataFilter) | Compute the distance function from one vtkPolyData to another. [ExternalContour](/PythonicAPI/PolyData/ExternalContour) | Get the external contour of a PolyData object. @@ -581,6 +578,9 @@ See [this tutorial](http://www.vtk.org/Wiki/VTK/Tutorials/3DDataTypes) for a bri [Cursor2D](/PythonicAPI/Visualization/Cursor2D) | [Cursor3D](/PythonicAPI/Visualization/Cursor3D) | [CurvatureBandsWithGlyphs](/PythonicAPI/Visualization/CurvatureBandsWithGlyphs) | Demonstrates the coloring of a surface by partitioning the gaussian curvature of a surface into bands and using arrows to display the normals on the surface. +[Curvatures](/PythonicAPI/Visualization/Curvatures) | Compute either Mean or Gaussian Curvatures with adjustments for edge effects. +[CurvaturesAdjustEdges](/PythonicAPI/Visualization/CurvaturesAdjustEdges) | Get the Gaussian and Mean curvatures of a surface with adjustments for edge effects. +[CurvaturesNormalsElevations](/PythonicAPI/Visualization/CurvaturesNormalsElevations) | Gaussian and Mean curvatures of a surface with arrows colored by elevation to display the normals. The surfaces are also adjusted for edge effects. [DataSetSurface](/PythonicAPI/VisualizationAlgorithms/DataSetSurface) | Cutting a hexahedron with a plane. The red line on the surface shows the cut. [DepthSortPolyData](/PythonicAPI/Visualization/DepthSortPolyData) | Poly Data Depth Sorting. [DisplacementPlot](/PythonicAPI/VisualizationAlgorithms/DisplacementPlot) | Show modal lines for a vibrating beam. diff --git a/src/PythonicAPI/PolyData/Curvatures.py b/src/PythonicAPI/Visualization/Curvatures.py old mode 100755 new mode 100644 similarity index 100% rename from src/PythonicAPI/PolyData/Curvatures.py rename to src/PythonicAPI/Visualization/Curvatures.py diff --git a/src/PythonicAPI/PolyData/CurvaturesAdjustEdges.md b/src/PythonicAPI/Visualization/CurvaturesAdjustEdges.md similarity index 100% rename from src/PythonicAPI/PolyData/CurvaturesAdjustEdges.md rename to src/PythonicAPI/Visualization/CurvaturesAdjustEdges.md diff --git a/src/PythonicAPI/PolyData/CurvaturesAdjustEdges.py b/src/PythonicAPI/Visualization/CurvaturesAdjustEdges.py old mode 100755 new mode 100644 similarity index 100% rename from src/PythonicAPI/PolyData/CurvaturesAdjustEdges.py rename to src/PythonicAPI/Visualization/CurvaturesAdjustEdges.py diff --git a/src/PythonicAPI/PolyData/CurvaturesNormalsElevations.md b/src/PythonicAPI/Visualization/CurvaturesNormalsElevations.md similarity index 100% rename from src/PythonicAPI/PolyData/CurvaturesNormalsElevations.md rename to src/PythonicAPI/Visualization/CurvaturesNormalsElevations.md diff --git a/src/PythonicAPI/PolyData/CurvaturesNormalsElevations.py b/src/PythonicAPI/Visualization/CurvaturesNormalsElevations.py similarity index 100% rename from src/PythonicAPI/PolyData/CurvaturesNormalsElevations.py rename to src/PythonicAPI/Visualization/CurvaturesNormalsElevations.py diff --git a/src/Testing/Baseline/Cxx/PolyData/TestCurvatures.png b/src/Testing/Baseline/Cxx/Visualization/TestCurvatures.png similarity index 100% rename from src/Testing/Baseline/Cxx/PolyData/TestCurvatures.png rename to src/Testing/Baseline/Cxx/Visualization/TestCurvatures.png diff --git a/src/Testing/Baseline/Cxx/PolyData/TestCurvaturesAdjustEdges.png b/src/Testing/Baseline/Cxx/Visualization/TestCurvaturesAdjustEdges.png similarity index 100% rename from src/Testing/Baseline/Cxx/PolyData/TestCurvaturesAdjustEdges.png rename to src/Testing/Baseline/Cxx/Visualization/TestCurvaturesAdjustEdges.png diff --git a/src/Testing/Baseline/Cxx/PolyData/TestCurvaturesDemo.png b/src/Testing/Baseline/Cxx/Visualization/TestCurvaturesDemo.png similarity index 100% rename from src/Testing/Baseline/Cxx/PolyData/TestCurvaturesDemo.png rename to src/Testing/Baseline/Cxx/Visualization/TestCurvaturesDemo.png diff --git a/src/Testing/Baseline/PythonicAPI/PolyData/TestCurvatures.png b/src/Testing/Baseline/Python/Visualization/TestCurvatures.png similarity index 100% rename from src/Testing/Baseline/PythonicAPI/PolyData/TestCurvatures.png rename to src/Testing/Baseline/Python/Visualization/TestCurvatures.png diff --git a/src/Testing/Baseline/Python/PolyData/TestCurvaturesAdjustEdges.png b/src/Testing/Baseline/Python/Visualization/TestCurvaturesAdjustEdges.png similarity index 100% rename from src/Testing/Baseline/Python/PolyData/TestCurvaturesAdjustEdges.png rename to src/Testing/Baseline/Python/Visualization/TestCurvaturesAdjustEdges.png diff --git a/src/Testing/Baseline/Python/PolyData/TestCurvaturesDemo.png b/src/Testing/Baseline/Python/Visualization/TestCurvaturesDemo.png similarity index 100% rename from src/Testing/Baseline/Python/PolyData/TestCurvaturesDemo.png rename to src/Testing/Baseline/Python/Visualization/TestCurvaturesDemo.png diff --git a/src/Testing/Baseline/PythonicAPI/Visualization/TestCurvatures.png b/src/Testing/Baseline/PythonicAPI/Visualization/TestCurvatures.png new file mode 100644 index 00000000000..5a7af1583d4 --- /dev/null +++ b/src/Testing/Baseline/PythonicAPI/Visualization/TestCurvatures.png @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:ede72d272f11b059d09af2dd5f15e4ed2ed6d9ecce27d2c7f716a21454502d61 +size 129918 diff --git a/src/Testing/Baseline/PythonicAPI/PolyData/TestCurvaturesAdjustEdges.png b/src/Testing/Baseline/PythonicAPI/Visualization/TestCurvaturesAdjustEdges.png similarity index 100% rename from src/Testing/Baseline/PythonicAPI/PolyData/TestCurvaturesAdjustEdges.png rename to src/Testing/Baseline/PythonicAPI/Visualization/TestCurvaturesAdjustEdges.png diff --git a/src/Testing/Baseline/PythonicAPI/PolyData/TestCurvaturesNormalsElevations.png b/src/Testing/Baseline/PythonicAPI/Visualization/TestCurvaturesNormalsElevations.png similarity index 100% rename from src/Testing/Baseline/PythonicAPI/PolyData/TestCurvaturesNormalsElevations.png rename to src/Testing/Baseline/PythonicAPI/Visualization/TestCurvaturesNormalsElevations.png -- GitLab From e4c9efa313cb8930bb6b69ae7d4eee1931f1d86d Mon Sep 17 00:00:00 2001 From: amaclean Date: Fri, 12 Jun 2026 16:57:59 +1000 Subject: [PATCH 4/5] Making the code more PythonicAPI --- .../Visualization/CurvatureBandsWithGlyphs.py | 57 +++++++-------- src/PythonicAPI/Visualization/Curvatures.py | 43 ++++++----- .../Visualization/CurvaturesAdjustEdges.py | 42 ++++++----- .../CurvaturesNormalsElevations.py | 71 +++++++++++-------- .../Visualization/ElevationBandsWithGlyphs.py | 16 +++-- 5 files changed, 127 insertions(+), 102 deletions(-) diff --git a/src/PythonicAPI/Visualization/CurvatureBandsWithGlyphs.py b/src/PythonicAPI/Visualization/CurvatureBandsWithGlyphs.py index 8952f5d6e6c..316ef18eabd 100644 --- a/src/PythonicAPI/Visualization/CurvatureBandsWithGlyphs.py +++ b/src/PythonicAPI/Visualization/CurvatureBandsWithGlyphs.py @@ -4,6 +4,7 @@ import copy import math from collections import namedtuple from dataclasses import dataclass +from pathlib import Path import numpy as np # noinspection PyUnresolvedReferences @@ -202,7 +203,7 @@ def main(argv): # ------------------------------------------------------------ ren = vtkRenderer(background=colors.GetColor3d('ParaViewBlueGrayBkg')) ren_win = vtkRenderWindow(size=(window_width, window_height), - window_name='CurvatureBandsWithGlyphs') + window_name=f'{Path(argv[0]).stem:s}') ren_win.AddRenderer(ren) iren = vtkRenderWindowInteractor() iren.render_window = ren_win @@ -692,11 +693,11 @@ def adjust_edge_curvatures(source, curvature_name, epsilon=1.0e-08): cell_ids = vtkIdList() source.GetPointCells(pt_id, cell_ids) neighbour = set() - for cell_idx in range(0, cell_ids.GetNumberOfIds()): + for cell_idx in range(0, cell_ids.number_of_ids): 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()): + for cell_pt_idx in range(0, cell_point_ids.number_of_ids): neighbour.add(cell_point_ids.GetId(cell_pt_idx)) return neighbour @@ -713,7 +714,7 @@ def adjust_edge_curvatures(source, curvature_name, epsilon=1.0e-08): return np.linalg.norm(pt_a - pt_b) # Get the active scalars - source.GetPointData().SetActiveScalars(curvature_name) + source.point_data.active_scalars = curvature_name np_source = dsa.WrapDataObject(source) curvatures = np_source.PointData[curvature_name] @@ -727,9 +728,9 @@ def adjust_edge_curvatures(source, curvature_name, epsilon=1.0e-08): (source >> id_filter >> edges).update() - edge_array = edges.output.GetPointData().GetArray(array_name) + edge_array = edges.output.point_data.GetArray(array_name) boundary_ids = [] - for i in range(edges.output.GetNumberOfPoints()): + for i in range(edges.output.number_of_points): boundary_ids.append(edge_array.GetValue(i)) # Remove duplicate Ids. p_ids_set = set(boundary_ids) @@ -767,10 +768,10 @@ def adjust_edge_curvatures(source, curvature_name, epsilon=1.0e-08): 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) + curv.name = curvature_name + source.point_data.RemoveArray(curvature_name) + source.point_data.AddArray(curv) + source.point_data.active_scalars = curvature_name def constrain_curvatures(source, curvature_name, lower_bound=0.0, upper_bound=0.0): @@ -795,7 +796,7 @@ def constrain_curvatures(source, curvature_name, lower_bound=0.0, upper_bound=0. bounds.append(lower_bound) # Get the active scalars. - source.GetPointData().SetActiveScalars(curvature_name) + source.point_data.active_scalars = curvature_name np_source = dsa.WrapDataObject(source) curvatures = np_source.PointData[curvature_name] @@ -806,10 +807,10 @@ def constrain_curvatures(source, curvature_name, lower_bound=0.0, upper_bound=0. 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) + curv.name = curvature_name + source.point_data.RemoveArray(curvature_name) + source.point_data.AddArray(curv) + source.point_data.active_scalars = curvature_name def adjust_ranges(bands, freq): @@ -872,7 +873,7 @@ def get_curvature_glyphs(surface_name, source, curvature, elif surface_name == 'sphere': scale_factor = 0.2 - source.output.GetPointData().SetActiveScalars(curvature) + source.output.point_data.SetActiveScalars(curvature) scalar_range = source.output.point_data.GetScalars(curvature).range max_lut_bands = 0 @@ -895,16 +896,16 @@ def get_curvature_glyphs(surface_name, source, curvature, lut = vtkLookupTable() color_series.BuildLookupTable(lut, color_series.CATEGORICAL) - lut.SetNanColor(0, 0, 0, 1) - lut.SetTableRange(scalar_range) + lut.nan_color = (0, 0, 0, 1) + lut.table_range = scalar_range if surface_name in ['plane', 'sphere']: lut.number_of_table_values = 1 lut1 = vtkLookupTable() color_series.BuildLookupTable(lut1, color_series.ORDINAL) - lut1.SetNanColor(0, 0, 0, 1) - lut1.SetTableRange(scalar_range) + lut1.nan_color = (0, 0, 0, 1) + lut1.table_range = scalar_range number_of_bands = lut.number_of_table_values lut1.number_of_table_values = number_of_bands @@ -944,9 +945,9 @@ def get_curvature_glyphs(surface_name, source, curvature, # Adjust the number of table values and scalar range. scalar_range = (bands[0][0], bands[len(bands) - 1][2]) - lut.TableRange = scalar_range + lut.table_range = scalar_range lut.number_of_table_values = len(bands) - lut1.TableRange = scalar_range + lut1.table_range = scalar_range lut1.number_of_table_values = len(bands) if frequency_table: @@ -1040,16 +1041,16 @@ def get_elevation_glyphs(surface_name, source, lut = vtkLookupTable() color_series.BuildLookupTable(lut, color_series.CATEGORICAL) - lut.SetNanColor(0, 0, 0, 1) - lut.SetTableRange(scalar_range) + lut.nan_color = (0, 0, 0, 1) + lut.table_range = scalar_range if surface_name in ['plane']: lut.number_of_table_values = 1 lut1 = vtkLookupTable() color_series.BuildLookupTable(lut1, color_series.ORDINAL) - lut1.SetNanColor(0, 0, 0, 1) - lut1.SetTableRange(scalar_range) + lut1.nan_color = (0, 0, 0, 1) + lut1.table_range = scalar_range number_of_bands = lut.number_of_table_values lut1.number_of_table_values = number_of_bands @@ -1074,9 +1075,9 @@ def get_elevation_glyphs(surface_name, source, # Adjust the number of table values and scalar range. scalar_range = (bands[0][0], bands[len(bands) - 1][2]) - lut.TableRange = scalar_range + lut.table_range = scalar_range lut.number_of_table_values = len(bands) - lut1.TableRange = scalar_range + lut1.table_range = scalar_range lut1.number_of_table_values = len(bands) if frequency_table: diff --git a/src/PythonicAPI/Visualization/Curvatures.py b/src/PythonicAPI/Visualization/Curvatures.py index aeaf1788771..d7d74df03da 100644 --- a/src/PythonicAPI/Visualization/Curvatures.py +++ b/src/PythonicAPI/Visualization/Curvatures.py @@ -116,14 +116,13 @@ def main(argv): window_height = 800 # Create a renderer, render window, and interactor - renderer = vtkRenderer(background=colors.GetColor3d('DarkSlateGray')) - ren_win = vtkRenderWindow(size=(window_width, window_height), window_name='Curvatures') + renderer = vtkRenderer(background=colors.GetColor3d('ParaViewBlueGrayBkg')) + ren_win = vtkRenderWindow(size=(window_width, window_height), + window_name=f'{Path(argv[0]).stem:s}') ren_win.AddRenderer(renderer) iren = vtkRenderWindowInteractor() iren.render_window = ren_win - # Important: The interactor must be set prior to enabling the widget. - iren.render_window = ren_win text_property = vtkTextProperty(color=colors.GetColor3d('AliceBlue'), bold=True, italic=True, shadow=True, font_size=16, @@ -145,7 +144,6 @@ def main(argv): # Add the actors to the scene renderer.AddActor(actor) - renderer.SetBackground(colors.GetColor3d('DarkSlateGray')) # Render and interact ren_win.Render() @@ -157,27 +155,34 @@ 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. - :param source: The vtkCurvatures object. + 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: + :return: The vtkPolyData object with the adjusted edge curvatures. """ def point_neighbourhood(pt_id): """ - Extract the topological neighbors for point. + Find the ids of the neighbors 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()): + for cell_idx in range(0, cell_ids.number_of_ids): 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()): + for cell_pt_idx in range(0, cell_point_ids.number_of_ids): neighbour.add(cell_point_ids.GetId(cell_pt_idx)) return neighbour @@ -185,33 +190,32 @@ def adjust_edge_curvatures(source, curvature_name, epsilon=1.0e-08): """ Compute the distance between two points given their ids. - :param pt_id_a: First point. - :param pt_id_b: Second point. - :return: The distance. + :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.point_data.SetActiveScalars(curvature_name) + source.point_data.active_scalars = curvature_name np_source = dsa.WrapDataObject(source) curvatures = np_source.PointData[curvature_name] # Get the boundary point IDs. array_name = 'ids' - id_filter = vtkGenerateIds(point_ids=True, cell_ids=False, - point_ids_array_name=array_name, - cell_ids_array_name=array_name) + id_filter = vtkGenerateIds(input_data=source, point_ids=True, cell_ids=False, + point_ids_array_name=array_name, cell_ids_array_name=array_name) edges = vtkFeatureEdges(boundary_edges=True, manifold_edges=False, non_manifold_edges=False, feature_edges=False) (source >> id_filter >> edges).update() - edge_array = edges.output.GetPointData().GetArray(array_name) + edge_array = edges.output.point_data.GetArray(array_name) boundary_ids = [] - for i in range(edges.output.GetNumberOfPoints()): + for i in range(edges.output.number_of_points): boundary_ids.append(edge_array.GetValue(i)) # Remove duplicate Ids. p_ids_set = set(boundary_ids) @@ -245,6 +249,7 @@ def adjust_edge_curvatures(source, curvature_name, epsilon=1.0e-08): # 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) diff --git a/src/PythonicAPI/Visualization/CurvaturesAdjustEdges.py b/src/PythonicAPI/Visualization/CurvaturesAdjustEdges.py index 81558568e04..3b2c4af3fca 100644 --- a/src/PythonicAPI/Visualization/CurvaturesAdjustEdges.py +++ b/src/PythonicAPI/Visualization/CurvaturesAdjustEdges.py @@ -367,27 +367,34 @@ 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. - :param source: The vtkCurvatures object. + 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: + :return: The vtkPolyData object with the adjusted edge curvatures. """ def point_neighbourhood(pt_id): """ - Extract the topological neighbors for point. + Find the ids of the neighbors 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()): + for cell_idx in range(0, cell_ids.number_of_ids): 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()): + for cell_pt_idx in range(0, cell_point_ids.number_of_ids): neighbour.add(cell_point_ids.GetId(cell_pt_idx)) return neighbour @@ -395,33 +402,32 @@ def adjust_edge_curvatures(source, curvature_name, epsilon=1.0e-08): """ Compute the distance between two points given their ids. - :param pt_id_a: First point. - :param pt_id_b: Second point. - :return: The distance. + :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.point_data.SetActiveScalars(curvature_name) + source.point_data.active_scalars = curvature_name np_source = dsa.WrapDataObject(source) curvatures = np_source.PointData[curvature_name] # Get the boundary point IDs. array_name = 'ids' - id_filter = vtkGenerateIds(point_ids=True, cell_ids=False, - point_ids_array_name=array_name, - cell_ids_array_name=array_name) + id_filter = vtkGenerateIds(input_data=source, point_ids=True, cell_ids=False, + point_ids_array_name=array_name, cell_ids_array_name=array_name) edges = vtkFeatureEdges(boundary_edges=True, manifold_edges=False, non_manifold_edges=False, feature_edges=False) (source >> id_filter >> edges).update() - edge_array = edges.output.GetPointData().GetArray(array_name) + edge_array = edges.output.point_data.GetArray(array_name) boundary_ids = [] - for i in range(edges.output.GetNumberOfPoints()): + for i in range(edges.output.number_of_points): boundary_ids.append(edge_array.GetValue(i)) # Remove duplicate Ids. p_ids_set = set(boundary_ids) @@ -455,6 +461,7 @@ def adjust_edge_curvatures(source, curvature_name, epsilon=1.0e-08): # 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) @@ -474,7 +481,7 @@ def constrain_curvatures(source, curvature_name, lower_bound=0.0, upper_bound=0. :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: + :return: The vtkPolyData object with the constrained curvatures. """ bounds = list() @@ -485,14 +492,15 @@ def constrain_curvatures(source, curvature_name, lower_bound=0.0, upper_bound=0. bounds.append(upper_bound) bounds.append(lower_bound) - # Get the active scalars - source.point_data.SetActiveScalars(curvature_name) + # Get the active scalars. + source.point_data.active_scalars = 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) diff --git a/src/PythonicAPI/Visualization/CurvaturesNormalsElevations.py b/src/PythonicAPI/Visualization/CurvaturesNormalsElevations.py index 2a5f64a3d68..82db6e20a86 100644 --- a/src/PythonicAPI/Visualization/CurvaturesNormalsElevations.py +++ b/src/PythonicAPI/Visualization/CurvaturesNormalsElevations.py @@ -182,7 +182,8 @@ def main(argv): # -------------------------------------------------- # Create the RenderWindow, Renderers and Interactor. # -------------------------------------------------- - ren_win = vtkRenderWindow(size=(window_width, window_height), window_name=f'{Path(argv[0]).name:s}') + ren_win = vtkRenderWindow(size=(window_width, window_height), + window_name=f'{Path(argv[0]).stem:s}') iren = vtkRenderWindowInteractor() iren.render_window = ren_win style = vtkInteractorStyleTrackballCamera() @@ -700,27 +701,34 @@ 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. - :param source: The vtkCurvatures object. + 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: + :return: The vtkPolyData object with the adjusted edge curvatures. """ def point_neighbourhood(pt_id): """ - Extract the topological neighbors for point. + Find the ids of the neighbors 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()): + for cell_idx in range(0, cell_ids.number_of_ids): 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()): + for cell_pt_idx in range(0, cell_point_ids.number_of_ids): neighbour.add(cell_point_ids.GetId(cell_pt_idx)) return neighbour @@ -728,33 +736,32 @@ def adjust_edge_curvatures(source, curvature_name, epsilon=1.0e-08): """ Compute the distance between two points given their ids. - :param pt_id_a: First point. - :param pt_id_b: Second point. - :return: The distance. + :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.point_data.SetActiveScalars(curvature_name) + source.point_data.active_scalars = curvature_name np_source = dsa.WrapDataObject(source) curvatures = np_source.PointData[curvature_name] # Get the boundary point IDs. array_name = 'ids' - id_filter = vtkGenerateIds(point_ids=True, cell_ids=False, - point_ids_array_name=array_name, - cell_ids_array_name=array_name) + id_filter = vtkGenerateIds(input_data=source, point_ids=True, cell_ids=False, + point_ids_array_name=array_name, cell_ids_array_name=array_name) edges = vtkFeatureEdges(boundary_edges=True, manifold_edges=False, non_manifold_edges=False, feature_edges=False) (source >> id_filter >> edges).update() - edge_array = edges.output.GetPointData().GetArray(array_name) + edge_array = edges.output.point_data.GetArray(array_name) boundary_ids = [] - for i in range(edges.output.GetNumberOfPoints()): + for i in range(edges.output.number_of_points): boundary_ids.append(edge_array.GetValue(i)) # Remove duplicate Ids. p_ids_set = set(boundary_ids) @@ -788,6 +795,7 @@ def adjust_edge_curvatures(source, curvature_name, epsilon=1.0e-08): # 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) @@ -807,7 +815,7 @@ def constrain_curvatures(source, curvature_name, lower_bound=0.0, upper_bound=0. :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: + :return: The vtkPolyData object with the constrained curvatures. """ bounds = list() @@ -818,14 +826,15 @@ def constrain_curvatures(source, curvature_name, lower_bound=0.0, upper_bound=0. bounds.append(upper_bound) bounds.append(lower_bound) - # Get the active scalars - source.point_data.SetActiveScalars(curvature_name) + # Get the active scalars. + source.point_data.active_scalars = 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) @@ -895,7 +904,7 @@ def get_curvature_glyphs(surface_name, source, curvature, elif surface_name == 'sphere': scale_factor = 0.2 - source.output.GetPointData().SetActiveScalars(curvature) + source.output.point_data.SetActiveScalars(curvature) scalar_range = source.output.point_data.GetScalars(curvature).range max_lut_bands = 0 @@ -918,16 +927,16 @@ def get_curvature_glyphs(surface_name, source, curvature, lut = vtkLookupTable() color_series.BuildLookupTable(lut, color_series.CATEGORICAL) - lut.SetNanColor(0, 0, 0, 1) - lut.SetTableRange(scalar_range) + lut.nan_color = (0, 0, 0, 1) + lut.table_range = scalar_range if surface_name in ['plane', 'sphere']: lut.number_of_table_values = 1 lut1 = vtkLookupTable() color_series.BuildLookupTable(lut1, color_series.ORDINAL) - lut1.SetNanColor(0, 0, 0, 1) - lut1.SetTableRange(scalar_range) + lut1.nan_color = (0, 0, 0, 1) + lut1.table_range = scalar_range number_of_bands = lut.number_of_table_values lut1.number_of_table_values = number_of_bands @@ -967,9 +976,9 @@ def get_curvature_glyphs(surface_name, source, curvature, # Adjust the number of table values and scalar range. scalar_range = (bands[0][0], bands[len(bands) - 1][2]) - lut.TableRange = scalar_range + lut.table_range = scalar_range lut.number_of_table_values = len(bands) - lut1.TableRange = scalar_range + lut1.table_range = scalar_range lut1.number_of_table_values = len(bands) if frequency_table: @@ -1063,16 +1072,16 @@ def get_elevation_glyphs(surface_name, source, lut = vtkLookupTable() color_series.BuildLookupTable(lut, color_series.CATEGORICAL) - lut.SetNanColor(0, 0, 0, 1) - lut.SetTableRange(scalar_range) + lut.nan_color = (0, 0, 0, 1) + lut.table_range = scalar_range if surface_name in ['plane']: lut.number_of_table_values = 1 lut1 = vtkLookupTable() color_series.BuildLookupTable(lut1, color_series.ORDINAL) - lut1.SetNanColor(0, 0, 0, 1) - lut1.SetTableRange(scalar_range) + lut1.nan_color = (0, 0, 0, 1) + lut1.table_range = scalar_range number_of_bands = lut.number_of_table_values lut1.number_of_table_values = number_of_bands @@ -1097,9 +1106,9 @@ def get_elevation_glyphs(surface_name, source, # Adjust the number of table values and scalar range. scalar_range = (bands[0][0], bands[len(bands) - 1][2]) - lut.TableRange = scalar_range + lut.table_range = scalar_range lut.number_of_table_values = len(bands) - lut1.TableRange = scalar_range + lut1.table_range = scalar_range lut1.number_of_table_values = len(bands) if frequency_table: diff --git a/src/PythonicAPI/Visualization/ElevationBandsWithGlyphs.py b/src/PythonicAPI/Visualization/ElevationBandsWithGlyphs.py index 2e1ed0141bf..9ef2e3d9da1 100644 --- a/src/PythonicAPI/Visualization/ElevationBandsWithGlyphs.py +++ b/src/PythonicAPI/Visualization/ElevationBandsWithGlyphs.py @@ -4,6 +4,7 @@ import copy import math from collections import namedtuple from dataclasses import dataclass +from pathlib import Path # noinspection PyUnresolvedReferences import vtkmodules.vtkRenderingOpenGL2 @@ -164,7 +165,8 @@ def main(argv): # ------------------------------------------------------------ ren = vtkRenderer(background=colors.GetColor3d('ParaViewBlueGrayBkg')) ren_win = vtkRenderWindow(size=(window_width, window_height), - window_name='ElevationBandsWithGlyphs') + window_name=f'{Path(argv[0]).stem:s}') + ren_win.AddRenderer(ren) iren = vtkRenderWindowInteractor() iren.render_window = ren_win @@ -686,16 +688,16 @@ def get_elevation_glyphs(surface_name, source, lut = vtkLookupTable() color_series.BuildLookupTable(lut, color_series.CATEGORICAL) - lut.SetNanColor(0, 0, 0, 1) - lut.SetTableRange(scalar_range) + lut.nan_color = (0, 0, 0, 1) + lut.table_range = scalar_range if surface_name in ['plane']: lut.number_of_table_values = 1 lut1 = vtkLookupTable() color_series.BuildLookupTable(lut1, color_series.ORDINAL) - lut1.SetNanColor(0, 0, 0, 1) - lut1.SetTableRange(scalar_range) + lut1.nan_color = (0, 0, 0, 1) + lut1.table_range = scalar_range number_of_bands = lut.number_of_table_values lut1.number_of_table_values = number_of_bands @@ -720,9 +722,9 @@ def get_elevation_glyphs(surface_name, source, # Adjust the number of table values and scalar range. scalar_range = (bands[0][0], bands[len(bands) - 1][2]) - lut.TableRange = scalar_range + lut.table_range = scalar_range lut.number_of_table_values = len(bands) - lut1.TableRange = scalar_range + lut1.table_range = scalar_range lut1.number_of_table_values = len(bands) if frequency_table: -- GitLab From 9b078af782bfe6c4dca48155dc607c71768c0080 Mon Sep 17 00:00:00 2001 From: amaclean Date: Fri, 12 Jun 2026 16:58:34 +1000 Subject: [PATCH 5/5] Adding initial version of CurvaturesDemo --- .../Visualization/CurvaturesDemo.md | 16 + .../Visualization/CurvaturesDemo.py | 677 ++++++++++++++++++ 2 files changed, 693 insertions(+) create mode 100644 src/PythonicAPI/Visualization/CurvaturesDemo.md create mode 100644 src/PythonicAPI/Visualization/CurvaturesDemo.py diff --git a/src/PythonicAPI/Visualization/CurvaturesDemo.md b/src/PythonicAPI/Visualization/CurvaturesDemo.md new file mode 100644 index 00000000000..d2c24859f9f --- /dev/null +++ b/src/PythonicAPI/Visualization/CurvaturesDemo.md @@ -0,0 +1,16 @@ +### Description + +**How to get the Gaussian and Mean curvatures of a surface.** + +Two different surfaces are used in this demonstration with each surface coloured according to its Gaussian and Mean curvatures. + +- The first surface is a superquadric surface, this demonstrates the use of extra filters that are needed to get a nice smooth surface. + +- The second surface is a parametric surface, in this case the surface has already been triangulated so no extra processing is necessary. + +In order to get a nice coloured image, a vtkColorTransferFunction has been used to generate a set of colours for the vtkLookupTable tables. We have used a diverging colour space. +Because of the symmetry of the ranges selected for the lookup tables, the white colouration represents a midpoint value whilst the blue represents values less than the midpoint value and orange represents colours greater than the midpoint value. + +In the case of the Random Hills Gaussian curvature surface, this colouration shows the nature of the surface quite nicely. The blue areas are saddle points (negative Gaussian curvature) and the orange areas have a positive Gaussian curvature. + +In the case of the mean curvature the blue colouration represents negative curvature perpendicular to one of the principal axes. diff --git a/src/PythonicAPI/Visualization/CurvaturesDemo.py b/src/PythonicAPI/Visualization/CurvaturesDemo.py new file mode 100644 index 00000000000..03f14537fea --- /dev/null +++ b/src/PythonicAPI/Visualization/CurvaturesDemo.py @@ -0,0 +1,677 @@ +#!/usr/bin/env python3 + +""" +The purpose of this is to demonstrate how to get the Gaussian and Mean curvatures of a surface. + +Two different surfaces are used in this demonstration with each surface coloured according + to its Gaussian and Mean curvatures. + +The first surface is a superquadric surface, this demonstrates the use of extra filters + that are needed to get a nice smooth surface. + +The second surface is a parametric surface, in this case the surface has already been triangulated +so no extra processing is necessary. + +In order to get a nice coloured image, a vtkColorTransferFunction has been used to generate + a set of colours for the vtkLookUp tables. We have used a diverging colour space. +Because of the symmetry of the ranges selected for the lookup tables, the white colouration + represents a midpoint value whilst the blue represents values less than the midopoint value + and orange represents colours greater than the midpoint value. + +In the case of the Random Hills Gaussian Curvature surface, this colouration shows the nature + of the surface quite nicely. The blue areas are saddle points (negative Gaussian curvature) + and the orange areas have a positive Gaussian curvature. +In the case of the mean curvature the blue colouration is representing negative curvature + perpendicular to one of the principal axes. + +This example also demonstrates the use of lists and the linking of the elements of the + lists together to form a pipeline. + +""" +import copy +from dataclasses import dataclass +from pathlib import Path + +import numpy as np +from pip._internal.index import sources +from vtk.util import numpy_support +from vtkmodules.numpy_interface import dataset_adapter as dsa +from vtkmodules.vtkCommonColor import vtkNamedColors +from vtkmodules.vtkCommonComputationalGeometry import vtkParametricRandomHills +from vtkmodules.vtkCommonCore import ( + VTK_DOUBLE, + vtkIdList, + vtkLookupTable +) +from vtkmodules.vtkCommonTransforms import vtkTransform +from vtkmodules.vtkFiltersCore import ( + vtkCleanPolyData, + vtkFeatureEdges, + vtkGenerateIds, + vtkTriangleFilter +) +from vtkmodules.vtkFiltersGeneral import ( + vtkCurvatures, + vtkTransformFilter +) +from vtkmodules.vtkFiltersSources import ( + vtkParametricFunctionSource, + vtkSuperquadricSource +) +from vtkmodules.vtkInteractionStyle import vtkInteractorStyleTrackballCamera +from vtkmodules.vtkInteractionWidgets import vtkScalarBarRepresentation, vtkScalarBarWidget, vtkTextRepresentation, \ + vtkTextWidget +from vtkmodules.vtkRenderingAnnotation import vtkScalarBarActor +from vtkmodules.vtkRenderingCore import ( + vtkActor, + vtkActor2D, + vtkColorTransferFunction, + vtkPolyDataMapper, + vtkRenderWindow, + vtkRenderWindowInteractor, + vtkRenderer, + vtkTextMapper, + vtkTextProperty, vtkTextActor +) + + +def main(argv): + colors = vtkNamedColors() + + # We are going to handle two different sources. + # The first source is a superquadric source. + torus = vtkSuperquadricSource(center=(0.0, 0.0, 0.0), scale=(1.0, 1.0, 1.0), + phi_resolution=64, + theta_resolution=64, theta_roundness=1, + thickness=0.5, size=0.5, toroidal=True) + + # Rotate the torus towards the observer (around the x-axis) + toroid_transform = vtkTransform() + toroid_transform.RotateX(55) + + toroid_transform_filter = vtkTransformFilter(transform=toroid_transform) + + # The quadric is made of strips, so pass it through a triangle filter as + # the curvature filter only operates on polys + tri = vtkTriangleFilter() + + # 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 = vtkCleanPolyData(tolerance=0.005) + + torus >> toroid_transform_filter >> tri >> cleaner + + # The next source will be a parametric function + rh_fn_src = vtkParametricFunctionSource(parametric_function=vtkParametricRandomHills()) + + sources = list() + for i in range(0, 4): + cc = vtkCurvatures() + if i < 2: + cc.input_connection = cleaner.output_port + else: + cc.input_connection = rh_fn_src.output_port + if i % 2 == 0: + cc.SetCurvatureTypeToGaussian() + curvature_name = 'Gauss_Curvature' + else: + cc.SetCurvatureTypeToMean() + curvature_name = 'Mean_Curvature' + cc.update() + adjust_edge_curvatures(cc.output, curvature_name) + sources.append(cc.output) + + curvatures = { + 0: 'Gauss_Curvature', + 1: 'Mean_Curvature', + 2: 'Gauss_Curvature', + 3: 'Mean_Curvature', + } + + # lut = get_diverging_lut() + lut = get_diverging_lut1() + + + # # Create a common text property. + # text_property = vtkTextProperty() + # text_property.SetFontSize(24) + # text_property.SetJustificationToCentered() + # Position the source name according to its length. + curvature_names = list(curvatures.values()) + curvature_names[:] = [x.replace('_','\n') for x in curvature_names] + + text_positions = get_text_positions(curvature_names, + justification=TextProperty.Justification.VTK_TEXT_CENTERED, + vertical_justification=TextProperty.VerticalJustification.VTK_TEXT_BOTTOM, + width=0.5) + + title_text_property = vtkTextProperty(color=colors.GetColor3d('AliceBlue'), bold=True, italic=True, shadow=True, + font_size=16, + justification=TextProperty.Justification.VTK_TEXT_LEFT) + label_text_property = vtkTextProperty(color=colors.GetColor3d('AliceBlue'), bold=False, italic=False, shadow=True, + font_size=12, + justification=TextProperty.Justification.VTK_TEXT_LEFT) + + + # RenderWindow Dimensions + # + renderer_size = 512 + grid_dimensions = 2 + window_width = renderer_size * grid_dimensions + window_height = renderer_size * grid_dimensions + + # Create the RenderWindow and interactor. + # + ren_win = vtkRenderWindow( + size=(renderer_size * grid_dimensions, renderer_size * grid_dimensions), + window_name=f'{Path(argv[0]).stem:s}') + iren = vtkRenderWindowInteractor() + iren.render_window = ren_win + style = vtkInteractorStyleTrackballCamera() + iren.interactor_style = style + + renderers = list() + mappers = list() + actors = list() + text_mappers = list() + text_actors = list() + scalar_bars = list() + sbws = list() + text_representations=list() + text_widgets=list() + + # Link the pipeline together. + for idx, source in enumerate(sources): + curvature_name = curvatures[idx].replace('_', '\n') + + source.point_data.active_scalars = curvatures[idx] + scalar_range = source.point_data.GetScalars(curvatures[idx]).range + + mapper = vtkPolyDataMapper(input_data=source, + lookup_table=lut, + scalar_range=scalar_range, + scalar_mode=Mapper.ScalarMode.VTK_SCALAR_MODE_USE_POINT_FIELD_DATA + ) + mapper.SelectColorArray(curvatures[idx]) + + actor = vtkActor(mapper=mapper) + actors.append(actor) + + ren = vtkRenderer() + # text_mapper = vtkTextMapper(input=curvature_name, text_property=text_property) + # + # text_actor = vtkActor2D(mapper=text_mapper, position=(250, 16)) + # text_actors.append(text_actor) + surface_name = curvature_name.replace('_', '\n') + text_actor = vtkTextActor(input=surface_name.title(), text_scale_mode=vtkTextActor.TEXT_SCALE_MODE_NONE, + text_property=title_text_property) + text_actors.append(text_actor) + # Create the text representation. Used for positioning the text actor. + text_representation = vtkTextRepresentation(enforce_normalized_viewport_bounds=True) + text_representation.position_coordinate.value = text_positions[surface_name]['p'] + text_representation.position2_coordinate.value = text_positions[surface_name]['p2'] + text_representations.append(text_representation) + + text_widget = vtkTextWidget(representation=text_representation, text_actor=text_actor, interactor=iren, + default_renderer = ren, current_renderer = ren, + selectable=False, On=True) + text_widgets.append(text_widget) + + sbp = ScalarBarProperties() + sbp.title_text = surface_name + '\n' + sbp.number_of_labels = 5 + sbp.lut = lut + sbp.orientation = True + sbw = make_scalar_bar_widget(sbp, title_text_property, + label_text_property, ren, iren) + sbws.append(sbw) + + renderers.append(ren) + + for idx in range(len(sources)): + if idx < grid_dimensions * grid_dimensions: + renderers.append(vtkRenderer) + + + # Add and position the renders to the render window. + viewport = list() + for row in range(grid_dimensions): + for col in range(grid_dimensions): + idx = row * grid_dimensions + col + + viewport[:] = [] + viewport.append(float(col) / grid_dimensions) + viewport.append(float(grid_dimensions - (row + 1)) / grid_dimensions) + viewport.append(float(col + 1) / grid_dimensions) + viewport.append(float(grid_dimensions - row) / grid_dimensions) + + if idx > (len(sources) - 1): + continue + + renderers[idx].SetViewport(viewport) + ren_win.AddRenderer(renderers[idx]) + + renderers[idx].AddActor(actors[idx]) + # renderers[idx].AddActor(text_actors[idx]) + # renderers[idx].AddActor(scalar_bars[idx]) + renderers[idx].SetBackground(colors.GetColor3d('ParaViewBlueGrayBkg')) + + # for idx in range(0,len(text_widgets)): + # text_widgets[idx].default_renderer = renderers[idx] + for sbw in sbws: + sbw.On() + for tw in text_widgets: + pass + + ren_win.Render() + + iren.Start() + + +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(color_space=ColorTransferFunction.ColorSpace.VTK_CTF_DIVERGING) + # Cool to warm. + ctf.AddRGBPoint(0.0, 0.230, 0.299, 0.754) + ctf.AddRGBPoint(0.5, 0.865, 0.865, 0.865) + ctf.AddRGBPoint(1.0, 0.706, 0.016, 0.150) + + 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 get_diverging_lut1(): + colors = vtkNamedColors() + + # Colur transfer function. + ctf = vtkColorTransferFunction(color_space=ColorTransferFunction.ColorSpace.VTK_CTF_DIVERGING) + p1 = [0.0] + list(colors.GetColor3d('MidnightBlue')) + p2 = [0.5] + list(colors.GetColor3d('Gainsboro')) + p3 = [1.0] + list(colors.GetColor3d('DarkOrange')) + ctf.AddRGBPoint(*p1) + ctf.AddRGBPoint(*p2) + ctf.AddRGBPoint(*p3) + + table_size = 256 + lut = vtkLookupTable(number_of_table_values=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 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: The vtkPolyData object with the adjusted edge curvatures. + """ + + def point_neighbourhood(pt_id): + """ + Find the ids of the neighbors 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.number_of_ids): + 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.number_of_ids): + 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.point_data.active_scalars = curvature_name + np_source = dsa.WrapDataObject(source) + curvatures = np_source.PointData[curvature_name] + + # Get the boundary point IDs. + array_name = 'ids' + id_filter = vtkGenerateIds(input_data=source, point_ids=True, cell_ids=False, + point_ids_array_name=array_name, cell_ids_array_name=array_name) + + edges = vtkFeatureEdges(boundary_edges=True, manifold_edges=False, + non_manifold_edges=False, feature_edges=False) + + (source >> id_filter >> edges).update() + + edge_array = edges.output.point_data.GetArray(array_name) + boundary_ids = [] + for i in range(edges.output.number_of_points): + 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.name = curvature_name + source.point_data.RemoveArray(curvature_name) + source.point_data.AddArray(curv) + source.point_data.active_scalars = curvature_name + +class ScalarBarProperties: + """ + The properties needed for scalar bars. + """ + named_colors = vtkNamedColors() + + lut = None + # These are in pixels + maximum_dimensions = {'width': 100, 'height': 260} + title_text = '', + number_of_labels: int = 5 + label_format = '{:0.2f}' + # Orientation vertical=True, horizontal=False. + orientation: bool = True + # Horizontal and vertical positioning. + # These are the default positions, don't change these. + default_v = {'p': (0.85, 0.1), 'p2': (0.1, 0.7)} + default_h = {'p': (0.1, 0.1), 'p2': (0.7, 0.1)} + # Modify these as needed. + position_v = copy.deepcopy(default_v) + position_h = copy.deepcopy(default_h) + + +def make_scalar_bar_widget(scalar_bar_properties, title_text_property, label_text_property, renderer, + interactor): + """ + Make a scalar bar widget. + + :param scalar_bar_properties: The lookup table, title name, maximum dimensions in pixels and position. + :param title_text_property: The properties for the title. + :param label_text_property: The properties for the labels. + :param renderer: The default renderer. + :param interactor: The vtkInteractor. + :return: The scalar bar widget. + """ + sb_actor = vtkScalarBarActor(lookup_table=scalar_bar_properties.lut, title=scalar_bar_properties.title_text, + unconstrained_font_size=True, + number_of_labels=scalar_bar_properties.number_of_labels, + title_text_property=title_text_property, label_text_property=label_text_property, + label_format=scalar_bar_properties.label_format, + ) + + sb_rep = vtkScalarBarRepresentation(enforce_normalized_viewport_bounds=True, + orientation=scalar_bar_properties.orientation) + + # Set the position. + sb_rep.position_coordinate.SetCoordinateSystemToNormalizedViewport() + sb_rep.position2_coordinate.SetCoordinateSystemToNormalizedViewport() + if scalar_bar_properties.orientation: + sb_rep.position_coordinate.value = scalar_bar_properties.position_v['p'] + sb_rep.position2_coordinate.value = scalar_bar_properties.position_v['p2'] + else: + sb_rep.position_coordinate.value = scalar_bar_properties.position_h['p'] + sb_rep.position2_coordinate.value = scalar_bar_properties.position_h['p2'] + + widget = vtkScalarBarWidget(representation=sb_rep, scalar_bar_actor=sb_actor, default_renderer=renderer, + interactor=interactor, enabled=True) + + return widget + + +def position_sbw_h(num_bands, max_bands): + """ + Position the vertical scalar bar widget. + :param: num_bands - the number of bands in the scalar bar. + :param: max_bands - the maximum number of bands. + :return: The scalar bar position. + """ + + max_bands = abs(max_bands) + num_bands = abs(num_bands) + if num_bands > max_bands: + num_bands = max_bands + if num_bands == 0: + num_bands = 1 + # Origin of the scalar bar. + xy0 = [0.125, 0.05] + # Width and height of the scalar bar. + dxy = [0.75, 0.1] + if num_bands >= max_bands: + return {'p': tuple(xy0), 'p2': tuple(dxy)} + + dx = dxy[0] - xy0[0] * num_bands / max_bands + dxy[0] = dxy[0] * num_bands / max_bands + if num_bands == 1: + xy0[0] = 0.5 - dx * num_bands / (max_bands * 2) + else: + xy0[0] = 0.5 - dx * (num_bands + 1) / (max_bands * 2) + return {'p': tuple(xy0), 'p2': tuple(dxy)} + + +def position_sbw_v(num_bands, max_bands): + """ + Position the vertical scalar bar widget. + :param: num_bands - the number of bands in the scalar bar. + :param: max_bands - the maximum number of bands. + :return: The scalar bar position. + """ + + max_bands = abs(max_bands) + num_bands = abs(num_bands) + if num_bands > max_bands: + num_bands = max_bands + if num_bands == 0: + num_bands = 1 + # Origin of the scalar bar. + xy0 = [0.9, 0.25] + # Width and height of the scalar bar. + dxy = [0.08, 0.5] + if num_bands >= max_bands: + return {'p': tuple(xy0), 'p2': tuple(dxy)} + + dy = dxy[1] - xy0[1] * num_bands / max_bands + dxy[1] = dxy[1] * num_bands / max_bands + if num_bands == 1: + xy0[1] = 0.5 - dy * num_bands / (max_bands * 2) + else: + xy0[1] = 0.5 - dy * (num_bands + 1) / (max_bands * 2) + return {'p': tuple(xy0), 'p2': tuple(dxy)} + + +def get_text_positions(names, justification=0, vertical_justification=0, width=0.96, height=0.1): + """ + Get viewport positioning information for a list of names. + + :param names: The list of names. + :param justification: Horizontal justification of the text, default is left. + :param vertical_justification: Vertical justification of the text, default is bottom. + :param width: Width of the bounding_box of the text in screen coordinates. + :param height: Height of the bounding_box of the text in screen coordinates. + :return: A dictionary of positioning information. + """ + # The gap between the left or right edge of the screen and the text. + dx = 0.02 + width = abs(width) + if width > 0.96: + width = 0.96 + + y0 = 0.01 + height = abs(height) + if height > 0.9: + height = 0.9 + dy = height + if vertical_justification == TextProperty.VerticalJustification.VTK_TEXT_TOP: + y0 = 1.0 - (dy + y0) + if vertical_justification == TextProperty.VerticalJustification.VTK_TEXT_CENTERED: + y0 = 0.5 - (dy / 2.0 + y0) + + name_len_min = 0 + name_len_max = 0 + first = True + for k in names: + sz = len(k) + if first: + name_len_max = sz + first = False + else: + name_len_min = min(name_len_min, sz) + name_len_max = max(name_len_max, sz) + text_positions = dict() + for k in names: + sz = len(k) + delta_sz = width * sz / name_len_max + if delta_sz > width: + delta_sz = width + + if justification == TextProperty.Justification.VTK_TEXT_CENTERED: + x0 = 0.5 - delta_sz / 2.0 + elif justification == TextProperty.Justification.VTK_TEXT_RIGHT: + x0 = 1.0 - dx - delta_sz + else: + # Default is left justification. + x0 = dx + + # For debugging! + # print( + # f'{k:16s}: (x0, y0) = ({x0:3.2f}, {y0:3.2f}), (x1, y1) = ({x0 + delta_sz:3.2f}, {y0 + dy:3.2f})' + # f', width={delta_sz:3.2f}, height={dy:3.2f}') + text_positions[k] = {'p': [x0, y0, 0], 'p2': [delta_sz, dy, 0]} + + return text_positions + + +@dataclass(frozen=True) +class ColorTransferFunction: + @dataclass(frozen=True) + class ColorSpace: + VTK_CTF_RGB: int = 0 + VTK_CTF_HSV: int = 1 + VTK_CTF_LAB: int = 2 + VTK_CTF_DIVERGING: int = 3 + VTK_CTF_LAB_CIEDE2000: int = 4 + VTK_CTF_STEP: int = 5 + + @dataclass(frozen=True) + class Scale: + VTK_CTF_LINEAR: int = 0 + VTK_CTF_LOG10: int = 1 + +@dataclass(frozen=True) +class Mapper: + @dataclass(frozen=True) + class ColorMode: + VTK_COLOR_MODE_DEFAULT: int = 0 + VTK_COLOR_MODE_MAP_SCALARS: int = 1 + VTK_COLOR_MODE_DIRECT_SCALARS: int = 2 + + @dataclass(frozen=True) + class ResolveCoincidentTopology: + VTK_RESOLVE_OFF: int = 0 + VTK_RESOLVE_POLYGON_OFFSET: int = 1 + VTK_RESOLVE_SHIFT_ZBUFFER: int = 2 + + @dataclass(frozen=True) + class ScalarMode: + VTK_SCALAR_MODE_DEFAULT: int = 0 + VTK_SCALAR_MODE_USE_POINT_DATA: int = 1 + VTK_SCALAR_MODE_USE_CELL_DATA: int = 2 + VTK_SCALAR_MODE_USE_POINT_FIELD_DATA: int = 3 + VTK_SCALAR_MODE_USE_CELL_FIELD_DATA: int = 4 + VTK_SCALAR_MODE_USE_FIELD_DATA: int = 5 + + +@dataclass(frozen=True) +class TextProperty: + @dataclass(frozen=True) + class Justification: + VTK_TEXT_LEFT: int = 0 + VTK_TEXT_CENTERED: int = 1 + VTK_TEXT_RIGHT: int = 2 + + @dataclass(frozen=True) + class VerticalJustification: + VTK_TEXT_BOTTOM: int = 0 + VTK_TEXT_CENTERED: int = 1 + VTK_TEXT_TOP: int = 2 + + +if __name__ == '__main__': + import sys + + main(sys.argv) -- GitLab