diff --git a/src/Cxx.md b/src/Cxx.md index a94b4fea5cf043b216987bdcfdcb2ae7a9044900..fe7fd051f6b08c77a97b6896bc1f817a35a822e9 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 91a5376d9ea22522b049a00833353c094b7eebb9..cb114e93380bea2ec5eeeb64d0abc53f1dab35cc 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 bb8b4c24a6cbb6ba2245a55df5717f994d266411..250d5df30db29ec544d7750a80bf5cbd26799ab8 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/Visualization/CurvatureBandsWithGlyphs.cxx b/src/Cxx/Visualization/CurvatureBandsWithGlyphs.cxx index 64f3f467fad405f7b6d053245a5f34ba98779754..cb6682e852de8377818e4f7f465c9381c68be847 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,100 +82,6 @@ vtkSmartPointer GetSource(std::string const& source); vtkNew ReverseLUT(vtkLookupTable* lut); -struct BandedGlyphs -{ - 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. * @@ -238,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. @@ -245,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; @@ -350,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. * @@ -601,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); @@ -691,1120 +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; - }; + return elevFilter->GetPolyDataOutput(); +} - 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) +vtkSmartPointer GetHills() +{ + // Create four hills on a plane. + // This will have regions of negative, zero and positive Gsaussian curvatures. + + auto xRes = 50; + auto yRes = 50; + auto xMin = -5.0; + auto xMax = 5.0; + auto dx = (xMax - xMin) / (xRes - 1.0); + auto yMin = -5.0; + auto yMax = 5.0; + auto dy = (yMax - yMin) / (xRes - 1.0); + + // Make a grid. + vtkNew points; + for (auto i = 0; i < xRes; ++i) { - curvatures.push_back(array->GetVariantValue(i).ToDouble()); + auto x = xMin + i * dx; + for (auto j = 0; j < yRes; ++j) + { + auto y = yMin + j * dy; + points->InsertNextPoint(x, y, 0); + } } - // 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(); + // Add the grid points to a polydata object. + vtkNew plane; + plane->SetPoints(points); - vtkNew edges; + // Triangulate the grid. + vtkNew delaunay; + delaunay->SetInputData(plane); + delaunay->Update(); - edges->SetInputConnection(idFilter->GetOutputPort()); - edges->BoundaryEdgesOn(); - edges->ManifoldEdgesOff(); - edges->NonManifoldEdgesOff(); - edges->FeatureEdgesOff(); - edges->Update(); + auto polydata = delaunay->GetOutput(); - 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) + 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) { - 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 = polydata->GetPoint(i); + for (size_t j = 0; j < hd.size(); ++j) { - // Corner case. - // countInvalid += 1; - // Assuming the curvature of the point is planar. - newCurv = 0.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); } - // Set the new curvature value. - curvatures[pId] = newCurv; + polydata->GetPoints()->SetPoint(i, x); + elevation->SetValue(i, x[2]); } - // Set small values to zero. - if (epsilon != 0.0) + vtkNew textures; + textures->SetNumberOfComponents(2); + textures->SetNumberOfTuples(2 * polydata->GetNumberOfPoints()); + textures->SetName("Textures"); + + for (auto i = 0; i < xRes; ++i) { - auto eps = std::abs(epsilon); - for (size_t i = 0; i < curvatures.size(); ++i) + float tc[2]; + tc[0] = i / (xRes - 1.0); + for (auto j = 0; j < yRes; ++j) { - if (std::abs(curvatures[i]) < eps) - { - curvatures[i] = 0.0; - } + // tc[1] = 1.0 - j / (yRes - 1.0); + tc[1] = j / (yRes - 1.0); + textures->SetTuple(static_cast(i) * yRes + j, tc); } } - 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()); + 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(); } -void ConstrainCurvatures(vtkPolyData* source, std::string const& curvatureName, - double const& lowerBound, double const& upperBound) +vtkSmartPointer GetParametricHills() { - std::array bounds{0.0, 0.0}; - if (lowerBound < upperBound) - { - bounds[0] = lowerBound; - bounds[1] = upperBound; - } - else - { - bounds[0] = upperBound; - bounds[1] = lowerBound; - } + vtkNew fn; + fn->AllowRandomGenerationOn(); + fn->SetRandomSeed(1); + fn->SetNumberOfHills(30); - 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) - { - if (curvatures[i] < bounds[0]) - { - curvatures[i] = bounds[0]; - } - else - { - if (curvatures[i] > bounds[1]) - { - curvatures[i] = bounds[1]; - } - } - } - 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; - } - - source->GetPointData()->SetActiveScalars(curvature.c_str()); - auto sr = source->GetPointData()->GetScalars(curvature.c_str())->GetRange(); - std::array scalarRange{sr[0], sr[1]}; - - 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 - { - colorSeries->SetColorScheme(colorSeries->BREWER_DIVERGING_SPECTRAL_5); - } + vtkNew source; + source->SetParametricFunction(fn); + source->SetUResolution(50); + source->SetVResolution(50); + source->SetScalarModeToZ(); + source->Update(); - vtkNew lut; - colorSeries->BuildLookupTable(lut, colorSeries->CATEGORICAL); - lut->SetNanColor(0, 0, 0, 1); - lut->SetTableRange(scalarRange.data()); + // 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 lut1; - colorSeries->BuildLookupTable(lut1, colorSeries->ORDINAL); - lut1->SetNanColor(0, 0, 0, 1); - lut1->SetTableRange(scalarRange.data()); + 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(); - auto numberOfBands = lut->GetNumberOfTableValues(); - lut1->SetNumberOfTableValues(numberOfBands); + return transformFilter->GetPolyDataOutput(); +} - auto bands = GetBands(scalarRange, numberOfBands, precision, nearestInteger); +vtkSmartPointer GetParametricTorus() +{ + vtkNew fn; + fn->SetRingRadius(5); + fn->SetCrossSectionRadius(2); - 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}}; + vtkNew source; + source->SetParametricFunction(fn); + source->SetUResolution(50); + source->SetVResolution(50); + source->SetScalarModeToZ(); + source->Update(); - // 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}, - }; + // 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"); - // Comment this out if you want to see how allocating - // equally spaced bands works. - bands = GetCustomBands(scalarRange, numberOfBands, myBands); - } - } + vtkNew transform; + transform->RotateX(-90.0); + vtkNew transformFilter; + transformFilter->SetInputConnection(source->GetOutputPort()); + transformFilter->SetTransform(transform); + transformFilter->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()); + return transformFilter->GetPolyDataOutput(); +} - if (frequencyTable) - { - auto freq = GetFrequencies(bands, source); - AdjustFrequencyRanges(bands, freq); - fmt::println("{:s} {:s}", Title(surfaceName), curvature); - PrintBandsFrequencies(bands, freq); - } +vtkSmartPointer GetPlane() +{ + 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(); - // 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])); - } + 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(); - // 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()); - } + // 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()); - // 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(); + // 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(); - // Generate the glyphs on the original surface. - auto glyphs = GetGlyphs(source, scaleFactor, false); + return cleaner->GetOutput(); +} - BandedGlyphs bg; - 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; +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 bg; + return source->GetOutput(); } -BandedGlyphs GetElevationGlyphs(std::string const& surfaceName, - vtkPolyData* source, unsigned int precision, - bool frequencyTable, bool nearestInteger) +vtkSmartPointer GetTorus() { - // 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->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(); + + // 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()); + + // 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(); - source->GetPointData()->SetActiveScalars("Elevation"); - auto sr = source->GetPointData()->GetScalars("Elevation")->GetRange(); - std::array scalarRange{sr[0], sr[1]}; + return cleaner->GetOutput(); +} - vtkNew colorSeries; - if (surfaceName == "hills" || 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") { - colorSeries->SetColorScheme( - colorSeries->BREWER_DIVERGING_BROWN_BLUE_GREEN_8); + return GetHills(); } - else + else if (surface == "parametric hills") { - colorSeries->SetColorScheme( - colorSeries->BREWER_DIVERGING_BROWN_BLUE_GREEN_5); + return GetParametricHills(); } - - 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") + else if (surface == "parametric torus") { - // 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 GetParametricTorus(); } - - // 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 == "plane") { - auto freq = GetFrequencies(bands, source); - AdjustFrequencyRanges(bands, freq); - fmt::println("{:s} Elevation", Title(surfaceName)); - PrintBandsFrequencies(bands, freq); + return GenerateElevations(GetPlane()); } - - // 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 == "sphere") { - labels.push_back(fmt::format("{:4.2f}", p->second[1])); + return GenerateElevations(GetSphere()); } - - // Annotate - vtkNew values; - for (size_t i = 0; i < labels.size(); ++i) + else if (surface == "torus") { - values->InsertNextValue(vtkVariant(labels[i])); + return GenerateElevations(GetTorus()); } - for (vtkIdType i = 0; i < values->GetNumberOfTuples(); ++i) + std::cout << "The surface is not available." << std::endl; + std::cout << "Using parametric hills instead." << std::endl; + return GetParametricHills(); +} + +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) { - lut->SetAnnotation(i, values->GetValue(i).ToString()); + 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()); } - - // 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) + t = lut->GetNumberOfAnnotatedValues() - 1; + for (vtkIdType i = t; i >= 0; --i) { - bcf->SetValue(i, p->second[2]); - ++i; + lutr->SetAnnotation(t - i, lut->GetAnnotation(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.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; + 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()); + 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; + }; - // 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 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; + }; - return cleaner->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()); + } -vtkSmartPointer GetSphere() -{ - vtkNew source; - source->SetCenter(0.0, 0.0, 0.0); - source->SetRadius(1.0); - source->SetThetaResolution(32); - source->SetPhiResolution(32); - source->Update(); + // 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 source->GetOutput(); -} + vtkNew edges; -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(); + edges->SetInputConnection(idFilter->GetOutputPort()); + edges->BoundaryEdgesOn(); + edges->ManifoldEdgesOff(); + edges->NonManifoldEdgesOff(); + edges->FeatureEdgesOff(); + edges->Update(); - // 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 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; + } - // 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(); + // 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; + } + } + } - return cleaner->GetOutput(); + 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()); } -vtkSmartPointer GetSource(std::string const& source) +void ConstrainCurvatures(vtkPolyData* source, std::string const& curvatureName, + double const& lowerBound, double const& upperBound) { - std::string surface = source; - std::transform(surface.begin(), surface.end(), surface.begin(), - [](unsigned char c) { return std::tolower(c); }); - if (surface == "hills") - { - return GetHills(); - } - else if (surface == "parametric hills") + std::array bounds{0.0, 0.0}; + if (lowerBound < upperBound) { - return GetParametricHills(); + bounds[0] = lowerBound; + bounds[1] = upperBound; } - else if (surface == "parametric torus") + else { - return GetParametricTorus(); + bounds[0] = upperBound; + bounds[1] = lowerBound; } - else if (surface == "plane") + + 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(GetPlane()); + curvatures.push_back(array->GetVariantValue(i).ToDouble()); } - else if (surface == "sphere") + // Set upper and lower bounds. + for (size_t i = 0; i < curvatures.size(); ++i) { - return GenerateElevations(GetSphere()); + if (curvatures[i] < bounds[0]) + { + curvatures[i] = bounds[0]; + } + else + { + if (curvatures[i] > bounds[1]) + { + curvatures[i] = bounds[1]; + } + } } - else if (surface == "torus") + vtkNew adjustedCurvatures; + for (auto curvature : curvatures) { - return GenerateElevations(GetTorus()); + adjustedCurvatures->InsertNextTuple1(curvature); } - std::cout << "The surface is not available." << std::endl; - std::cout << "Using parametric hills instead." << std::endl; - return GetParametricHills(); + 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) + colorSeries->SetColorScheme(colorSeries->BREWER_DIVERGING_SPECTRAL_6); + maxBands = 6; + } + + vtkNew lut; + colorSeries->BuildLookupTable(lut, colorSeries->CATEGORICAL); + lut->SetNanColor(0, 0, 0, 1); + lut->SetTableRange(scalarRange.data()); + + if (surfaceName == "plane" || surfaceName == "sphere") + { + 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 (curvature == "Gaussian_Curvature") + { + if (surfaceName == "random hills") { - for (std::vector::iterator p = b.begin(); p != b.end(); ++p) - { - *p = RoundOff(*p); - } - b[0] = x[0]; + // 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); } - bands[i] = b; - for (std::vector::iterator p = b.begin(); p != b.end(); ++p) + else if (surfaceName == "hills") { - *p = RoundOff(*p + dx); + 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); } } - 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)) + // 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) { - return bands; + auto freq = GetFrequencies(bands, source); + AdjustRanges(bands, freq); + fmt::println("{:s} {:s}", Title(surfaceName), curvature); + PrintBandsFrequencies(bands, freq); } - std::vector> x; - std::copy(myBands.begin(), myBands.end(), std::back_inserter(x)); + // 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])); + } - // Determine the index of the range minimum and range maximum. - int idxMin = 0; - for (auto idx = 0; idx < static_cast(myBands.size()); ++idx) + // Annotate + vtkNew values; + for (size_t i = 0; i < labels.size(); ++i) { - if (scalarRange[0] < myBands[idx][1] && scalarRange[0] >= myBands[idx][0]) - { - idxMin = idx; - break; - } + values->InsertNextValue(vtkVariant(labels[i])); } - int idxMax = static_cast(myBands.size()) - 1; - for (int idx = static_cast(myBands.size()) - 1; idx >= 0; --idx) + for (vtkIdType i = 0; i < values->GetNumberOfTuples(); ++i) { - if (scalarRange[1] < myBands[idx][1] && scalarRange[1] >= myBands[idx][0]) - { - idxMax = static_cast(idx); - break; - } + lut->SetAnnotation(i, values->GetValue(i).ToString()); } - // 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) + // 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) { - 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; + bcf->SetValue(i, p->second[2]); + ++i; } - return bands; + // 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> -GetIntegralBands(std::array const& scalarRange) +BandedGlyphs GetElevationGlyphs(std::string const& surfaceName, + vtkPolyData* source, unsigned int precision, + bool frequencyTable, bool nearestInteger) { - if (scalarRange[1] < scalarRange[0]) + // The length of the normal arrow glyphs. + auto scaleFactor = 1.0; + if (surfaceName == "hills") { - std::map> bands; - return bands; + scaleFactor = 0.5; } - std::array x{0, 0}; - for (int i = 0; i < 2; ++i) + if (surfaceName == "sphere") { - x[i] = scalarRange[i]; + scaleFactor = 0.2; } - 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); -} -std::map GetFrequencies(std::map>& bands, - vtkPolyData* src) -{ - std::map freq; - for (auto i = 0; i < static_cast(bands.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") { - freq[i] = 0; + colorSeries->SetColorScheme( + colorSeries->BREWER_DIVERGING_BROWN_BLUE_GREEN_8); + maxBands = 8; } - vtkIdType tuples = src->GetPointData()->GetScalars()->GetNumberOfTuples(); - for (int i = 0; i < tuples; ++i) + else { - 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; - } - } + colorSeries->SetColorScheme( + colorSeries->BREWER_DIVERGING_BROWN_BLUE_GREEN_6); + maxBands = 6; } - 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) + vtkNew lut; + colorSeries->BuildLookupTable(lut, colorSeries->CATEGORICAL); + lut->SetNanColor(0, 0, 0, 1); + lut->SetTableRange(scalarRange.data()); + + if (surfaceName == "plane") { - if (freq[i] != 0) - { - first = i; - break; - } + lut->SetNumberOfTableValues(1); } - std::vector keys; - for (std::map::iterator it = freq.begin(); it != freq.end(); ++it) + 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") { - keys.push_back(it->first); + // 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::reverse(keys.begin(), keys.end()); - auto last = keys[0]; - for (size_t i = 0; i < keys.size(); ++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) { - if (freq[keys[i]] != 0) - { - last = keys[i]; - break; - } + auto freq = GetFrequencies(bands, source); + AdjustRanges(bands, freq); + fmt::println("{:s} Elevation", Title(surfaceName)); + PrintBandsFrequencies(bands, freq); } - // 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) + + // 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) { - adjFreq[idx] = p.second; - ++idx; + labels.push_back(fmt::format("{:4.2f}", p->second[1])); } - std::map> adjBands; - idx = 0; - for (auto const& p : bands) + + // Annotate + vtkNew values; + for (size_t i = 0; i < labels.size(); ++i) { - adjBands[idx] = p.second; - ++idx; + values->InsertNextValue(vtkVariant(labels[i])); } - 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; - - if (bands.size() != freq.size()) + for (vtkIdType i = 0; i < values->GetNumberOfTuples(); ++i) { - std::cout << "Bands and frequencies must be the same size." << std::endl; - return; + lut->SetAnnotation(i, values->GetValue(i).ToString()); } - std::ostringstream os; - os << "Bands & Frequencies:\n"; - size_t idx = 0; - auto total = 0; - auto width = prec + 6; - std::string res{""}; + + // 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 @@ -2006,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/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/Cxx/Visualization/ElevationBandsWithGlyphs.cxx b/src/Cxx/Visualization/ElevationBandsWithGlyphs.cxx index dedca86b48b35798f5941c7572e36e8845d24bcb..af41d8c560c30f38a22385c05517ef52340849f2 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,46 +74,6 @@ vtkSmartPointer GetSource(std::string const& source); vtkSmartPointer ReverseLUT(vtkLookupTable* lut); -struct BandedGlyphs -{ - 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. * @@ -182,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; @@ -287,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. * @@ -476,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); @@ -534,137 +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.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}; @@ -1169,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; @@ -1227,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 @@ -1471,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/Python.md b/src/Python.md index 10174d7df757306f1a5f8f5ed9248c1a2e7a00b6..8eb1ae74fe7fe7af86b06f9ddde05216f742ffb9 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 5081f1fb643e3cb5122aa45eb937b2c39dae7149..59856f07380f132055c793bb3011dec94c98f997 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/CurvaturesNormalsElevations.py b/src/PythonicAPI/PolyData/CurvaturesNormalsElevations.py deleted file mode 100644 index 40bc7d2b8e22675a4b9e477185379cabae120952..0000000000000000000000000000000000000000 --- a/src/PythonicAPI/PolyData/CurvaturesNormalsElevations.py +++ /dev/null @@ -1,1397 +0,0 @@ -#!/usr/bin/env python3 - -import copy -import math -from collections import namedtuple, OrderedDict -from dataclasses import dataclass -from pathlib import Path - -import numpy as np -# noinspection PyUnresolvedReferences -import vtkmodules.vtkInteractionStyle -# noinspection PyUnresolvedReferences -import vtkmodules.vtkRenderingOpenGL2 -from vtk.util import numpy_support -from vtkmodules.numpy_interface import dataset_adapter as dsa -from vtkmodules.vtkCommonColor import ( - vtkColorSeries, - vtkNamedColors -) -from vtkmodules.vtkCommonComputationalGeometry import ( - vtkParametricRandomHills, - vtkParametricTorus -) -from vtkmodules.vtkCommonCore import ( - VTK_DOUBLE, - vtkFloatArray, - vtkIdList, - vtkLookupTable, - vtkPoints, - vtkVariant, - vtkVariantArray -) -from vtkmodules.vtkCommonDataModel import ( - vtkCellArray, - vtkTriangle -) -from vtkmodules.vtkCommonDataModel import vtkPolyData -from vtkmodules.vtkCommonTransforms import vtkTransform -from vtkmodules.vtkFiltersCore import ( - vtkCleanPolyData, - vtkElevationFilter, - vtkFeatureEdges, - vtkGlyph3D, - vtkGenerateIds, - vtkMaskPoints, - vtkPolyDataNormals, - vtkReverseSense, - vtkTriangleFilter -) -from vtkmodules.vtkFiltersGeneral import ( - vtkCurvatures, - vtkTransformFilter -) -from vtkmodules.vtkFiltersModeling import vtkBandedPolyDataContourFilter -from vtkmodules.vtkFiltersSources import ( - vtkArrowSource, - vtkParametricFunctionSource, - vtkPlaneSource, - vtkSphereSource, - vtkSuperquadricSource -) -from vtkmodules.vtkInteractionStyle import vtkInteractorStyleTrackballCamera -from vtkmodules.vtkInteractionWidgets import ( - vtkCameraOrientationWidget, - vtkOrientationMarkerWidget, - vtkScalarBarRepresentation, - vtkScalarBarWidget, - vtkTextRepresentation, - vtkTextWidget -) -from vtkmodules.vtkRenderingAnnotation import vtkAxesActor, vtkScalarBarActor -from vtkmodules.vtkRenderingCore import ( - vtkActor, - vtkColorTransferFunction, - vtkPolyDataMapper, - vtkRenderWindow, - vtkRenderWindowInteractor, - vtkRenderer, - vtkTextActor, - vtkTextProperty -) - - -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 - 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', - 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', - help='Use an OrientationMarkerWidget instead of a CameraOrientationWidget.') - - args = parser.parse_args() - return args.surface_name, args.frequency_table, args.omw - - -def main(argv): - surface_name, frequency_table, use_camera_omw = get_program_parameters() - - available_surfaces = ['hills', 'parametric torus', 'plane', 'random hills', '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'] - - surface_name = ' '.join(surface_name.lower().replace('_', ' ').split()) - if surface_name == 'randomhills': - surface_name = 'random hills' - if surface_name == 'parametrictorus': - surface_name = 'parametric torus' - if surface_name.lower() not in available_surfaces: - print('Nonexistent surface:', surface_name) - print('Available surfaces are:') - asl = sorted(available_surfaces) - asl = [asl[i].title() for i in range(0, len(asl))] - asl = [asl[i:i + 5] for i in range(0, len(asl), 5)] - for i in range(0, len(asl)): - s = ', '.join(asl[i]) - 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"') - return - - Surface = namedtuple('Surface', 'name source') - surface = Surface(surface_name, get_source(surface_name, available_surfaces)) - - # -------------------------------------------------------------------------------------- - # 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) - - 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() - viewports['Gauss_Curvature'] = [0.0, 0.0, 0.5, 1.0] - viewports['Mean_Curvature'] = [0.5, 0.0, 1.0, 1.0] - - window_height = 800 - window_width = 2 * window_height - - # -------------------------------------------------- - # Create the RenderWindow, Renderers and Interactor. - # -------------------------------------------------- - ren_win = vtkRenderWindow(size=(window_width, window_height), window_name=f'{Path(argv[0]).name:s}') - iren = vtkRenderWindowInteractor() - iren.render_window = ren_win - style = vtkInteractorStyleTrackballCamera() - iren.interactor_style = style - - 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) - text_actor = vtkTextActor(input=surface_name.title(), text_scale_mode=vtkTextActor.TEXT_SCALE_MODE_NONE, - text_property=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_widget = vtkTextWidget(representation=text_representation, text_actor=text_actor, interactor=iren, - selectable=False) - - first = True - 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) - - 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_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, - scalar_visibility=True, - color_mode=Mapper.ColorMode.VTK_COLOR_MODE_MAP_SCALARS) - glyph_mapper.SelectColorArray('Elevation') - - glyph_actor = vtkActor(mapper=glyph_mapper) - v['glyph'] >> glyph_mapper - - renderer = 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) - - 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) - - for k in curvatures.keys(): - if k == 'Gauss_Curvature': - contour_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: - 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, - interactive=True) - - 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): - """ - - :param source: The name of the source. - :param available_surfaces: The surfaces - :return: - """ - 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 - - -def get_hills(): - """ - Create four hills on a plane. - This will have regions of negative, zero and positive Gaussian curvatures. - - :return: - """ - - x_res = 50 - y_res = 50 - x_min = -5.0 - x_max = 5.0 - dx = (x_max - x_min) / (x_res - 1) - y_min = -5.0 - y_max = 5.0 - dy = (y_max - y_min) / (x_res - 1) - - # Make a grid. - # We define the parameters for the hills here. - # There are four hills. - # [[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 - - # 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) - - for i in range(0, x_res): - tc = [i / (x_res - 1.0), 0.0] - for j in range(0, y_res): - # tc[1] = 1.0 - j / (y_res - 1.0) - tc[1] = j / (y_res - 1.0) - textures.SetTuple(i * y_res + j, tc) - - bounds = poly_data.bounds - elevation_filter = vtkElevationFilter(low_point=(0.0, 0.0, bounds[4]), high_point=(0.0, 0.0, bounds[5])) - - normals = vtkPolyDataNormals(feature_angle=30, splitting=False) - - tr = vtkTransform() - tr.RotateX(-90.0) - tf = vtkTransformFilter(transform=tr) - - return poly_data >> elevation_filter >> normals >> tf - - -def get_parametric_hills(): - fn = vtkParametricRandomHills(random_seed=1, number_of_hills=30) - fn.AllowRandomGenerationOn() - - source = vtkParametricFunctionSource(parametric_function=fn, u_resolution=50, v_resolution=50, - scalar_mode=vtkParametricFunctionSource.SCALAR_Z) - source.SetScalarModeToZ() - src = source.update().output - - # Rename the scalars to 'Elevation' since we are using the Z-scalars as elevations. - src.point_data.scalars.SetName('Elevation') - - transform = vtkTransform() - transform.Translate(0.0, 5.0, 15.0) - transform.RotateX(-90.0) - transform_filter = vtkTransformFilter(transform=transform) - - return src >> transform_filter - - -def get_parametric_torus(): - fn = vtkParametricTorus(ring_radius=5, cross_section_radius=2) - - source = vtkParametricFunctionSource(parametric_function=fn, u_resolution=50, v_resolution=50, - scalar_mode=vtkParametricFunctionSource.SCALAR_Z) - src = source.update().output - - # Rename the scalars to 'Elevation' since we are using the Z-scalars as elevations. - src.point_data.scalars.SetName('Elevation') - - transform = vtkTransform() - transform.Translate(0.0, 0.0, 0.0) - transform.RotateX(-90.0) - transform_filter = vtkTransformFilter(transform=transform) - - return src >> transform_filter - - -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 - - transform = vtkTransform() - transform.Translate(0.0, 0.0, 0.0) - transform.RotateX(-90.0) - - transform_filter = vtkTransformFilter(transform=transform) - - # We have an 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() - - # Pass it though a CleanPolyDataFilter and merge any points which - # 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 - - -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 - - 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 - - -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 - - # 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) - - 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 - - -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 reverse_lut(lut): - """ - Create a lookup table with the colors reversed. - :param: lut - An indexed lookup table. - :return: The reversed indexed lookup table. - """ - lutr = vtkLookupTable() - lutr.DeepCopy(lut) - t = lut.GetNumberOfTableValues() - 1 - rev_range = reversed(list(range(t + 1))) - for i in rev_range: - rgba = [0.0] * 3 - v = float(i) - lut.GetColor(v, rgba) - rgba.append(lut.GetOpacity(v)) - lutr.SetTableValue(t - i, rgba) - 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)) - return lutr - - -def get_glyphs(surface, arrow_scale=None, scale_factor=None, reverse_normals=False): - """ - Glyph 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) - - # 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. - if reverse_normals: - reverse = vtkReverseSense(reverse_cells=True, reverse_normals=True) - source >> reverse >> mask_pts - else: - source >> 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 - - 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() - - return mask_pts >> glyph - - -def get_bands(d_r, number_of_bands, precision=2, nearest_integer=False): - """ - Divide a range into bands. - - :param: d_r - [min, max] the range that is to be covered by the bands. - :param: number_of_bands - The number of bands, a positive integer. - :param: precision - The decimal precision of the bounds. - :param: nearest_integer - If True then [floor(min), ceil(max)] is used. - :return: A dictionary consisting of the band number and [min, midpoint, max] for each band. - """ - prec = abs(precision) - if prec > 14: - prec = 14 - - bands = dict() - if (d_r[1] < d_r[0]) or (number_of_bands <= 0): - return bands - x = list(d_r) - if nearest_integer: - x[0] = math.floor(x[0]) - x[1] = math.ceil(x[1]) - dx = (x[1] - x[0]) / float(number_of_bands) - b = [x[0], x[0] + dx / 2.0, x[0] + dx] - i = 0 - while i < number_of_bands: - b = list(map(lambda ele_b: round(ele_b, prec), b)) - if i == 0: - b[0] = x[0] - bands[i] = b - b = [b[0] + dx, b[1] + dx, b[2] + dx] - i += 1 - return bands - - -def get_custom_bands(d_r, number_of_bands, my_bands): - """ - Divide a range into custom bands. - - You need to specify each band as 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: number_of_bands - the number of bands, a positive integer. - :return: A dictionary consisting of band number and [min, midpoint, max] for each band. - """ - bands = dict() - if (d_r[1] < d_r[0]) or (number_of_bands <= 0): - return bands - x = my_bands - # Determine the index of the range minimum and range maximum. - idx_min = 0 - for idx in range(0, len(my_bands)): - if my_bands[idx][1] > d_r[0] >= my_bands[idx][0]: - idx_min = idx - break - - idx_max = len(my_bands) - 1 - for idx in range(len(my_bands) - 1, -1, -1): - if my_bands[idx][1] > d_r[1] >= my_bands[idx][0]: - idx_max = idx - break - - # Set the minimum to match the range minimum. - x[idx_min][0] = d_r[0] - x[idx_max][1] = d_r[1] - x = x[idx_min: idx_max + 1] - for idx, e in enumerate(x): - bands[idx] = [e[0], e[0] + (e[1] - e[0]) / 2, e[1]] - return bands - - -def get_frequencies(bands, src): - """ - Count the number of scalars in each band. - The scalars used are the active scalars in the polydata. - - :param: bands - The bands. - :param: src - The vtkPolyData source. - :return: The frequencies of the scalars in each band. - """ - freq = dict() - for i in range(len(bands)): - freq[i] = 0 - tuples = src.GetPointData().GetScalars().GetNumberOfTuples() - for i in range(tuples): - x = src.GetPointData().GetScalars().GetTuple1(i) - for j in range(len(bands)): - if x <= bands[j][2]: - freq[j] += 1 - break - return freq - - -def adjust_ranges(bands, freq): - """ - The bands and frequencies are adjusted so that the first and last - frequencies in the range are non-zero. - - :param bands: The bands dictionary. - :param freq: The frequency dictionary. - :return: Adjusted bands and frequencies. - """ - # Get the indices of the first and last non-zero frequency elements. - first = 0 - for k, v in freq.items(): - if v != 0: - first = k - break - rev_keys = list(freq.keys())[::-1] - last = rev_keys[0] - for idx in list(freq.keys())[::-1]: - if freq[idx] != 0: - last = idx - break - # Now adjust the ranges. - min_key = min(freq.keys()) - max_key = max(freq.keys()) - for idx in range(min_key, first): - freq.pop(idx) - bands.pop(idx) - for idx in range(last + 1, max_key + 1): - freq.popitem() - bands.popitem() - old_keys = freq.keys() - adj_freq = dict() - adj_bands = dict() - - for idx, k in enumerate(old_keys): - adj_freq[idx] = freq[k] - adj_bands[idx] = bands[k] - - return adj_bands, adj_freq - - -def print_bands_frequencies(curvature, bands, freq, precision=2): - prec = abs(precision) - if prec > 14: - prec = 14 - - if len(bands) != len(freq): - print('Bands and Frequencies must be the same size.') - 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) - - -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}') - - # 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()) - - -def generate_gaussian_curvatures(surface, needs_adjusting, frequency_table=False): - """ - Generate the filters for calculating Gaussian curvatures on the surface. - - :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)) - - 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) - - curv_lut.SetTableRange(curv_scalar_range) - curv_lut.SetNumberOfTableValues(len(curv_bands)) - - annotate_lut(curv_lut, curv_bands) - - # 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) - - # 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} - - -def generate_mean_curvatures(surface, needs_adjusting, frequency_table=False): - """ - Generate the filters for calculating mean curvatures on the surface. - - :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) - - curv_lut.SetTableRange(curv_scalar_range) - curv_lut.SetNumberOfTableValues(len(curv_bands)) - - annotate_lut(curv_lut, curv_bands) - - # 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) - - # 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} - - -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 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 BandedPolyDataContourFilter: - @dataclass(frozen=True) - class ScalarMode: - VTK_SCALAR_MODE_INDEX: int = 0 - VTK_SCALAR_MODE_VALUE: int = 1 - - -@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 Curvatures: - @dataclass(frozen=True) - class CurvatureType: - VTK_CURVATURE_GAUSS: int = 0 - VTK_CURVATURE_MEAN: int = 1 - VTK_CURVATURE_MAXIMUM: int = 2 - VTK_CURVATURE_MINIMUM: int = 3 - - -@dataclass(frozen=True) -class Glyph3D: - @dataclass(frozen=True) - class ColorMode: - VTK_COLOR_BY_SCALE: int = 0 - VTK_COLOR_BY_SCALAR: int = 1 - VTK_COLOR_BY_VECTOR: int = 2 - - @dataclass(frozen=True) - class IndexMode: - VTK_INDEXING_OFF: int = 0 - VTK_INDEXING_BY_SCALAR: int = 1 - VTK_INDEXING_BY_VECTOR: int = 2 - - @dataclass(frozen=True) - class ScaleMode: - VTK_SCALE_BY_SCALAR: int = 0 - VTK_SCALE_BY_VECTOR: int = 1 - VTK_SCALE_BY_VECTORCOMPONENTS: int = 2 - VTK_DATA_SCALING_OFF: int = 3 - - @dataclass(frozen=True) - class VectorMode: - VTK_USE_VECTOR: int = 0 - VTK_USE_NORMAL: int = 1 - VTK_VECTOR_ROTATION_OFF: int = 2 - VTK_FOLLOW_CAMERA_DIRECTION: int = 3 - - -@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) diff --git a/src/PythonicAPI/Visualization/CurvatureBandsWithGlyphs.py b/src/PythonicAPI/Visualization/CurvatureBandsWithGlyphs.py index 99cc5d54f333c3db33c5d0aa4ba1a5ae9986fe2f..316ef18eabd51f97d91c5335f6b6f3350fddbbcd 100644 --- a/src/PythonicAPI/Visualization/CurvatureBandsWithGlyphs.py +++ b/src/PythonicAPI/Visualization/CurvatureBandsWithGlyphs.py @@ -2,7 +2,9 @@ import copy import math +from collections import namedtuple from dataclasses import dataclass +from pathlib import Path import numpy as np # noinspection PyUnresolvedReferences @@ -201,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 @@ -241,36 +243,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) @@ -296,154 +288,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. @@ -821,6 +665,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.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 + + +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.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) + 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 @@ -881,28 +873,39 @@ 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 + 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) + 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 @@ -942,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: @@ -983,11 +986,19 @@ 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_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 @@ -996,18 +1007,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): """ @@ -1030,21 +1029,28 @@ 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) + 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 @@ -1069,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: @@ -1104,11 +1110,19 @@ 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_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 @@ -1117,18 +1131,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/PolyData/Curvatures.py b/src/PythonicAPI/Visualization/Curvatures.py old mode 100755 new mode 100644 similarity index 91% rename from src/PythonicAPI/PolyData/Curvatures.py rename to src/PythonicAPI/Visualization/Curvatures.py index aeaf1788771ff6372e1969ebc1fc2dd27d4da228..d7d74df03dacb834289e953cba6426b07e69ece9 --- a/src/PythonicAPI/PolyData/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/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 96% rename from src/PythonicAPI/PolyData/CurvaturesAdjustEdges.py rename to src/PythonicAPI/Visualization/CurvaturesAdjustEdges.py index 81558568e04a0430d323f46c05bee03e79de8616..3b2c4af3fcac385dca29602dd336d972cd7431af --- a/src/PythonicAPI/PolyData/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/CurvaturesDemo.md b/src/PythonicAPI/Visualization/CurvaturesDemo.md new file mode 100644 index 0000000000000000000000000000000000000000..d2c24859f9fdc9b97e46e9a286b554aced02806f --- /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 0000000000000000000000000000000000000000..03f14537fead055e6d647ef9cecddf32e195644e --- /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) 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/Visualization/CurvaturesNormalsElevations.py b/src/PythonicAPI/Visualization/CurvaturesNormalsElevations.py new file mode 100644 index 0000000000000000000000000000000000000000..82db6e20a866314f832db957e24104673656d306 --- /dev/null +++ b/src/PythonicAPI/Visualization/CurvaturesNormalsElevations.py @@ -0,0 +1,1531 @@ +#!/usr/bin/env python3 + +import copy +import math +from collections import namedtuple, OrderedDict +from dataclasses import dataclass +from pathlib import Path + +import numpy as np +# noinspection PyUnresolvedReferences +import vtkmodules.vtkInteractionStyle +# noinspection PyUnresolvedReferences +import vtkmodules.vtkRenderingOpenGL2 +from vtk.util import numpy_support +from vtkmodules.numpy_interface import dataset_adapter as dsa +from vtkmodules.vtkCommonColor import ( + vtkColorSeries, + vtkNamedColors +) +from vtkmodules.vtkCommonComputationalGeometry import ( + vtkParametricRandomHills, + vtkParametricTorus +) +from vtkmodules.vtkCommonCore import ( + VTK_DOUBLE, + vtkDoubleArray, + vtkFloatArray, + vtkIdList, + vtkLookupTable, + vtkPoints, + vtkVariant, + vtkVariantArray +) +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 +) +from vtkmodules.vtkFiltersGeneral import ( + vtkCurvatures, + vtkTransformFilter +) +from vtkmodules.vtkFiltersModeling import vtkBandedPolyDataContourFilter +from vtkmodules.vtkFiltersSources import ( + vtkArrowSource, + vtkParametricFunctionSource, + vtkPlaneSource, + vtkSphereSource, + vtkSuperquadricSource +) +from vtkmodules.vtkInteractionStyle import vtkInteractorStyleTrackballCamera +from vtkmodules.vtkInteractionWidgets import ( + vtkCameraOrientationWidget, + vtkOrientationMarkerWidget, + vtkScalarBarRepresentation, + vtkScalarBarWidget, + vtkTextRepresentation, + vtkTextWidget +) +from vtkmodules.vtkRenderingAnnotation import vtkAxesActor, vtkScalarBarActor +from vtkmodules.vtkRenderingCore import ( + vtkActor, + vtkPolyDataMapper, + vtkRenderWindow, + vtkRenderWindowInteractor, + vtkRenderer, + vtkTextActor, + vtkTextProperty +) + + +def get_program_parameters(): + import argparse + description = 'Demonstrates Gaussian and Mean curvatures on a surface, along with normals colored by elevation.' + epilogue = ''' + 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='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('-o', '--omw', action='store_true', + help='Use an OrientationMarkerWidget instead of a CameraOrientationWidget.') + + args = parser.parse_args() + return args.surface_name, args.frequency_table, args.omw + + +def main(argv): + surface_name, frequency_table, use_omw = get_program_parameters() + + 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 hills', 'parametric torus', 'plane'] + + surface_name = ' '.join(surface_name.lower().replace('_', ' ').split()) + 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: + print('Nonexistent surface:', surface_name) + print('Available surfaces are:') + asl = sorted(available_surfaces) + asl = [asl[i].title() for i in range(0, len(asl))] + asl = [asl[i:i + 5] for i in range(0, len(asl), 5)] + for i in range(0, len(asl)): + s = ', '.join(asl[i]) + if i < len(asl) - 1: + s += ',' + print(f' {s}') + print('If a name has spaces in it, delineate the name with quotes e.g. "parametric hills"') + return + + source = get_source(surface_name) + if not source: + print('The surface is not available.') + return + + # Use an ordered dictionary as we want the keys in a specific order. + curvatures = OrderedDict() + # 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() + + # Define viewport ranges [x_min, y_min, x_max, y_max] + viewports = dict() + viewports['Gauss_Curvature'] = [0.0, 0.0, 0.5, 1.0] + viewports['Mean_Curvature'] = [0.5, 0.0, 1.0, 1.0] + + window_height = 800 + window_width = 2 * window_height + + # -------------------------------------------------- + # Create the RenderWindow, Renderers and Interactor. + # -------------------------------------------------- + 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() + iren.interactor_style = style + + renderers = list() + contour_widgets = dict() + elevation_widgets = dict() + + # 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.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=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_widget = vtkTextWidget(representation=text_representation, text_actor=text_actor, interactor=iren, + selectable=False) + + # Scalar bar properties for the curvatures and elevation. + curv_sbps = dict() + for k, v in curvatures.items(): + 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) + + # Create contour edges + 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') + + 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, + scalar_mode=Mapper.ScalarMode.VTK_SCALAR_MODE_USE_POINT_FIELD_DATA) + + glyph_mapper.SelectColorArray('Elevation') + glyph_actor = vtkActor(mapper=glyph_mapper) + + ren = vtkRenderer(background=colors.GetColor3d('ParaViewBlueGrayBkg')) + + contour_widgets[k] = make_scalar_bar_widget(curv_sbps[k], title_text_property, + label_text_property, ren, iren) + contour_widgets[k].Off() + + 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() + + # 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=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() + + ren_win.Render() + iren.Start() + + +def generate_elevations(src): + """ + Generate elevations over the surface. + :param: src - the vtkPolyData source. + :return: - vtkPolyData source with elevations. + """ + 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. + + x_res = 50 + y_res = 50 + x_min = -5.0 + x_max = 5.0 + dx = (x_max - x_min) / (x_res - 1) + y_min = -5.0 + y_max = 5.0 + dy = (y_max - y_min) / (x_res - 1) + + # Make a grid. + points = vtkPoints() + for i in range(0, x_res): + x = x_min + i * dx + for j in range(0, y_res): + y = y_min + j * dy + points.InsertNextPoint(x, y, 0) + + # Add the grid points to a polydata object. + plane = vtkPolyData(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]] + 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]) + + 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] + for j in range(0, y_res): + # tc[1] = 1.0 - j / (y_res - 1.0) + tc[1] = j / (y_res - 1.0) + textures.SetTuple(i * y_res + j, tc) + + polydata.GetPointData().SetScalars(elevation) + polydata.GetPointData().GetScalars().name = 'Elevation' + polydata.GetPointData().TCoords = textures + + normals = vtkPolyDataNormals(feature_angle=30, splitting=False) + + tr = vtkTransform() + tr.RotateX(-90) + + tf = vtkTransformFilter(transform=tr) + + return (polydata >> normals >> tf).update().output + + +def get_parametric_hills(): + """ + 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() + + 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() + source.update() + # Rename the scalars to 'Elevation' since we are using the Z-scalars as elevations. + source.output.GetPointData().GetScalars().SetName('Elevation') + + # Build the tangents. + tangents = vtkPolyDataTangents() + + tr = vtkTransform() + tr.Translate(0.0, 5.0, 15.0) + tr.RotateX(-90.0) + + tf = vtkTransformFilter(transform=tr) + + return (source >> tangents >> tf).update().output + + +def get_parametric_torus(): + """ + Make a parametric torus as the source. + :return: vtkPolyData with normal and scalar data. + """ + + 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. + 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) + + t = vtkTransformFilter(transform=tr) + + return (source >> tangents >> t).update().output + + +def get_plane(): + """ + 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) + + tr = vtkTransform() + tr.Translate(0.0, 0.0, 0.0) + tr.RotateX(-90.0) + + tf = vtkTransformFilter(transform=tr) + + # 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() + + # Pass it though a CleanPolyDataFilter and merge any points which + # are coincident, or very close + cleaner = vtkCleanPolyData(tolerance=0.005) + + return (source >> tf >> tri >> cleaner).update().output + + +def get_sphere(): + source = vtkSphereSource(center=(0.0, 0.0, 0.0), radius=1.0, + theta_resolution=32, phi_resolution=32) + + return source.update().output + + +def get_torus(): + """ + 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 + 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) + + return (source >> tri >> cleaner).update().output + + +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): + """ + Create a lookup table with the colors reversed. + :param: lut - An indexed lookup table. + :return: The reversed indexed lookup table. + """ + lutr = vtkLookupTable() + lutr.DeepCopy(lut) + t = lut.number_of_table_values - 1 + rev_range = reversed(list(range(t + 1))) + for i in rev_range: + rgba = [0.0] * 3 + v = float(i) + lut.GetColor(v, rgba) + rgba.append(lut.GetOpacity(v)) + lutr.SetTableValue(t - i, rgba) + t = lut.number_of_annotated_values - 1 + rev_range = reversed(list(range(t + 1))) + for i in rev_range: + lutr.SetAnnotation(t - i, lut.GetAnnotation(i)) + return lutr + + +def get_glyphs(src, scale_factor=1.0, reverse_normals=False): + """ + Glyph the normals on the surface. + + You may need to adjust the parameters for mask_pts, arrow and glyph for a + nice appearance. + + :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) + # Choose a random subset of points. + mask_pts = vtkMaskPoints(on_ratio=5, random_mode=True) + src >> reverse >> mask_pts + else: + # 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(tip_resolution=16, tip_length=0.3, tip_radius=0.1) + + # 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 + + +def get_bands(scalar_range, number_of_bands, precision=2, nearest_integer=False): + """ + 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. + :return: A dictionary consisting of the band number and [min, midpoint, max] for each band. + """ + prec = abs(precision) + if prec > 14: + prec = 14 + + bands = dict() + if (scalar_range[1] < scalar_range[0]) or (number_of_bands <= 0): + return bands + x = list(scalar_range) + if nearest_integer: + x[0] = math.floor(x[0]) + x[1] = math.ceil(x[1]) + dx = (x[1] - x[0]) / float(number_of_bands) + b = [x[0], x[0] + dx / 2.0, x[0] + dx] + i = 0 + while i < number_of_bands: + b = list(map(lambda ele_b: round(ele_b, prec), b)) + if i == 0: + b[0] = x[0] + bands[i] = b + b = [b[0] + dx, b[1] + dx, b[2] + dx] + i += 1 + return 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]...] + + :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 (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] > 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] > scalar_range[1] >= my_bands[idx][0]: + idx_max = idx + break + + # Set the minimum to match the range minimum. + 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]] + return bands + + +def get_frequencies(bands, src): + """ + Count the number of scalars in each band. + The scalars used are the active scalars in the polydata. + + :param: bands - The bands. + :param: src - The vtkPolyData source. + :return: The frequencies of the scalars in each band. + """ + freq = dict() + for i in range(len(bands)): + freq[i] = 0 + tuples = src.GetPointData().GetScalars().GetNumberOfTuples() + for i in range(tuples): + x = src.GetPointData().GetScalars().GetTuple1(i) + for j in range(len(bands)): + if x <= bands[j][2]: + freq[j] += 1 + break + return freq + + +def adjust_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 + + +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.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) + 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 dictionary containing the bands. + :param freq: The frequency dictionary. + :return: Adjusted bands and frequencies. + """ + # Get the indices of the first and last non-zero elements. + first = 0 + for k, v in freq.items(): + if v != 0: + first = k + break + rev_keys = list(freq.keys())[::-1] + last = rev_keys[0] + for idx in list(freq.keys())[::-1]: + if freq[idx] != 0: + last = idx + break + # Now adjust the ranges. + min_key = min(freq.keys()) + max_key = max(freq.keys()) + for idx in range(min_key, first): + freq.pop(idx) + bands.pop(idx) + for idx in range(last + 1, max_key + 1): + freq.popitem() + bands.popitem() + old_keys = freq.keys() + adj_freq = dict() + adj_bands = dict() + + for idx, k in enumerate(old_keys): + adj_freq[idx] = freq[k] + adj_bands[idx] = bands[k] + + return adj_bands, adj_freq + + +def 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 + + """ + # 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.point_data.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.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.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 + + 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.table_range = scalar_range + lut.number_of_table_values = len(bands) + lut1.table_range = 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) + + # 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 + 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=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): + """ + 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 + + """ + # 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.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_6 + max_lut_bands = 6 + + lut = vtkLookupTable() + color_series.BuildLookupTable(lut, color_series.CATEGORICAL) + 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.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 + + 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) + + # Adjust the number of table values and scalar range. + scalar_range = (bands[0][0], bands[len(bands) - 1][2]) + lut.table_range = scalar_range + lut.number_of_table_values = len(bands) + lut1.table_range = scalar_range + lut1.number_of_table_values = len(bands) + + if frequency_table: + 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) + + # 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 + 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=source, + 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 = '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: + """ + 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 + + +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) + class ScalarMode: + VTK_SCALAR_MODE_INDEX: int = 0 + VTK_SCALAR_MODE_VALUE: int = 1 + + +@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 Curvatures: + @dataclass(frozen=True) + class CurvatureType: + VTK_CURVATURE_GAUSS: int = 0 + VTK_CURVATURE_MEAN: int = 1 + VTK_CURVATURE_MAXIMUM: int = 2 + VTK_CURVATURE_MINIMUM: int = 3 + + +@dataclass(frozen=True) +class Glyph3D: + @dataclass(frozen=True) + class ColorMode: + VTK_COLOR_BY_SCALE: int = 0 + VTK_COLOR_BY_SCALAR: int = 1 + VTK_COLOR_BY_VECTOR: int = 2 + + @dataclass(frozen=True) + class IndexMode: + VTK_INDEXING_OFF: int = 0 + VTK_INDEXING_BY_SCALAR: int = 1 + VTK_INDEXING_BY_VECTOR: int = 2 + + @dataclass(frozen=True) + class ScaleMode: + VTK_SCALE_BY_SCALAR: int = 0 + VTK_SCALE_BY_VECTOR: int = 1 + VTK_SCALE_BY_VECTORCOMPONENTS: int = 2 + VTK_DATA_SCALING_OFF: int = 3 + + @dataclass(frozen=True) + class VectorMode: + VTK_USE_VECTOR: int = 0 + VTK_USE_NORMAL: int = 1 + VTK_VECTOR_ROTATION_OFF: int = 2 + VTK_FOLLOW_CAMERA_DIRECTION: int = 3 + + +@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) diff --git a/src/PythonicAPI/Visualization/ElevationBandsWithGlyphs.py b/src/PythonicAPI/Visualization/ElevationBandsWithGlyphs.py index 4dd06aa74303e88ce9901569649e264e09b181f3..9ef2e3d9da1aed062e397d822b33621c76c3a884 100644 --- a/src/PythonicAPI/Visualization/ElevationBandsWithGlyphs.py +++ b/src/PythonicAPI/Visualization/ElevationBandsWithGlyphs.py @@ -2,7 +2,9 @@ import copy import math +from collections import namedtuple from dataclasses import dataclass +from pathlib import Path # noinspection PyUnresolvedReferences import vtkmodules.vtkRenderingOpenGL2 @@ -163,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 @@ -194,20 +197,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) @@ -672,26 +671,33 @@ 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) + 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 @@ -716,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: @@ -751,11 +757,19 @@ 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_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 @@ -764,18 +778,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/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 0000000000000000000000000000000000000000..5a7af1583d464816b115f4db0ac13bfe382312ad --- /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