CurvatureBandsWithGlyphs
Repository source: CurvatureBandsWithGlyphs
Description¶
In this example we are coloring the surface by partitioning the gaussian curvature into bands and using arrows to display the normals on the surface.
Rather beautiful surfaces are generated.
The banded contour filter and an indexed lookup table are used to generate the curvature bands on the surface. To further enhance the surface, the surface normals are glyphed and colored by elevation using a diverging lookup table.
Note that:
-
If the regions on a surface have zero Gaussian curvature, then they can be flattened into a plane with no distortion, and the geometry of the region is Euclidean geometry.
-
If the regions on a surface have positive Gaussian curvature, then the geometry of the surface is spherical geometry.
-
If the regions on the surface have a negative Gaussian curvature, then the geometry of the surface is hyperbolic geometry.
In the above image you can see that the random hills incorporate all of these geometries.
The surface selected is the parametric random hills surface. The problem with the random hills surface is:
- Most of the gaussian curvatures will lie in the range -1 to 0.2 (say) with a few large values say 20 to 40 at the peaks of the hills.
- The edges of the random hills surface also have large irregular values so we need to handle these also. In order to fix this, a function is provided to adjust the edges.
So we need to manually generate custom bands to group the curvatures. The bands selected in the examples show that the surface is mostly planar with some hyperbolic regions (saddle points) and some spherical regions.
Feel free to experiment with different color schemes and/or the other sources from the parametric function group or the torus etc.
You will usually need to adjust the parameters for maskPts, arrow and glyph for a nice appearance.
A histogram of the frequencies is also output to the console. This is useful if you want to get an idea of the distribution of the scalars in each band.
Other languages
See (Python)
Question
If you have a question about this example, please use the VTK Discourse Forum
Code¶
CurvatureBandsWithGlyphs.cxx
#include <vtkActor.h>
#include <vtkArrowSource.h>
#include <vtkBandedPolyDataContourFilter.h>
#include <vtkCamera.h>
#include <vtkCellArray.h>
#include <vtkCleanPolyData.h>
#include <vtkColorSeries.h>
#include <vtkColorTransferFunction.h>
#include <vtkCurvatures.h>
#include <vtkDelaunay2D.h>
#include <vtkDoubleArray.h>
#include <vtkElevationFilter.h>
#include <vtkFeatureEdges.h>
#include <vtkFloatArray.h>
#include <vtkGlyph3D.h>
#include <vtkIdList.h>
#include <vtkInteractorStyleTrackballCamera.h>
#include <vtkLookupTable.h>
#include <vtkMaskPoints.h>
#include <vtkNamedColors.h>
#include <vtkNew.h>
#include <vtkParametricFunctionSource.h>
#include <vtkParametricRandomHills.h>
#include <vtkParametricTorus.h>
#include <vtkPlaneSource.h>
#include <vtkPointData.h>
#include <vtkPolyData.h>
#include <vtkPolyDataMapper.h>
#include <vtkPolyDataNormals.h>
#include <vtkProperty.h>
#include <vtkRenderWindow.h>
#include <vtkRenderWindowInteractor.h>
#include <vtkRenderer.h>
#include <vtkReverseSense.h>
#include <vtkScalarBarActor.h>
#include <vtkScalarBarRepresentation.h>
#include <vtkScalarBarWidget.h>
#include <vtkSmartPointer.h>
#include <vtkSphereSource.h>
#include <vtkSuperquadricSource.h>
#include <vtkSystemIncludes.h>
#include <vtkTextActor.h>
#include <vtkTextProperty.h>
#include <vtkTextRepresentation.h>
#include <vtkTextWidget.h>
#include <vtkTransform.h>
#include <vtkTransformPolyDataFilter.h>
#include <vtkTriangle.h>
#include <vtkTriangleFilter.h>
#include <vtkVariant.h>
#include <vtkVariantArray.h>
#include <vtkVersion.h>
#include <vtk_cli11.h>
#include <vtk_fmt.h>
// clang-format off
#include VTK_FMT(fmt/format.h)
// clang-format on
#if VTK_VERSION_NUMBER >= 90020210809ULL
#define HAS_COW
#include <vtkCameraOrientationWidget.h>
#endif
// vtkGenerateIds was introduced in VTK build version 20240504
#if VTK_BUILD_VERSION >= 20240504
#define USE_USE_GENERATE_IDS
#include <vtkGenerateIds.h>
#else
#include <vtkIdFilter.h>
#endif
// #include <algorithm>
// #include <array>
// #include <cctype>
// #include <cmath>
// #include <cstdlib>
// #include <cstring>
// #include <filesystem>
// #include <functional>
// #include <iomanip>
// #include <iostream>
// #include <iterator>
// #include <numeric>
// #include <set>
// #include <sstream>
// #include <string>
// #include <vector>
namespace fs = std::filesystem;
namespace {
struct CurvatureData
{
// Holds the necessary data to build the view.
vtkSmartPointer<vtkBandedPolyDataContourFilter> bcf;
vtkSmartPointer<vtkGlyph3D> glyph;
std::array<double, 2> curvScalarRange;
std::array<double, 2> elevScalarRange;
vtkSmartPointer<vtkLookupTable> curvLut;
vtkSmartPointer<vtkLookupTable> elevLut;
int elevLabels = 5;
};
/** Generate the filters for calculating Gaussian curvatures on 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 - True if a frequency table is to be displayed.
* @return Return the filters, scalar ranges of curvatures and elevation
* along with the lookup tables.
*/
CurvatureData
GenerateGaussianCurvatures(std::string const& surfaceName, vtkPolyData* source,
std::vector<std::string> const& needsAdjusting,
bool frequencyTable = false);
/** Generate the filters for calculating mean curvatures on 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 - True if a frequency table is to be displayed.
* @return Return the filters, scalar ranges of curvatures and elevation
* along with the lookup tables.
*/
CurvatureData
GenerateMeanCurvatures(std::string const& surfaceName, vtkPolyData* source,
std::vector<std::string> const& needsAdjusting,
bool frequencyTable = false);
/** 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);
/** Generate elevations over the surface.
* @param src - the vtkPolyData source.
* @return elev - the elevations.
*/
vtkSmartPointer<vtkPolyData> GetElevations(vtkPolyData* src);
vtkSmartPointer<vtkPolyData> GetHills();
vtkSmartPointer<vtkPolyData> GetParametricHills();
vtkSmartPointer<vtkPolyData> GetParametricTorus();
vtkSmartPointer<vtkPolyData> GetPlane();
vtkSmartPointer<vtkPolyData> GetSphere();
vtkSmartPointer<vtkPolyData> GetTorus();
vtkSmartPointer<vtkPolyData> GetSource(std::string const& source);
// vtkNew<vtkColorSeries> GetColorSeries();
vtkNew<vtkLookupTable> GetCategoricalLUT();
vtkNew<vtkLookupTable> GetOrdinaLUT();
vtkNew<vtkLookupTable> GetDivergingLUT();
vtkNew<vtkLookupTable> ReverseLUT(vtkLookupTable* lut);
/** Glyph the normals on the surface.
*
* @param name - the name of the vtkPolyData source.
* @param src - the vtkPolyData source.
* @param arrowScale - scaling for the arrows, default is [1, 1, 1].
* @param scaleFactor - the scale factor for the glyphs.
* @param reverseNormals - if True the normals on the surface are reversed.
* @return The glyphs.
*/
vtkNew<vtkGlyph3D> GetGlyphs(std::string const& name, vtkPolyData* src,
std::array<double, 3> const& arrowScale,
double const& scaleFactor = 1.0,
bool const& reverseNormals = false);
std::map<int, std::vector<double>> GetBands(double const dR[2],
int const& numberOfBands,
int const& precision = 2,
bool const& nearestInteger = false);
/** Divide a range into custom bands
*
* You need to specify each band as an array [r1, r2] where r1 < r2 and append
these to a vector.
* The vector should ultimately look like this: [[r1, r2], [r2, r3], [r3,
r4]...]
*
* @param dR - [min, max] the range that is to be covered by the bands.
* @param numberOfBands - the number of bands, a positive integer.
* @param myBands - the bands.
* @return A map consisting of the band inxex and [min, midpoint, max] for
each band.
*/
std::map<int, std::vector<double>>
GetCustomBands(double const dR[2], int const& numberOfBands,
std::vector<std::array<double, 2>> const& myBands);
/** Divide a range into integral bands
*
* Divide a range into bands
*
* @param dR - [min, max] the range that is to be covered by the bands.
* @return A map consisting of the band inxex and [min, midpoint, max] for each
* band.
*/
std::map<int, std::vector<double>> GetIntegralBands(double const dR[2]);
/** 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.
*/
std::map<int, int> GetFrequencies(std::map<int, std::vector<double>>& bands,
vtkPolyData* src);
/** Adjust the frequency ranges.
*
* The bands and frequencies are adjusted so that the first and last
* frequencies in the range are non-zero.
*
* @param bands - The bands.
* @param freq - The frequencies.
*/
void AdjustFrequencyRanges(std::map<int, std::vector<double>>& bands,
std::map<int, int>& freq);
void PrintBandsFrequencies(std::map<int, std::vector<double>> const& bands,
std::map<int, int>& freq, int const& precision = 2);
void AnnotateLut(vtkLookupTable* lut, std::map<int, std::vector<double>> bands);
typedef std::map<std::string, std::array<double, 2>> TTextPosition;
typedef std::map<std::string, TTextPosition> TTextPositions;
struct ScalarBarProperties
{
vtkSmartPointer<vtkNamedColors> colors;
// The properties needed for scalar bars.
vtkSmartPointer<vtkLookupTable> lut;
// These are in pixels.
std::map<std::string, int> maximumDimensions{{"width", 100}, {"height", 260}};
std::string titleText{""};
int number_of_labels{5};
std::string labelFormat{"{:0.2f}"};
// Orientation vertical=true, horizontal=false.
bool orientation{true};
// Horizontal and vertical positioning.
// These are the defaults, don't change these.
TTextPosition defaultV = {{"p", {0.85, 0.1}}, {"p2", {0.1, 0.7}}};
TTextPosition defaultH = {{"p", {0.10, 0.1}}, {"p2", {0.7, 0.1}}};
// Modify these as needed.
TTextPosition positionV = {{"p", {0.85, 0.1}}, {"p2", {0.1, 0.7}}};
TTextPosition positionH = {{"p", {0.10, 0.1}}, {"p2", {0.7, 0.1}}};
};
/** Make a scalar bar widget.
*
* @param scalar_bar_properties - The lookup table, title name, maximum
* dimensions in pixels and position.
* @param textProperty - The properties for the title.
* @param labelTextProperty - The properties for the labels.
* @param ren - The vtkRenderer.
* @param iren - The vtkInteractor.
* @return The scalar bar widget.
*/
vtkNew<vtkScalarBarWidget>
MakeScalarBarWidget(ScalarBarProperties& scalarBarProperties,
vtkTextProperty* textProperty,
vtkTextProperty* labelTextProperty, vtkRenderer* ren,
vtkRenderWindowInteractor* iren);
/**
* Get viewport positioning information for a vector of names.
*
* Note: You must include vtkSystemIncludes.h to get these defines:
* VTK_TEXT_LEFT 0, VTK_TEXT_CENTERED 1, VTK_TEXT_RIGHT 2,
* VTK_TEXT_BOTTOM 0, VTK_TEXT_TOP 2
*
* @param names - The vector of names.
* @param justification - Horizontal justification of the text.
* @param vertical_justification - Vertical justification of the text.
* @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 map of positioning information.
*/
TTextPositions
GetTextPositions(std::vector<std::string> const& names,
int const justification = VTK_TEXT_LEFT,
int const vertical_justification = VTK_TEXT_BOTTOM,
double const width = 0.96, double const height = 0.1);
/** Convert a string to title case.
*
* @param input - The string to convert.
* @return The string converted to title case.
*/
std::string Title(const std::string& input);
} // namespace
int main(int argc, char* argv[])
{
CLI::App app{"Color a surface using curvatures, adding normal vectors "
"colored by elevation."};
// Define options.
std::string surfaceName{""};
app.add_option(
"surfaceName", surfaceName,
"The name of the surface - enclose the name in quotes if it has spaces.");
auto meanCurvature{false};
app.add_flag("-m, --meanCurvature", meanCurvature,
"Use mean curvature instead of gaussian curvature.");
auto frequencyTable{false};
app.add_flag("-f, --frequencyTable", frequencyTable,
"Display the frequency table.");
CLI11_PARSE(app, argc, argv);
if (surfaceName.empty())
{
surfaceName = "RandomHills";
}
// Remove all whitespace characters.
surfaceName.erase(
std::remove_if(surfaceName.begin(), surfaceName.end(), ::isspace),
surfaceName.end());
// Convert to lowercase
std::transform(surfaceName.begin(), surfaceName.end(), surfaceName.begin(),
[](unsigned char c) { return std::tolower(c); });
if (surfaceName == "randomhills")
{
surfaceName = "random hills";
}
if (surfaceName == "parametrictorus")
{
surfaceName = "parametric torus";
}
// Get the surface
std::vector<std::string> availableSurfaces{
"hills", "parametric torus", "plane", "random hills", "sphere", "torus",
};
// Surfaces whose curvatures need to be adjusted along the edges of the
// surface or constrained.
std::vector<std::string> needsAdjusting{"hills", "parametric torus", "plane",
"random hills", "torus"};
auto it = std::find(availableSurfaces.begin(), availableSurfaces.end(),
surfaceName);
if (it == availableSurfaces.end())
{
std::string res{"The surface is not available,"
" available surfaces are:\n "};
for (const auto& element : availableSurfaces)
{
res += fmt::format("{:s}, ", element);
}
if (res.length() >= 2)
{
auto pos = res.length() - 2;
res.replace(pos, 2, "");
}
std::cout << res << std::endl;
return EXIT_FAILURE;
}
// Set this to false if you want mean curvatures.
auto gaussianCurvature = !meanCurvature;
std::string curvature =
(meanCurvature) ? "Mean_Curvature" : "Gauss_Curvature";
auto source = GetSource(surfaceName);
std::map<std::string, CurvatureData> curvatures;
if (gaussianCurvature)
{
curvatures[curvature] = GenerateGaussianCurvatures(
surfaceName, source, needsAdjusting, frequencyTable);
}
else
{
curvatures[curvature] = GenerateMeanCurvatures(
surfaceName, source, needsAdjusting, frequencyTable);
}
// ------------------------------------------------------------
// Create the mappers and actors
// ------------------------------------------------------------
vtkNew<vtkNamedColors> colors;
colors->SetColor("ParaViewBlueGrayBkg",
std::array<unsigned char, 4>{84, 89, 109, 255}.data());
colors->SetColor("ParaViewWarmGrayBkg",
std::array<unsigned char, 4>{98, 93, 90, 255}.data());
vtkNew<vtkPolyDataMapper> srcMapper;
srcMapper->SetInputConnection(curvatures[curvature].bcf->GetOutputPort());
srcMapper->SetScalarRange(curvatures[curvature].curvScalarRange.data());
srcMapper->SetLookupTable(curvatures[curvature].curvLut);
srcMapper->SetScalarModeToUseCellData();
vtkNew<vtkActor> srcActor;
srcActor->SetMapper(srcMapper);
// Create contour edges
vtkNew<vtkPolyDataMapper> edgeMapper;
edgeMapper->SetInputData(curvatures[curvature].bcf->GetContourEdgesOutput());
edgeMapper->SetResolveCoincidentTopologyToPolygonOffset();
vtkNew<vtkActor> edgeActor;
edgeActor->SetMapper(edgeMapper);
edgeActor->GetProperty()->SetColor(colors->GetColor3d("Black").GetData());
vtkNew<vtkPolyDataMapper> glyphMapper;
glyphMapper->SetInputConnection(curvatures[curvature].glyph->GetOutputPort());
glyphMapper->SetScalarModeToUsePointFieldData();
glyphMapper->SetColorModeToMapScalars();
glyphMapper->ScalarVisibilityOn();
glyphMapper->SelectColorArray("Elevation");
// Colour by scalars.
glyphMapper->SetLookupTable(curvatures[curvature].elevLut);
glyphMapper->SetScalarRange(curvatures[curvature].elevScalarRange.data());
vtkNew<vtkActor> glyphActor;
glyphActor->SetMapper(glyphMapper);
auto windowWidth = 800;
auto windowHeight = 800;
vtkNew<vtkTextProperty> textProperty;
textProperty->SetColor(colors->GetColor3d("AliceBlue").GetData());
textProperty->SetJustificationToLeft();
textProperty->SetFontSize(16);
textProperty->BoldOn();
textProperty->ItalicOn();
textProperty->ShadowOn();
std::string curvatureType = curvature;
// std::replace(curvatureType.begin(), curvatureType.end(), '_', '\n');
// ------------------------------------------------------------
// Create the RenderWindow, Renderer and Interactor
// ------------------------------------------------------------
vtkNew<vtkRenderer> ren;
vtkNew<vtkRenderWindow> renWin;
vtkNew<vtkRenderWindowInteractor> iren;
vtkNew<vtkInteractorStyleTrackballCamera> style;
iren->SetInteractorStyle(style);
renWin->AddRenderer(ren);
// Important: The interactor must be set prior to enabling the widget.
iren->SetRenderWindow(renWin);
#ifdef HAS_COW
vtkNew<vtkCameraOrientationWidget> cow;
cow->SetParentRenderer(ren);
// Enable the widget.
cow->On();
#endif
// Add actors.
ren->AddViewProp(srcActor);
ren->AddViewProp(edgeActor);
ren->AddViewProp(glyphActor);
ren->SetBackground(colors->GetColor3d("ParaViewWarmGrayBkg").GetData());
renWin->SetSize(windowWidth, windowHeight);
auto appFn = fs::path((app.get_name())).stem().string();
renWin->SetWindowName(appFn.c_str());
auto namePositions =
GetTextPositions(availableSurfaces, VTK_TEXT_LEFT, VTK_TEXT_TOP, 0.45);
vtkNew<vtkTextProperty> textProp;
// textProp->SetFontFamilyToCourier();
textProp->SetFontSize(16);
textProp->SetVerticalJustificationToTop();
textProp->SetJustificationToLeft();
textProp->BoldOn();
textProp->ItalicOn();
textProp->ShadowOn();
textProp->SetColor(colors->GetColor3d("AliceBlue").GetData());
auto name = Title(surfaceName);
vtkNew<vtkTextActor> textActor;
textActor->SetInput(name.c_str());
textActor->SetTextScaleModeToNone();
textActor->SetTextProperty(textProp);
TTextPosition pos = namePositions[surfaceName];
vtkNew<vtkTextRepresentation> textRep;
textRep->EnforceNormalizedViewportBoundsOn();
textRep->SetPosition(pos["p"][0], pos["p"][1]);
textRep->SetPosition2(pos["p2"][0], pos["p2"][1]);
vtkNew<vtkTextWidget> textWidget;
textWidget->SetRepresentation(textRep);
textWidget->SetTextActor(textActor);
textWidget->SetInteractor(iren);
textWidget->SelectableOff();
textWidget->ResizableOn();
textWidget->On();
// Add scalar bars.
auto sbProperties = ScalarBarProperties();
if (surfaceName == "plane" || surfaceName == "sphere")
{
sbProperties.positionH = {{"p", {0.3, 0.1}}, {"p2", {0.35, 0.1}}};
}
sbProperties.lut = curvatures[curvature].curvLut;
sbProperties.orientation = false;
if (gaussianCurvature)
{
sbProperties.titleText = "Gaussian Curvature\n";
}
else
{
sbProperties.titleText = "Mean Curvature\n";
}
sbProperties.number_of_labels =
curvatures[curvature].curvLut->GetNumberOfTableValues();
auto sbWidgetCurv =
MakeScalarBarWidget(sbProperties, textProperty, textProperty, ren, iren);
if (surfaceName == "plane" || surfaceName == "sphere")
{
sbProperties.positionV = sbProperties.defaultV;
}
if (surfaceName == "plane")
{
sbProperties.positionV = {{"p", {0.85, 0.5}}, {"p2", {0.1, 0.1}}};
}
sbProperties.lut = curvatures[curvature].elevLut;
sbProperties.orientation = true;
sbProperties.titleText = "Elevation\n";
sbProperties.number_of_labels = curvatures[curvature].elevLabels;
auto sbWidgetElev =
MakeScalarBarWidget(sbProperties, textProperty, textProperty, ren, iren);
if (surfaceName == "plane")
{
sbProperties.positionV = sbProperties.defaultV;
}
auto camera = ren->GetActiveCamera();
camera->Elevation(60);
// This moves the window center slightly to ensure that
// the whole surface is not obscured by the scalar bars.
camera->SetWindowCenter(0.0, -0.15);
ren->ResetCamera();
if (surfaceName == "plane")
{
camera->Zoom(0.8);
}
renWin->Render();
iren->Start();
return EXIT_SUCCESS;
}
namespace {
CurvatureData
GenerateGaussianCurvatures(std::string const& surfaceName, vtkPolyData* source,
std::vector<std::string> const& needsAdjusting,
bool frequencyTable)
{
std::string curvature = "Gauss_Curvature";
vtkNew<vtkCurvatures> cc;
cc->SetInputData(source);
cc->SetCurvatureTypeToGaussian();
cc->Update();
if (std::find(needsAdjusting.begin(), needsAdjusting.end(), surfaceName) !=
needsAdjusting.end())
{
AdjustEdgeCurvatures(cc->GetOutput(), curvature);
}
if (surfaceName == "plane")
{
ConstrainCurvatures(cc->GetOutput(), curvature, 0.0, 0.0);
}
if (surfaceName == "sphere")
{
// Gaussian curvature is 1/r^2
auto radius = 10;
auto gauss_curvature = 1.0 / (radius * radius);
ConstrainCurvatures(cc->GetOutput(), curvature, gauss_curvature,
gauss_curvature);
}
cc->GetOutput()->GetPointData()->SetActiveScalars("Elevation");
auto elevSR =
cc->GetOutput()->GetPointData()->GetScalars("Elevation")->GetRange();
std::array<double, 2> elevScalarRange{elevSR[0], elevSR[1]};
// auto elevLut = GetOrdinaLUT();
// auto elevLabels = elevLut->GetNumberOfTableValues();
auto elevLut = GetDivergingLUT();
auto elevLabels = 5;
if (surfaceName == "plane")
{
elevLabels = 1;
}
elevLut->SetTableRange(elevScalarRange.data());
cc->GetOutput()->GetPointData()->SetActiveScalars(curvature.c_str());
auto curvSR = cc->GetOutput()
->GetPointData()
->GetScalars(curvature.c_str())
->GetRange();
std::array<double, 2> curvScalarRange{curvSR[0], curvSR[1]};
// Generate the curvature.
auto precision = 10;
auto curvLut = GetCategoricalLUT();
curvLut->SetTableRange(curvScalarRange.data());
auto curvNumBands = curvLut->GetNumberOfTableValues();
auto curvBands =
GetBands(curvScalarRange.data(), curvNumBands, precision, false);
if (surfaceName == "random hills")
{
// These are my custom bands.
// Generated by first running:
// bands = GetBands(curvScalarRanges, curvNumBands, precision, false);
// then:
// std::vector<int> freq = Frequencies(bands, src);
// PrintBandsFrequencies(bands, freq);
// Finally using the output to create this table:
// std::vector<std::array<double, 2>> myBands = {
// {-0.630, -0.190}, {-0.190, -0.043}, {-0.043, -0.0136},
// {-0.0136, 0.0158}, {0.0158, 0.0452}, {0.0452, 0.0746},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<std::array<double, 2>> 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.
curvBands = GetCustomBands(curvScalarRange.data(), curvNumBands, myBands);
// Adjust the number of table values
curvLut->SetNumberOfTableValues(static_cast<vtkIdType>(curvBands.size()));
}
else if (surfaceName == "hills")
{
std::vector<std::array<double, 2>> 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.
curvBands = GetCustomBands(curvScalarRange.data(), curvNumBands, myBands);
// Adjust the number of table values
curvLut->SetNumberOfTableValues(static_cast<vtkIdType>(curvBands.size()));
}
// Let's do a frequency table.
auto curvFreq = GetFrequencies(curvBands, cc->GetOutput());
AdjustFrequencyRanges(curvBands, curvFreq);
if (frequencyTable)
{
PrintBandsFrequencies(curvBands, curvFreq);
}
// curvScalarRange[0] = curvBands.begin()->second[0];
// curvScalarRange[1] = std::prev(curvBands.end())->second[2];
curvLut->SetTableRange(curvScalarRange.data());
curvLut->SetNumberOfTableValues(curvBands.size());
AnnotateLut(curvLut, curvBands);
// Create the contour bands.
vtkNew<vtkBandedPolyDataContourFilter> bcf;
bcf->SetInputData(cc->GetOutput());
// We will use an indexed lookup table.
bcf->SetScalarModeToIndex();
bcf->GenerateContourEdgesOn();
// Use either the minimum or maximum value for each band.
int i = 0;
for (std::map<int, std::vector<double>>::const_iterator p = curvBands.begin();
p != curvBands.end(); ++p)
{
bcf->SetValue(i, p->second[2]);
++i;
}
// Generate the glyphs on the original surface.
std::array<double, 3> arrowScale{1, 1, 1};
if (surfaceName == "plane")
{
arrowScale = std::array<double, 3>{1.25, 1.25, 1.25};
}
auto scaleFactor = 1.0;
if (surfaceName == "hills")
{
scaleFactor = 0.5;
}
if (surfaceName == "sphere")
{
scaleFactor = 2.0;
}
auto glyph =
GetGlyphs(surfaceName, cc->GetOutput(), arrowScale, scaleFactor, false);
CurvatureData ret;
ret.bcf = bcf;
ret.glyph = glyph;
ret.curvScalarRange = curvScalarRange;
ret.elevScalarRange = elevScalarRange;
ret.curvLut = curvLut;
ret.elevLut = elevLut;
ret.elevLabels = elevLabels;
return ret;
}
CurvatureData
GenerateMeanCurvatures(std::string const& surfaceName, vtkPolyData* source,
std::vector<std::string> const& needsAdjusting,
bool frequencyTable)
{
std::string curvature = "Mean_Curvature";
vtkNew<vtkCurvatures> cc;
cc->SetInputData(source);
cc->SetCurvatureTypeToMean();
cc->Update();
if (std::find(needsAdjusting.begin(), needsAdjusting.end(), surfaceName) !=
needsAdjusting.end())
{
AdjustEdgeCurvatures(cc->GetOutput(), curvature);
}
if (surfaceName == "plane")
{
ConstrainCurvatures(cc->GetOutput(), curvature, 0.0, 0.0);
}
if (surfaceName == "sphere")
{
// Gaussian curvature is 1/r^2
auto radius = 10;
auto mean_curvature = 1.0 / radius;
ConstrainCurvatures(cc->GetOutput(), curvature, mean_curvature,
mean_curvature);
}
cc->GetOutput()->GetPointData()->SetActiveScalars("Elevation");
auto elevSR =
cc->GetOutput()->GetPointData()->GetScalars("Elevation")->GetRange();
std::array<double, 2> elevScalarRange{elevSR[0], elevSR[1]};
// auto elevLut = GetOrdinaLUT();
// auto elevLabels = elevLut->GetNumberOfTableValues();
auto elevLut = GetDivergingLUT();
auto elevLabels = 5;
if (surfaceName == "plane")
{
elevLabels = 1;
}
elevLut->SetTableRange(elevScalarRange.data());
cc->GetOutput()->GetPointData()->SetActiveScalars(curvature.c_str());
auto curvSR = cc->GetOutput()
->GetPointData()
->GetScalars(curvature.c_str())
->GetRange();
std::array<double, 2> curvScalarRange{curvSR[0], curvSR[1]};
// Generate the curvature.
auto precision = 10;
auto curvLut = GetCategoricalLUT();
curvLut->SetTableRange(curvScalarRange.data());
auto curvNumBands = curvLut->GetNumberOfTableValues();
auto curvBands =
GetBands(curvScalarRange.data(), curvNumBands, precision, false);
// If any bands need adjusting, we would do it here.
// Let's do a frequency table.
auto curvFreq = GetFrequencies(curvBands, cc->GetOutput());
AdjustFrequencyRanges(curvBands, curvFreq);
if (frequencyTable)
{
PrintBandsFrequencies(curvBands, curvFreq);
}
// curvScalarRange[0] = curvBands.begin()->second[0];
// curvScalarRange[1] = curvBands.rbegin()->second[2];
curvLut->SetTableRange(curvScalarRange.data());
curvLut->SetNumberOfTableValues(curvBands.size());
AnnotateLut(curvLut, curvBands);
// Create the contour bands.
vtkNew<vtkBandedPolyDataContourFilter> bcf;
bcf->SetInputData(cc->GetOutput());
// We will use an indexed lookup table.
bcf->SetScalarModeToIndex();
bcf->GenerateContourEdgesOn();
// Use either the minimum or maximum value for each band.
int i = 0;
for (std::map<int, std::vector<double>>::const_iterator p = curvBands.begin();
p != curvBands.end(); ++p)
{
bcf->SetValue(i, p->second[2]);
++i;
}
// Generate the glyphs on the original surface.
std::array<double, 3> arrowScale{1, 1, 1};
if (surfaceName == "plane")
{
arrowScale = std::array<double, 3>{1.25, 1.25, 1.25};
}
auto scaleFactor = 1.0;
if (surfaceName == "hills")
{
scaleFactor = 0.5;
}
if (surfaceName == "sphere")
{
scaleFactor = 2.0;
}
// Generate the glyphs on the original surface.
auto glyph =
GetGlyphs(surfaceName, cc->GetOutput(), arrowScale, scaleFactor, false);
CurvatureData ret;
ret.bcf = bcf;
ret.glyph = glyph;
ret.curvScalarRange = curvScalarRange;
ret.elevScalarRange = elevScalarRange;
ret.curvLut = curvLut;
ret.elevLut = elevLut;
ret.elevLabels = elevLabels;
return ret;
}
void AdjustEdgeCurvatures(vtkPolyData* source, std::string const& curvatureName,
double const& epsilon)
{
auto PointNeighbourhood =
[&source](vtkIdType const& pId) -> std::set<vtkIdType> {
// 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<vtkIdList> cellIds;
source->GetPointCells(pId, cellIds);
std::set<vtkIdType> neighbours;
for (vtkIdType i = 0; i < cellIds->GetNumberOfIds(); ++i)
{
auto cellId = cellIds->GetId(i);
vtkNew<vtkIdList> cellPointIds;
source->GetCellPoints(cellId, cellPointIds);
for (vtkIdType j = 0; j < cellPointIds->GetNumberOfIds(); ++j)
{
neighbours.insert(cellPointIds->GetId(j));
}
}
return neighbours;
};
auto ComputeDistance = [&source](vtkIdType const& ptIdA,
vtkIdType const& ptIdB) {
std::array<double, 3> ptA{0.0, 0.0, 0.0};
std::array<double, 3> ptB{0.0, 0.0, 0.0};
std::array<double, 3> 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<double>());
// Calculate the norm.
auto result = std::sqrt(std::inner_product(std::begin(ptC), std::end(ptC),
std::begin(ptC), 0.0));
return result;
};
source->GetPointData()->SetActiveScalars(curvatureName.c_str());
// Curvature as a vector.
auto array = source->GetPointData()->GetAbstractArray(curvatureName.c_str());
std::vector<double> curvatures;
for (vtkIdType i = 0; i < source->GetNumberOfPoints(); ++i)
{
curvatures.push_back(array->GetVariantValue(i).ToDouble());
}
// Get the boundary point IDs.
std::string name = "Ids";
#ifdef USE_USE_GENERATE_IDS
vtkNew<vtkGenerateIds> idFilter;
#else
vtkNew<vtkIdFilter> idFilter;
#endif
idFilter->SetInputData(source);
idFilter->SetPointIds(true);
idFilter->SetCellIds(false);
idFilter->SetPointIdsArrayName(name.c_str());
idFilter->SetCellIdsArrayName(name.c_str());
idFilter->Update();
vtkNew<vtkFeatureEdges> edges;
edges->SetInputConnection(idFilter->GetOutputPort());
edges->BoundaryEdgesOn();
edges->ManifoldEdgesOff();
edges->NonManifoldEdgesOff();
edges->FeatureEdgesOff();
edges->Update();
auto edgeAarray =
edges->GetOutput()->GetPointData()->GetAbstractArray(name.c_str());
std::vector<vtkIdType> boundaryIds;
for (vtkIdType i = 0; i < edges->GetOutput()->GetNumberOfPoints(); ++i)
{
boundaryIds.push_back(edgeAarray->GetVariantValue(i).ToInt());
}
// Remove duplicate Ids.
std::set<vtkIdType> pIdsSet(boundaryIds.begin(), boundaryIds.end());
for (auto const pId : boundaryIds)
{
auto pIdsNeighbors = PointNeighbourhood(pId);
std::set<vtkIdType> 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<double> curvs;
std::vector<double> dists;
for (auto const pIdN : pIdsNeighborsInterior)
{
curvs.push_back(curvatures[pIdN]);
dists.push_back(ComputeDistance(pIdN, pId));
}
std::vector<vtkIdType> nonZeroDistIds;
for (size_t i = 0; i < dists.size(); ++i)
{
if (dists[i] > 0)
{
nonZeroDistIds.push_back(i);
}
}
std::vector<double> curvsNonZero;
std::vector<double> 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<double> 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;
}
// 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;
}
}
}
if (static_cast<size_t>(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<vtkDoubleArray> adjustedCurvatures;
adjustedCurvatures->SetName(curvatureName.c_str());
for (auto curvature : curvatures)
{
adjustedCurvatures->InsertNextTuple1(curvature);
}
source->GetPointData()->AddArray(adjustedCurvatures);
source->GetPointData()->SetActiveScalars(curvatureName.c_str());
}
void ConstrainCurvatures(vtkPolyData* source, std::string const& curvatureName,
double const& lowerBound, double const& upperBound)
{
std::array<double, 2> bounds{0.0, 0.0};
if (lowerBound < upperBound)
{
bounds[0] = lowerBound;
bounds[1] = upperBound;
}
else
{
bounds[0] = upperBound;
bounds[1] = lowerBound;
}
source->GetPointData()->SetActiveScalars(curvatureName.c_str());
// Curvature as a vector.
auto array = source->GetPointData()->GetAbstractArray(curvatureName.c_str());
std::vector<double> 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<vtkDoubleArray> 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());
}
vtkSmartPointer<vtkPolyData> GetElevations(vtkPolyData* src)
{
double bounds[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
src->GetBounds(bounds);
if (std::abs(bounds[2]) < 1.0e-8 && std::abs(bounds[3]) < 1.0e-8)
{
bounds[3] = bounds[2] + 1;
}
vtkNew<vtkElevationFilter> elevFilter;
elevFilter->SetInputData(src);
elevFilter->SetLowPoint(0, bounds[2], 0);
elevFilter->SetHighPoint(0, bounds[3], 0);
elevFilter->SetScalarRange(bounds[2], bounds[3]);
elevFilter->Update();
return elevFilter->GetPolyDataOutput();
}
vtkSmartPointer<vtkPolyData> 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.
// We define the parameters for the hills here.
// There are four hills.
// [[0: x0, 1: y0, 2: x variance, 3: y variance, 4: amplitude]...]
std::vector<std::array<double, 5>> hd{{-2.5, -2.5, 2.5, 6.5, 3.5},
{2.5, 2.5, 2.5, 2.5, 2},
{5.0, -2.5, 1.5, 1.5, 2.5},
{-5.0, 5, 2.5, 3.0, 3}};
// Three hills.
// std::vector<std::array<double, 5>> 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<std::array<double, 5>> hd{{5.0, -2.5, 1.5, 1.5, 2.5},
// {-5.0, 5, 2.5, 3.0, 3}};
std::array<double, 2> xx{0.0, 0.0};
vtkNew<vtkPoints> points;
points->SetDataTypeToDouble();
auto y = yMin;
for (auto r = 0; r < xRes; ++r)
{
auto x = xMin;
for (auto c = 0; c < yRes; ++c)
{
auto z = 0.0;
for (size_t j = 0; j < hd.size(); ++j)
{
xx[0] = std::pow(x - hd[j][0] / hd[j][2], 2.0);
xx[1] = std::pow(y - hd[j][1] / hd[j][3], 2.0);
z += hd[j][4] * std::exp(-(xx[0] + xx[1]) / 2.0);
}
x += dx;
points->InsertNextPoint(x, y, z);
}
y += dy;
}
// Triangulate each quad.
vtkNew<vtkCellArray> triangles;
auto tNum = 0;
for (auto r = 0; r < xRes - 1; ++r)
{
for (auto c = 0; c < yRes - 1; ++c)
{
std::vector<int> indices;
// We select the index of each point so that
// the ordering is counterclockwise.
auto rc = r * yRes;
// Triangle 1.
indices.push_back(rc + c);
indices.push_back(indices[0] + 1);
indices.push_back(rc + yRes + c);
vtkNew<vtkTriangle> triangle1;
for (auto i = 0; i < indices.size(); ++i)
{
triangle1->GetPointIds()->SetId(i, indices[i]);
}
triangles->InsertNextCell(triangle1);
indices.clear();
tNum += 1;
// Triangle 2.
indices.push_back(rc + c + 1);
indices.push_back(rc + yRes + c + 1);
indices.push_back(rc + yRes + c);
vtkNew<vtkTriangle> triangle2;
for (auto i = 0; i < indices.size(); ++i)
{
triangle2->GetPointIds()->SetId(i, indices[i]);
}
triangles->InsertNextCell(triangle2);
indices.clear();
tNum += 1;
}
}
vtkNew<vtkPolyData> polydata;
polydata->SetPoints(points);
polydata->SetPolys(triangles);
vtkNew<vtkDoubleArray> elevation;
elevation->SetNumberOfTuples(points->GetNumberOfPoints());
vtkNew<vtkFloatArray> textures;
textures->SetNumberOfComponents(2);
textures->SetNumberOfTuples(2 * polydata->GetNumberOfPoints());
textures->SetName("Textures");
for (auto i = 0; i < xRes; ++i)
{
float tc[2]{0.0, 0.0};
tc[0] = i / (xRes - 1.0);
for (auto j = 0; j < yRes; ++j)
{
// tc[1] = 1.0 - j / (yRes - 1.0);
tc[1] = j / (yRes - 1.0);
textures->SetTuple(static_cast<vtkIdType>(i) * yRes + j, tc);
}
}
polydata->GetPointData()->SetTCoords(textures);
auto bounds = polydata->GetBounds();
vtkNew<vtkElevationFilter> elevationFilter;
elevationFilter->SetLowPoint(0.0, 0.0, bounds[4]);
elevationFilter->SetHighPoint(0.0, 0.0, bounds[5]);
elevationFilter->SetInputData(polydata);
vtkNew<vtkPolyDataNormals> normals;
normals->SetInputConnection(elevationFilter->GetOutputPort());
normals->SetFeatureAngle(30);
normals->SplittingOff();
vtkNew<vtkTransform> tr;
tr->RotateX(-90);
vtkNew<vtkTransformPolyDataFilter> tf;
tf->SetInputConnection(normals->GetOutputPort());
tf->SetTransform(tr);
tf->Update();
return tf->GetOutput();
}
vtkSmartPointer<vtkPolyData> GetParametricHills()
{
vtkNew<vtkParametricRandomHills> fn;
fn->AllowRandomGenerationOn();
fn->SetRandomSeed(1);
fn->SetNumberOfHills(30);
vtkNew<vtkParametricFunctionSource> 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<vtkTransform> transform;
transform->Translate(0.0, 5.0, 15.0);
transform->RotateX(-90.0);
vtkNew<vtkTransformPolyDataFilter> transformFilter;
transformFilter->SetInputConnection(source->GetOutputPort());
transformFilter->SetTransform(transform);
transformFilter->Update();
return transformFilter->GetOutput();
}
vtkSmartPointer<vtkPolyData> GetParametricTorus()
{
vtkNew<vtkParametricTorus> fn;
fn->SetRingRadius(5);
fn->SetCrossSectionRadius(2);
vtkNew<vtkParametricFunctionSource> 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<vtkTransform> transform;
transform->RotateX(-90.0);
vtkNew<vtkTransformPolyDataFilter> transformFilter;
transformFilter->SetInputConnection(source->GetOutputPort());
transformFilter->SetTransform(transform);
transformFilter->Update();
return transformFilter->GetOutput();
}
vtkSmartPointer<vtkPolyData> GetPlane()
{
vtkNew<vtkPlaneSource> source;
source->SetOrigin(-10.0, -10.0, 0.0);
source->SetPoint1(10.0, -10.0, 0.0);
source->SetPoint2(-10.0, 10.0, 0.0);
source->SetXResolution(5);
source->SetYResolution(5);
source->Update();
vtkNew<vtkTransform> transform;
transform->Translate(0.0, 0.0, 0.0);
transform->RotateX(-90.0);
vtkNew<vtkTransformPolyDataFilter> 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<vtkTriangleFilter> tri;
tri->SetInputConnection(transformFilter->GetOutputPort());
// Pass it though a CleanPolyDataFilter and merge any points which
// are coincident, or very close
vtkNew<vtkCleanPolyData> cleaner;
cleaner->SetInputConnection(tri->GetOutputPort());
cleaner->SetTolerance(0.005);
cleaner->Update();
return cleaner->GetOutput();
}
vtkSmartPointer<vtkPolyData> GetSphere()
{
vtkNew<vtkSphereSource> source;
source->SetCenter(0.0, 0.0, 0.0);
source->SetRadius(10.0);
source->SetThetaResolution(32);
source->SetPhiResolution(32);
source->Update();
return source->GetOutput();
}
vtkSmartPointer<vtkPolyData> GetTorus()
{
vtkNew<vtkSuperquadricSource> source;
source->SetCenter(0.0, 0.0, 0.0);
source->SetScale(1.0, 1.0, 1.0);
source->SetPhiResolution(64);
source->SetThetaResolution(64);
source->SetThetaRoundness(1);
source->SetThickness(0.5);
source->SetSize(10);
source->SetToroidal(1);
source->Update();
auto bounds = source->GetOutput()->GetBounds();
// The quadric is made of strips, so pass it through a triangle filter as
// the curvature filter only operates on polys
vtkNew<vtkTriangleFilter> 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<vtkCleanPolyData> cleaner;
cleaner->SetInputConnection(tri->GetOutputPort());
cleaner->SetTolerance(0.005);
cleaner->Update();
return cleaner->GetOutput();
}
vtkSmartPointer<vtkPolyData> GetSource(std::string const& source)
{
std::string surface = source;
std::transform(surface.begin(), surface.end(), surface.begin(),
[](unsigned char c) { return std::tolower(c); });
std::map<std::string, int> available_surfaces = {
{"hills", 0}, {"parametric torus", 1}, {"plane", 2},
{"random hills", 3}, {"sphere", 4}, {"torus", 5}};
if (available_surfaces.find(surface) == available_surfaces.end())
{
std::cout << "The surface is not available." << std::endl;
std::cout << "Using Random Hills instead." << std::endl;
surface = "random hills";
}
switch (available_surfaces[surface])
{
case 0:
return GetHills();
break;
case 1:
return GetParametricTorus();
break;
case 2:
return GetElevations(GetPlane());
break;
case 3:
return GetParametricHills();
break;
case 4:
return GetElevations(GetSphere());
break;
case 5:
return GetElevations(GetTorus());
break;
}
return GetParametricHills();
}
// vtkSmartPointer<vtkColorSeries> GetColorSeries()
// {
// vtkNew<vtkColorSeries> colorSeries;
// // Select a color scheme.
// int colorSeriesEnum;
// // colorSeriesEnum = colorSeries->BREWER_DIVERGING_BROWN_BLUE_GREEN_9;
// // colorSeriesEnum = colorSeries->BREWER_DIVERGING_SPECTRAL_10;
// // colorSeriesEnum = colorSeries->BREWER_DIVERGING_SPECTRAL_3;
// // colorSeriesEnum = colorSeries->BREWER_DIVERGING_PURPLE_ORANGE_9;
// // colorSeriesEnum = colorSeries->BREWER_SEQUENTIAL_BLUE_PURPLE_9;
// // colorSeriesEnum = colorSeries->BREWER_SEQUENTIAL_BLUE_GREEN_9;
// colorSeriesEnum = colorSeries->BREWER_QUALITATIVE_SET3;
// // colorSeriesEnum = colorSeries->CITRUS;
// colorSeries->SetColorScheme(colorSeriesEnum);
// return colorSeries;
// }
vtkNew<vtkLookupTable> GetCategoricalLUT()
{
vtkNew<vtkColorSeries> colorSeries;
colorSeries->SetColorScheme(colorSeries->BREWER_QUALITATIVE_SET3);
// Make the lookup table.
vtkNew<vtkLookupTable> lut;
colorSeries->BuildLookupTable(lut, vtkColorSeries::CATEGORICAL);
lut->SetNanColor(0, 0, 0, 1);
return lut;
}
vtkNew<vtkLookupTable> GetOrdinaLUT()
{
vtkNew<vtkColorSeries> colorSeries;
colorSeries->SetColorScheme(
colorSeries->BREWER_DIVERGING_BROWN_BLUE_GREEN_11);
// Make the lookup table.
vtkNew<vtkLookupTable> lut;
colorSeries->BuildLookupTable(lut, vtkColorSeries::ORDINAL);
lut->SetNanColor(0, 0, 0, 1);
return lut;
}
// clang-format off
/**
* 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
*
*/
// clang-format on
vtkNew<vtkLookupTable> GetDivergingLUT()
{
vtkNew<vtkColorTransferFunction> ctf;
ctf->SetColorSpaceToDiverging();
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);
auto tableSize = 256;
vtkNew<vtkLookupTable> lut;
lut->SetNumberOfTableValues(tableSize);
lut->Build();
for (auto i = 0; i < lut->GetNumberOfColors(); ++i)
{
std::array<double, 3> rgb{0.0, 0.0, 0.0};
ctf->GetColor(static_cast<double>(i) / lut->GetNumberOfColors(),
rgb.data());
std::array<double, 4> rgba{0.0, 0.0, 0.0, 1.0};
std::copy(std::begin(rgb), std::end(rgb), std::begin(rgba));
lut->SetTableValue(static_cast<vtkIdType>(i), rgba.data());
}
return lut;
}
vtkNew<vtkLookupTable> ReverseLUT(vtkLookupTable* lut)
{
// First do a deep copy just to get the whole structure
// and then reverse the colors and annotations.
vtkNew<vtkLookupTable> lutr;
lutr->DeepCopy(lut);
vtkIdType t = lut->GetNumberOfTableValues() - 1;
for (vtkIdType i = t; i >= 0; --i)
{
std::array<double, 3> rgb{0.0, 0.0, 0.0};
std::array<double, 4> rgba{0.0, 0.0, 0.0, 1.0};
lut->GetColor(i, rgb.data());
std::copy(std::begin(rgb), std::end(rgb), std::begin(rgba));
rgba[3] = lut->GetOpacity(i);
lutr->SetTableValue(t - i, rgba.data());
}
t = lut->GetNumberOfAnnotatedValues() - 1;
for (vtkIdType i = t; i >= 0; --i)
{
lutr->SetAnnotation(t - i, lut->GetAnnotation(i));
}
return lutr;
}
vtkNew<vtkGlyph3D> GetGlyphs(std::string const& name, vtkPolyData* src,
std::array<double, 3> const& arrowScale,
double const& scaleFactor,
bool const& reverseNormals)
{
// Choose a random subset of points.
vtkNew<vtkMaskPoints> maskPts;
maskPts->RandomModeOn();
if (name == "plane")
{
maskPts->SetOnRatio(1);
}
else
{
maskPts->SetOnRatio(5);
}
// Sometimes the contouring algorithm can create a volume whose gradient
// vector and ordering of polygon(using the right hand rule) are
// inconsistent. vtkReverseSense cures this problem.
vtkNew<vtkReverseSense> reverse;
if (reverseNormals)
{
reverse->SetInputData(src);
reverse->ReverseCellsOn();
reverse->ReverseNormalsOn();
maskPts->SetInputConnection(reverse->GetOutputPort());
}
else
{
maskPts->SetInputData(src);
}
// Source for the glyph filter.
vtkNew<vtkArrowSource> arrow;
arrow->SetShaftResolution(16);
arrow->SetShaftRadius(0.03);
arrow->SetTipResolution(16);
arrow->SetTipLength(0.3);
arrow->SetTipRadius(0.1);
vtkNew<vtkTransform> tr;
tr->Scale(arrowScale.data());
vtkNew<vtkTransformPolyDataFilter> tf;
tf->SetTransform(tr);
tf->SetInputConnection(arrow->GetOutputPort());
tf->Update();
vtkNew<vtkGlyph3D> glyph;
glyph->SetSourceConnection(tf->GetOutputPort());
glyph->SetInputConnection(maskPts->GetOutputPort());
glyph->SetVectorModeToUseNormal();
glyph->SetScaleFactor(scaleFactor);
glyph->SetColorModeToColorByVector();
glyph->SetScaleModeToScaleByVector();
glyph->OrientOn();
glyph->Update();
return glyph;
}
std::map<int, std::vector<double>> GetBands(double const dR[2],
int const& numberOfBands,
int const& precision,
bool const& nearestInteger)
{
auto prec = abs(precision);
prec = (prec > 14) ? 14 : prec;
auto RoundOff = [&prec](const double& x) {
auto pow_10 = std::pow(10.0, prec);
return std::round(x * pow_10) / pow_10;
};
std::map<int, std::vector<double>> bands;
if ((dR[1] < dR[0]) || (numberOfBands <= 0))
{
return bands;
}
double x[2]{0.0, 0.0};
for (int i = 0; i < 2; ++i)
{
x[i] = dR[i];
}
if (nearestInteger)
{
x[0] = std::floor(x[0]);
x[1] = std::ceil(x[1]);
}
double dx = (x[1] - x[0]) / static_cast<double>(numberOfBands);
std::vector<double> b;
b.push_back(x[0]);
b.push_back(x[0] + dx / 2.0);
b.push_back(x[0] + dx);
for (int i = 0; i < numberOfBands; ++i)
{
if (i == 0)
{
for (std::vector<double>::iterator p = b.begin(); p != b.end(); ++p)
{
*p = RoundOff(*p);
}
b[0] = x[0];
}
bands[i] = b;
for (std::vector<double>::iterator p = b.begin(); p != b.end(); ++p)
{
*p = RoundOff(*p + dx);
}
}
return bands;
}
std::map<int, std::vector<double>>
GetCustomBands(double const dR[2], int const& numberOfBands,
std::vector<std::array<double, 2>> const& myBands)
{
std::map<int, std::vector<double>> bands;
if ((dR[1] < dR[0]) || (numberOfBands <= 0))
{
return bands;
}
std::vector<std::array<double, 2>> x;
std::copy(myBands.begin(), myBands.end(), std::back_inserter(x));
// Determine the index of the range minimum and range maximum.
int idxMin = 0;
for (auto idx = 0; idx < static_cast<int>(myBands.size()); ++idx)
{
if (dR[0] < myBands[idx][1] && dR[0] >= myBands[idx][0])
{
idxMin = idx;
break;
}
}
int idxMax = static_cast<int>(myBands.size()) - 1;
for (int idx = static_cast<int>(myBands.size()) - 1; idx >= 0; --idx)
{
if (dR[1] < myBands[idx][1] && dR[1] >= myBands[idx][0])
{
idxMax = static_cast<int>(idx);
break;
}
}
// Set the minimum to match the range minimum.
x[idxMin][0] = dR[0];
x[idxMax][1] = dR[1];
for (int i = idxMin; i < idxMax + 1; ++i)
{
std::vector<double> b(3);
b[0] = x[i][0];
b[1] = x[i][0] + (x[i][1] - x[i][0]) / 2.0;
b[2] = x[i][1];
bands[i] = b;
}
return bands;
}
std::map<int, std::vector<double>> GetIntegralBands(double const dR[2])
{
if (dR[1] < dR[0])
{
std::map<int, std::vector<double>> bands;
return bands;
}
double x[2]{0.0, 0.0};
for (int i = 0; i < 2; ++i)
{
x[i] = dR[i];
}
x[0] = std::floor(x[0]);
x[1] = std::ceil(x[1]);
int numberOfBands = static_cast<int>(std::abs(x[1]) + std::abs(x[0]));
return GetBands(x, numberOfBands, false);
}
std::map<int, int> GetFrequencies(std::map<int, std::vector<double>>& bands,
vtkPolyData* src)
{
std::map<int, int> freq;
for (auto i = 0; i < static_cast<int>(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<int>(bands.size()); ++j)
{
if (*x <= bands[j][2])
{
freq[j] = freq[j] + 1;
break;
}
}
}
return freq;
}
void AdjustFrequencyRanges(std::map<int, std::vector<double>>& bands,
std::map<int, int>& freq)
{
// Get the indices of the first and last non-zero elements.
auto first = 0;
for (auto i = 0; i < static_cast<int>(freq.size()); ++i)
{
if (freq[i] != 0)
{
first = i;
break;
}
}
std::vector<int> keys;
for (std::map<int, int>::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<int, int>::iterator freqItr;
freqItr = freq.find(first);
freq.erase(freq.begin(), freqItr);
freqItr = ++freq.find(last);
freq.erase(freqItr, freq.end());
std::map<int, std::vector<double>>::iterator bandItr;
bandItr = bands.find(first);
bands.erase(bands.begin(), bandItr);
bandItr = ++bands.find(last);
bands.erase(bandItr, bands.end());
// Reindex freq and bands.
std::map<int, int> adjFreq;
int idx = 0;
for (auto p : freq)
{
adjFreq[idx] = p.second;
++idx;
}
std::map<int, std::vector<double>> adjBands;
idx = 0;
for (auto const& p : bands)
{
adjBands[idx] = p.second;
++idx;
}
bands = adjBands;
freq = adjFreq;
}
void PrintBandsFrequencies(std::map<int, std::vector<double>> const& bands,
std::map<int, int>& freq, int const& precision)
{
auto prec = abs(precision);
prec = (prec > 14) ? 14 : prec;
if (bands.size() != freq.size())
{
std::cout << "Bands and frequencies must be the same size." << std::endl;
return;
}
std::ostringstream os;
os << "Bands & Frequencies:\n";
size_t idx = 0;
auto total = 0;
auto width = prec + 6;
for (std::map<int, std::vector<double>>::const_iterator p = bands.begin();
p != bands.end(); ++p)
{
total += freq[p->first];
for (std::vector<double>::const_iterator q = p->second.begin();
q != p->second.end(); ++q)
{
if (q == p->second.begin())
{
os << std::setw(4) << idx << " [";
}
if (q == std::prev(p->second.end()))
{
os << std::fixed << std::setw(width) << std::setprecision(prec) << *q
<< "]: " << std::setw(8) << freq[p->first] << "\n";
}
else
{
os << std::fixed << std::setw(width) << std::setprecision(prec) << *q
<< ", ";
}
}
++idx;
}
width = 3 * width + 13;
os << std::left << std::setw(width) << "Total" << std::right << std::setw(8)
<< total << std::endl;
std::cout << os.str() << endl;
}
void AnnotateLut(vtkLookupTable* lut, std::map<int, std::vector<double>> bands)
{
// We will use the midpoint of the band as the label.
std::vector<std::string> labels;
for (std::map<int, std::vector<double>>::const_iterator p = bands.begin();
p != bands.end(); ++p)
{
labels.push_back(fmt::format("{:4.2f}", p->second[1]));
}
// Annotate
vtkNew<vtkVariantArray> values;
for (size_t i = 0; i < labels.size(); ++i)
{
values->InsertNextValue(vtkVariant(labels[i]));
}
lut->ResetAnnotations();
for (vtkIdType i = 0; i < values->GetNumberOfTuples(); ++i)
{
lut->SetAnnotation(i, values->GetValue(i).ToString());
}
}
vtkNew<vtkScalarBarWidget>
MakeScalarBarWidget(ScalarBarProperties& sbProperties,
vtkTextProperty* textProperty,
vtkTextProperty* labelTextProperty, vtkRenderer* ren,
vtkRenderWindowInteractor* iren)
{
vtkNew<vtkScalarBarActor> sbActor;
sbActor->SetLookupTable(sbProperties.lut);
sbActor->SetTitle(sbProperties.titleText.c_str());
sbActor->UnconstrainedFontSizeOn();
sbActor->SetNumberOfLabels(sbProperties.number_of_labels);
sbActor->SetTitleTextProperty(textProperty);
sbActor->SetTitleTextProperty(textProperty);
sbActor->SetLabelTextProperty(labelTextProperty);
sbActor->SetLabelFormat(sbProperties.labelFormat.c_str());
vtkNew<vtkScalarBarRepresentation> sbRep;
sbRep->EnforceNormalizedViewportBoundsOn();
sbRep->SetOrientation(sbProperties.orientation);
// Set the position.
sbRep->GetPositionCoordinate()->SetCoordinateSystemToNormalizedViewport();
sbRep->GetPosition2Coordinate()->SetCoordinateSystemToNormalizedViewport();
if (sbProperties.orientation)
{
auto p1 = sbProperties.positionV["p"];
auto p2 = sbProperties.positionV["p2"];
sbRep->GetPositionCoordinate()->SetValue(p1.data());
sbRep->GetPosition2Coordinate()->SetValue(p2.data());
}
else
{
auto p1 = sbProperties.positionH["p"];
auto p2 = sbProperties.positionH["p2"];
sbRep->GetPositionCoordinate()->SetValue(p1.data());
sbRep->GetPosition2Coordinate()->SetValue(p2.data());
}
vtkNew<vtkScalarBarWidget> widget;
widget->SetRepresentation(sbRep);
widget->SetScalarBarActor(sbActor);
widget->SetDefaultRenderer(ren);
widget->SetInteractor(iren);
widget->EnabledOn();
return widget;
}
TTextPositions GetTextPositions(std::vector<std::string> const& names,
int const justification,
int const vertical_justification,
double const width, double const height)
{
// The gap between the left or right edge of the screen and the text.
auto dx = 0.02;
auto w = abs(width);
if (w > 0.96)
{
w = 0.96;
}
auto y0 = 0.01;
auto h = abs(height);
if (h > 0.9)
{
h = 0.9;
}
auto dy = h;
if (vertical_justification == VTK_TEXT_TOP)
{
y0 = 1.0 - (dy + y0);
}
if (vertical_justification == VTK_TEXT_CENTERED)
{
y0 = 0.5 - (dy / 2.0 + y0);
}
auto minmaxIt =
std::minmax_element(names.begin(), names.end(),
[](const std::string& a, const std::string& b) {
return a.length() < b.length();
});
// auto nameLenMin = minmaxIt.first->size();
auto nameLenMax = minmaxIt.second->size();
TTextPositions textPositions;
for (const auto& k : names)
{
auto sz = k.size();
auto delta_sz = w * sz / nameLenMax;
if (delta_sz > w)
{
delta_sz = w;
}
double x0 = 0;
if (justification == VTK_TEXT_CENTERED)
{
x0 = 0.5 - delta_sz / 2.0;
}
else if (justification == VTK_TEXT_RIGHT)
{
x0 = 1.0 - dx - delta_sz;
}
else
{
// Default is left justification.
x0 = dx;
}
textPositions[k] = {{"p", {x0, y0}}, {"p2", {delta_sz, dy}}};
// For testing.
// std::cout << k << std::endl;
// std::cout << " p: " << textPositions[k]["p"][0] << ", "
// << textPositions[k]["p"][1] << std::endl;
// std::cout << " p2: " << textPositions[k]["p2"][0] << ", "
// << textPositions[k]["p2"][1] << std::endl;
}
return textPositions;
}
std::string Title(const std::string& input)
{
std::string result = input;
bool capitalizeNext = true;
for (char& c : result)
{
if (std::isspace(static_cast<unsigned char>(c)))
{
capitalizeNext = true;
}
else if (capitalizeNext)
{
c = static_cast<char>(std::toupper(static_cast<unsigned char>(c)));
capitalizeNext = false;
}
else
{
c = static_cast<char>(std::tolower(static_cast<unsigned char>(c)));
}
}
return result;
}
} // namespace
CMakeLists.txt¶
cmake_minimum_required(VERSION 3.12 FATAL_ERROR)
project(CurvatureBandsWithGlyphs)
find_package(VTK COMPONENTS
CommonColor
CommonComputationalGeometry
CommonCore
CommonDataModel
CommonTransforms
FiltersCore
FiltersGeneral
FiltersModeling
FiltersSources
InteractionStyle
InteractionWidgets
RenderingAnnotation
RenderingContextOpenGL2
RenderingCore
RenderingFreeType
RenderingGL2PSOpenGL2
RenderingOpenGL2
)
if (NOT VTK_FOUND)
message(FATAL_ERROR "CurvatureBandsWithGlyphs: Unable to find the VTK build folder.")
endif()
# Prevent a "command line is too long" failure in Windows.
set(CMAKE_NINJA_FORCE_RESPONSE_FILE "ON" CACHE BOOL "Force Ninja to use response files.")
add_executable(CurvatureBandsWithGlyphs MACOSX_BUNDLE CurvatureBandsWithGlyphs.cxx )
target_link_libraries(CurvatureBandsWithGlyphs PRIVATE ${VTK_LIBRARIES}
)
# vtk_module_autoinit is needed
vtk_module_autoinit(
TARGETS CurvatureBandsWithGlyphs
MODULES ${VTK_LIBRARIES}
)
Download and Build CurvatureBandsWithGlyphs¶
Click here to download CurvatureBandsWithGlyphs and its CMakeLists.txt file. Once the tarball CurvatureBandsWithGlyphs.tar has been downloaded and extracted,
cd CurvatureBandsWithGlyphs/build
If VTK is installed:
cmake ..
If VTK is not installed but compiled on your system, you will need to specify the path to your VTK build:
cmake -DVTK_DIR:PATH=/home/me/vtk_build ..
Build the project:
make
and run it:
./CurvatureBandsWithGlyphs
WINDOWS USERS
Be sure to add the VTK bin directory to your path. This will resolve the VTK dll's at run time.