CurvaturesAdjustEdges

#include <vtkActor.h>
#include <vtkActor2D.h>
#include <vtkCamera.h>
#include <vtkColorTransferFunction.h>
#include <vtkCubeSource.h>
#include <vtkCurvatures.h>
#include <vtkDelaunay2D.h>
#include <vtkDiscretizableColorTransferFunction.h>
#include <vtkDoubleArray.h>
#include <vtkElevationFilter.h>
#include <vtkFeatureEdges.h>
#include <vtkFloatArray.h>
#include <vtkIdList.h>
#include <vtkInteractorStyleTrackballCamera.h>
#include <vtkLinearSubdivisionFilter.h>
#include <vtkLookupTable.h>
#include <vtkNamedColors.h>
#include <vtkNew.h>
#include <vtkParametricBour.h>
#include <vtkParametricEnneper.h>
#include <vtkParametricFunctionSource.h>
#include <vtkParametricMobius.h>
#include <vtkParametricRandomHills.h>
#include <vtkParametricTorus.h>
#include <vtkPointData.h>
#include <vtkPolyData.h>
#include <vtkPolyDataMapper.h>
#include <vtkPolyDataNormals.h>
#include <vtkPolyDataTangents.h>
#include <vtkRenderWindow.h>
#include <vtkRenderWindowInteractor.h>
#include <vtkRenderer.h>
#include <vtkScalarBarActor.h>
#include <vtkScalarBarRepresentation.h>
#include <vtkScalarBarWidget.h>
#include <vtkSmartPointer.h>
#include <vtkTextActor.h>
#include <vtkTextProperty.h>
#include <vtkTextRepresentation.h>
#include <vtkTextWidget.h>
#include <vtkTexturedSphereSource.h>
#include <vtkTransform.h>
#include <vtkTransformPolyDataFilter.h>
#include <vtkTriangle.h>
#include <vtkTriangleFilter.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 <filesystem>
// #include <functional>
// #include <iomanip>
// #include <iostream>
// #include <iterator>
// #include <map>
// #include <numeric>
// #include <set>
// #include <sstream>

namespace fs = std::filesystem;

namespace {
struct CurvatureData
{
  // Holds the necessary data to build the view.
  std::string curvatureName;
  std::array<double, 2> scalarRange;
  int scalarBarLabels = 5;
  vtkSmartPointer<vtkLookupTable> lut;
};

/** 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);

// Some sample surfaces to try.
vtkSmartPointer<vtkPolyData> GetBour();
vtkSmartPointer<vtkPolyData> GetCube();
vtkSmartPointer<vtkPolyData> GetEnneper();
vtkSmartPointer<vtkPolyData> GetHills();
vtkSmartPointer<vtkPolyData> GetMobius();
vtkSmartPointer<vtkPolyData> GetRandomHills();
vtkSmartPointer<vtkPolyData> GetSphere();
vtkSmartPointer<vtkPolyData> GetTorus();

vtkSmartPointer<vtkPolyData> GetSource(std::string const& soource);

vtkSmartPointer<vtkLookupTable> GetDivergingLut();
vtkSmartPointer<vtkLookupTable> GetDivergingLut1();
vtkSmartPointer<vtkDiscretizableColorTransferFunction> GetCTF();

std::map<int, std::vector<double>> GetBands(double const dR[2],
                                            int const& numberOfBands,
                                            int const& precision = 2,
                                            bool const& nearestInteger = false);
//! 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);
//!
/*
 * 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::string const& curvature,
                           std::map<int, std::vector<double>> const& bands,
                           std::map<int, int>& freq);

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;
  // vtkSmartPointer<vtkDiscretizableColorTransferFunction> 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.05}}, {"p2", {0.1, 0.7}}};
  TTextPosition defaultH = {{"p", {0.125, 0.05}}, {"p2", {0.75, 0.1}}};
  // Modify these as needed.
  TTextPosition positionV = {{"p", {0.85, 0.05}}, {"p2", {0.1, 0.7}}};
  TTextPosition positionH = {{"p", {0.125, 0.05}}, {"p2", {0.75, 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{"Display the Gaussian and Mean curvatures of a"
               "surface with adjustments for edge effects."};

  // 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 frequencyTable{false};
  app.add_flag("-f, --frequencyTable", frequencyTable,
               "Display the frequency table.");

  CLI11_PARSE(app, argc, argv);

  if (surfaceName.empty())
  {
    surfaceName = "Random Hills";
  }
  // 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";
  }

  // Get the surface
  std::vector<std::string> availableSurfaces{
      "bour",   "cube",         "hills",  "enneper",
      "mobius", "random hills", "sphere", "torus",
  };
  // Surfaces whose curvatures need to be adjusted along the edges of the
  // surface or constrained.
  std::vector<std::string> needsAdjusting{"bour", "enneper", "hills",
                                          "random hills", "torus"};

  auto it = std::find(availableSurfaces.begin(), availableSurfaces.end(),
                      surfaceName);
  if (it == availableSurfaces.end())
  {
    std::string res = fmt::format(
        "Nonexistent surface {}\nAvailable surfaces are:\n   ", surfaceName);
    for (const auto& element : availableSurfaces)
    {
      res += fmt::format("{:s}, ", element);
    }
    if (res.length() >= 2)
    {
      auto pos = res.length() - 2;
      res.replace(pos, 2, "");
    }
    res += "\nIf a name has spaces in it, delineate the name with quotes e.g. "
           "\"random hills \"";
    std::cout << res << std::endl;

    return EXIT_FAILURE;
  }

  auto source = GetSource(surfaceName);

  std::map<std::string, CurvatureData> curvatures;

  curvatures["Gauss_Curvature"] = GenerateGaussianCurvatures(
      surfaceName, source, needsAdjusting, frequencyTable);
  curvatures["Mean_Curvature"] = GenerateMeanCurvatures(
      surfaceName, source, needsAdjusting, frequencyTable);

  // Let's visualise what we have done.

  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());

  auto lut = GetDivergingLut();
  // auto lut = GetDivergingLut1();
  // auto lut = GetCTF();

  // Define viewport ranges
  std::array<double, 2> xmins{0, 0.5};
  std::array<double, 2> xmaxs{0.5, 1};
  std::array<double, 2> ymins{0, 0};
  std::array<double, 2> ymaxs{1.0, 1.0};

  auto windowHeight = 512;
  auto windowWidth = 2 * windowHeight;

  // Create the RenderWindow, Renderers and Interactor.
  vtkNew<vtkRenderWindow> renWin;
  renWin->SetSize(windowWidth, windowHeight);
  vtkNew<vtkRenderWindowInteractor> iRen;
  iRen->SetRenderWindow(renWin);
  vtkNew<vtkInteractorStyleTrackballCamera> style;
  iRen->SetInteractorStyle(style);
  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();

  // These surfaces have just one curvature value.
  std::vector<std::string> singleGaussCurv = {"bour", "enneper", "cube",
                                              "mobius", "sphere"};
  std::vector<std::string> singleMeanCurv = {"bour", "enneper", "cube",
                                             "sphere"};
  auto first = true;
  std::map<size_t, vtkSmartPointer<vtkScalarBarWidget>> sbWidgets;
  std::map<size_t, vtkSmartPointer<vtkRenderer>> renderers;

  for (const auto& kv : curvatures)
  {
    std::string curvatureTitle = kv.first;
    std::replace(curvatureTitle.begin(), curvatureTitle.end(), '_', '\n');

    CurvatureData curvatureData = kv.second;
    source->GetPointData()->SetActiveScalars(
        curvatureData.curvatureName.c_str());

    vtkNew<vtkPolyDataMapper> mapper;
    mapper->SetInputData(source);
    mapper->SetScalarModeToUsePointFieldData();
    mapper->SelectColorArray(curvatureData.curvatureName.c_str());
    mapper->SetScalarRange(curvatureData.scalarRange.data());
    mapper->SetLookupTable(curvatureData.lut);

    vtkNew<vtkActor> actor;
    actor->SetMapper(mapper);

    vtkNew<vtkRenderer> renderer;
    renderer->SetBackground(
        colors->GetColor3d("ParaViewWarmGrayBkg").GetData());
    if (first)
    {
      textWidget->SetDefaultRenderer(renderer);
      first = false;
    }

    auto sbProperties = ScalarBarProperties();
    sbProperties.titleText = curvatureData.curvatureName + "\n";
    std::replace(sbProperties.titleText.begin(), sbProperties.titleText.end(),
                 '_', ' ');
    sbProperties.number_of_labels =
        std::min<int>(5, curvatureData.scalarBarLabels);
    sbProperties.lut = curvatureData.lut;
    sbProperties.orientation = false;
    if (curvatureData.curvatureName == "Gauss_Curvature")
    {
      auto it = std::find(singleGaussCurv.begin(), singleGaussCurv.end(),
                          surfaceName);
      if (it != singleGaussCurv.end())
      {
        sbProperties.positionH = {{"p", {0.3, 0.05}}, {"p2", {0.35, 0.1}}};
      }
    }
    if (curvatureData.curvatureName == "Mean_Curvature")
    {
      auto it =
          std::find(singleMeanCurv.begin(), singleMeanCurv.end(), surfaceName);
      if (it != singleMeanCurv.end())
      {
        sbProperties.positionH = {{"p", {0.3, 0.05}}, {"p2", {0.35, 0.1}}};
      }
    }

    size_t idx = 0;
    if (curvatureData.curvatureName == "Gauss_Curvature")
    {
      idx = 0;
    }
    else
    {
      idx = 1;
    }
    sbWidgets[idx] =
        MakeScalarBarWidget(sbProperties, textProp, textProp, renderer, iRen);
    // sbProperties.positionH = sbProperties.defaultH;

    renderer->AddActor(actor);

    renderer->SetViewport(xmins[idx], ymins[idx], xmaxs[idx], ymaxs[idx]);

    renderers[idx] = renderer;
  }

  vtkCamera* camera = nullptr;
  for (const auto& kv : curvatures)
  {
    size_t idx = 0;
    if (kv.second.curvatureName == "Gauss_Curvature")
    {
      idx = 0;
    }
    else
    {
      idx = 1;
    }

    renWin->AddRenderer(renderers[idx]);
    sbWidgets[idx]->On();
    if (idx == 0)
    {
      camera = renderers[idx]->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);
    }
    else
    {
      renderers[idx]->SetActiveCamera(camera);
    }
    renderers[idx]->ResetCamera();
  }
  textWidget->On();

#ifdef HAS_COW
  vtkNew<vtkCameraOrientationWidget> cow;
  cow->SetParentRenderer(renderers[0]);
  cow->On();
#endif

  renWin->Render();
  iRen->Start();

  return EXIT_SUCCESS;
}

namespace {
CurvatureData
GenerateGaussianCurvatures(std::string const& surfaceName, vtkPolyData* source,
                           std::vector<std::string> const& needsAdjusting,
                           bool frequencyTable)
{
  auto scalarBarLabels = 5;

  vtkNew<vtkCurvatures> cc;
  cc->SetInputData(source);
  cc->SetCurvatureTypeToGaussian();
  cc->Update();

  // These will be the active scalars.
  std::string curvature =
      cc->GetOutput()->GetPointData()->GetScalars()->GetName();

  if (std::find(needsAdjusting.begin(), needsAdjusting.end(), surfaceName) !=
      needsAdjusting.end())
  {
    AdjustEdgeCurvatures(cc->GetOutput(), curvature);
  }

  if (surfaceName == "bour")
  {
    // Gaussian curvature is -1/(r(r+1)^4)
    auto r = 1.0;
    auto gauss_curvature = -1.0 / (r * std::pow(r + 1, 4));
    ConstrainCurvatures(cc->GetOutput(), curvature, gauss_curvature,
                        gauss_curvature);
    scalarBarLabels = 1;
  }
  if (surfaceName == "enneper")
  {
    // Gaussian curvature is -4/(1 + r^2)^4
    auto r = 1.0;
    auto gauss_curvature = -4.0 / std::pow((1.0 + r * r), 4);
    ConstrainCurvatures(cc->GetOutput(), curvature, gauss_curvature,
                        gauss_curvature);
    scalarBarLabels = 1;
  }
  if (surfaceName == "cube")
  {
    ConstrainCurvatures(cc->GetOutput(), curvature, 0.0, 0.0);
    scalarBarLabels = 1;
  }
  if (surfaceName == "mobius")
  {
    ConstrainCurvatures(cc->GetOutput(), curvature, 0.0, 0.0);
    scalarBarLabels = 1;
  }
  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);
    scalarBarLabels = 1;
  }

  source->GetPointData()->AddArray(
      cc->GetOutput()->GetPointData()->GetAbstractArray(curvature.c_str()));

  source->GetPointData()->SetActiveScalars(curvature.c_str());
  auto sr = source->GetPointData()->GetScalars(curvature.c_str())->GetRange();
  std::array<double, 2> scalarRange{sr[0], sr[1]};

  // Generate the curvature.
  auto precision = 10;
  auto bands = GetBands(scalarRange.data(), scalarBarLabels, precision, false);

  // Let's do a frequency table.
  auto curvFreq = GetFrequencies(bands, cc->GetOutput());
  AdjustFrequencyRanges(bands, curvFreq);
  if (frequencyTable)
  {
    PrintBandsFrequencies(curvature, bands, curvFreq);
  }

  auto lut = GetDivergingLut();
  lut->SetTableRange(scalarRange.data());

  CurvatureData ret;
  ret.curvatureName = curvature;
  ret.scalarRange = scalarRange;
  ret.lut = lut;
  ret.scalarBarLabels = scalarBarLabels;

  return ret;
}

CurvatureData
GenerateMeanCurvatures(std::string const& surfaceName, vtkPolyData* source,
                       std::vector<std::string> const& needsAdjusting,
                       bool frequencyTable)
{
  // std::string curvature = "Mean_Curvature";

  auto scalarBarLabels = 5;

  vtkNew<vtkCurvatures> cc;
  cc->SetInputData(source);
  cc->SetCurvatureTypeToMean();
  cc->Update();

  // These will be the active scalars.
  std::string curvature =
      cc->GetOutput()->GetPointData()->GetScalars()->GetName();

  if (std::find(needsAdjusting.begin(), needsAdjusting.end(), surfaceName) !=
      needsAdjusting.end())
  {
    AdjustEdgeCurvatures(cc->GetOutput(), curvature);
  }

  if (surfaceName == "bour")
  {
    // Mean curvature is 0
    ConstrainCurvatures(cc->GetOutput(), curvature, 0.0, 0.0);
    scalarBarLabels = 1;
  }
  if (surfaceName == "enneper")
  {
    // Mean curvature is 0
    ConstrainCurvatures(cc->GetOutput(), curvature, 0.0, 0.0);
    scalarBarLabels = 1;
  }
  if (surfaceName == "cube")
  {
    ConstrainCurvatures(cc->GetOutput(), curvature, 0.0, 0.0);
    scalarBarLabels = 1;
  }
  if (surfaceName == "mobius")
  {
    ConstrainCurvatures(cc->GetOutput(), curvature, -0.6, 0.6);
  }
  if (surfaceName == "sphere")
  {
    // Gaussian curvature is 1/r
    auto radius = 10;
    auto gauss_curvature = 1.0 / radius;
    ConstrainCurvatures(cc->GetOutput(), curvature, gauss_curvature,
                        gauss_curvature);
    scalarBarLabels = 1;
  }

  source->GetPointData()->AddArray(
      cc->GetOutput()->GetPointData()->GetAbstractArray(curvature.c_str()));

  source->GetPointData()->SetActiveScalars(curvature.c_str());
  auto sr = source->GetPointData()->GetScalars(curvature.c_str())->GetRange();
  std::array<double, 2> scalarRange{sr[0], sr[1]};

  // Generate the curvature.
  auto precision = 10;
  auto bands = GetBands(scalarRange.data(), scalarBarLabels, precision, false);

  // Let's do a frequency table.
  auto curvFreq = GetFrequencies(bands, cc->GetOutput());
  AdjustFrequencyRanges(bands, curvFreq);
  if (frequencyTable)
  {
    PrintBandsFrequencies(curvature, bands, curvFreq);
  }

  auto lut = GetDivergingLut1();
  lut->SetTableRange(scalarRange.data());

  CurvatureData ret;
  ret.curvatureName = curvature;
  ret.scalarRange = scalarRange;
  ret.lut = lut;
  ret.scalarBarLabels = scalarBarLabels;

  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());
}

// 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
vtkSmartPointer<vtkLookupTable> GetDivergingLut()
{

  vtkNew<vtkColorTransferFunction> ctf;
  ctf->SetColorSpaceToDiverging();
  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);

  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;
}

vtkSmartPointer<vtkLookupTable> GetDivergingLut1()
{
  vtkNew<vtkNamedColors> colors;
  // Colour transfer function.
  vtkNew<vtkColorTransferFunction> ctf;
  ctf->SetColorSpaceToDiverging();
  ctf->AddRGBPoint(0.0, colors->GetColor3d("MidnightBlue").GetRed(),
                   colors->GetColor3d("MidnightBlue").GetGreen(),
                   colors->GetColor3d("MidnightBlue").GetBlue());
  ctf->AddRGBPoint(0.5, colors->GetColor3d("Gainsboro").GetRed(),
                   colors->GetColor3d("Gainsboro").GetGreen(),
                   colors->GetColor3d("Gainsboro").GetBlue());
  ctf->AddRGBPoint(1.0, colors->GetColor3d("DarkOrange").GetRed(),
                   colors->GetColor3d("DarkOrange").GetGreen(),
                   colors->GetColor3d("DarkOrange").GetBlue());

  // Lookup table.
  vtkNew<vtkLookupTable> lut;
  lut->SetNumberOfColors(256);
  for (auto i = 0; i < lut->GetNumberOfColors(); ++i)
  {
    std::array<double, 3> rgb{0.0, 0.0, 0.0};
    ctf->GetColor(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(i, rgba.data());
  }

  return lut;
}

vtkSmartPointer<vtkDiscretizableColorTransferFunction> GetCTF()
{
  // name: Fast, creator: Francesca Samsel, and Alan W. Scott
  // interpolationspace: Lab, space: rgb
  // file name: Fast.json

  vtkNew<vtkDiscretizableColorTransferFunction> ctf;

  ctf->SetColorSpaceToLab();
  ctf->SetScaleToLinear();

  ctf->SetNanColor(0.0, 1.0, 0.0);

  ctf->AddRGBPoint(0, 0.05639999999999999, 0.05639999999999999, 0.47);
  ctf->AddRGBPoint(0.17159223942480895, 0.24300000000000013, 0.4603500000000004,
                   0.81);
  ctf->AddRGBPoint(0.2984914818394138, 0.3568143826543521, 0.7450246485363142,
                   0.954367702893722);
  ctf->AddRGBPoint(0.4321287371255907, 0.6882, 0.93, 0.9179099999999999);
  ctf->AddRGBPoint(0.5, 0.8994959551205902, 0.944646394975174,
                   0.7686567142818399);
  ctf->AddRGBPoint(0.5882260353170073, 0.957107977357604, 0.8338185108985666,
                   0.5089156299842102);
  ctf->AddRGBPoint(0.7061412605695164, 0.9275207599610714, 0.6214389091739178,
                   0.31535705838676426);
  ctf->AddRGBPoint(0.8476395308725272, 0.8, 0.3520000000000001,
                   0.15999999999999998);
  ctf->AddRGBPoint(1, 0.59, 0.07670000000000013, 0.11947499999999994);

  ctf->SetNumberOfValues(9);
  ctf->DiscretizeOff();

  return ctf;
}

vtkSmartPointer<vtkPolyData> GetBour()
{
  auto uResolution = 50;
  auto vResolution = 50;
  vtkNew<vtkParametricBour> surface;

  vtkNew<vtkParametricFunctionSource> source;
  source->SetUResolution(uResolution);
  source->SetVResolution(vResolution);
  source->GenerateTextureCoordinatesOn();
  source->SetParametricFunction(surface);
  source->Update();

  // Build the tangents
  vtkNew<vtkPolyDataTangents> tangents;
  tangents->SetInputConnection(source->GetOutputPort());
  tangents->Update();

  return tangents->GetOutput();
}

vtkSmartPointer<vtkPolyData> GetCube()
{
  vtkNew<vtkCubeSource> surface;

  // Triangulate
  vtkNew<vtkTriangleFilter> triangulation;
  triangulation->SetInputConnection(surface->GetOutputPort());
  // Subdivide the triangles
  vtkNew<vtkLinearSubdivisionFilter> subdivide;
  subdivide->SetInputConnection(triangulation->GetOutputPort());
  subdivide->SetNumberOfSubdivisions(3);
  // Now the tangents
  vtkNew<vtkPolyDataTangents> tangents;
  tangents->SetInputConnection(subdivide->GetOutputPort());
  tangents->Update();
  return tangents->GetOutput();
}

vtkSmartPointer<vtkPolyData> GetEnneper()
{
  auto uResolution = 50;
  auto vResolution = 50;
  vtkNew<vtkParametricEnneper> surface;

  vtkNew<vtkParametricFunctionSource> source;
  source->SetUResolution(uResolution);
  source->SetVResolution(vResolution);
  source->GenerateTextureCoordinatesOn();
  source->SetParametricFunction(surface);
  source->Update();

  // Build the tangents
  vtkNew<vtkPolyDataTangents> tangents;
  tangents->SetInputConnection(source->GetOutputPort());
  tangents->Update();

  return tangents->GetOutput();
}

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> GetMobius()
{
  auto uResolution = 50;
  auto vResolution = 50;
  vtkNew<vtkParametricMobius> surface;
  surface->SetMinimumV(-0.25);
  surface->SetMaximumV(0.25);

  vtkNew<vtkParametricFunctionSource> source;
  source->SetUResolution(uResolution);
  source->SetVResolution(vResolution);
  source->GenerateTextureCoordinatesOn();
  source->SetParametricFunction(surface);
  source->Update();

  // Build the tangents
  vtkNew<vtkPolyDataTangents> tangents;
  tangents->SetInputConnection(source->GetOutputPort());
  tangents->Update();

  vtkNew<vtkTransform> transform;
  transform->RotateX(-90.0);
  vtkNew<vtkTransformPolyDataFilter> transformFilter;
  transformFilter->SetInputConnection(tangents->GetOutputPort());
  transformFilter->SetTransform(transform);
  transformFilter->Update();

  return transformFilter->GetOutput();
}

vtkSmartPointer<vtkPolyData> GetRandomHills()
{
  auto uResolution = 50;
  auto vResolution = 50;
  vtkNew<vtkParametricRandomHills> surface;
  surface->SetRandomSeed(1);
  surface->SetNumberOfHills(30);
  // If you want a plane
  // surface->SetHillAmplitude(0);

  vtkNew<vtkParametricFunctionSource> source;
  source->SetUResolution(uResolution);
  source->SetVResolution(vResolution);
  source->GenerateTextureCoordinatesOn();
  source->SetParametricFunction(surface);
  source->Update();

  // Build the tangents
  vtkNew<vtkPolyDataTangents> tangents;
  tangents->SetInputConnection(source->GetOutputPort());
  tangents->Update();

  vtkNew<vtkTransform> transform;
  transform->Translate(0.0, 5.0, 15.0);
  transform->RotateX(-90.0);
  vtkNew<vtkTransformPolyDataFilter> transformFilter;
  transformFilter->SetInputConnection(tangents->GetOutputPort());
  transformFilter->SetTransform(transform);
  transformFilter->Update();

  return transformFilter->GetOutput();
}

vtkSmartPointer<vtkPolyData> GetSphere()
{
  auto thetaResolution = 32;
  auto phiResolution = 32;
  vtkNew<vtkTexturedSphereSource> surface;
  surface->SetThetaResolution(thetaResolution);
  surface->SetPhiResolution(phiResolution);

  // Now the tangents
  vtkNew<vtkPolyDataTangents> tangents;
  tangents->SetInputConnection(surface->GetOutputPort());
  tangents->Update();
  return tangents->GetOutput();
}

vtkSmartPointer<vtkPolyData> GetTorus()
{
  auto uResolution = 50;
  auto vResolution = 50;
  vtkNew<vtkParametricTorus> surface;

  vtkNew<vtkParametricFunctionSource> source;
  source->SetUResolution(uResolution);
  source->SetVResolution(vResolution);
  source->GenerateTextureCoordinatesOn();
  source->SetParametricFunction(surface);
  source->Update();

  // Build the tangents
  vtkNew<vtkPolyDataTangents> tangents;
  tangents->SetInputConnection(source->GetOutputPort());
  tangents->Update();

  vtkNew<vtkTransform> transform;
  transform->RotateX(-90.0);
  vtkNew<vtkTransformPolyDataFilter> transformFilter;
  transformFilter->SetInputConnection(tangents->GetOutputPort());
  transformFilter->SetTransform(transform);
  transformFilter->Update();

  return transformFilter->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 = {
      {"bour", 0},   {"cube", 1},         {"enneper", 2}, {"hills", 3},
      {"mobius", 4}, {"random hills", 5}, {"sphere", 6},  {"torus", 7}};
  if (available_surfaces.find(surface) == available_surfaces.end())
  {
    std::cout << "The surface is not available." << std::endl;
    std::cout << "Using RandomHills instead." << std::endl;
    surface = "random hills";
  }
  switch (available_surfaces[surface])
  {
  case 0:
    return GetBour();
    break;
  case 1:
    return GetCube();
    break;
  case 2:
    return GetEnneper();
    break;
  case 3:
    return GetHills();
    break;
  case 4:
    return GetMobius();
    break;
  case 5:
    return GetRandomHills();
    break;
  case 6:
    return GetSphere();
    break;
  case 7:
    return GetTorus();
    break;
  }
  return GetRandomHills();
}

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, 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::string const& curvature,
                           std::map<int, std::vector<double>> const& bands,
                           std::map<int, int>& freq)
{

  if (bands.size() != freq.size())
  {
    std::cout << "Bands and frequencies must be the same size." << std::endl;
    return;
  }
  std::string name = curvature;
  std::replace(name.begin(), name.end(), '_', ' ');
  std::string s = fmt::format("Bands & Frequencies:\n{}\n", name);

  size_t idx = 0;
  auto total = 0;
  // 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())
      {
        s += fmt::format("{:4d} [", idx);
      }
      if (q == std::prev(p->second.end()))
      {
        s += fmt::format("{:8.2f}]: {:8d}\n", *q, freq[p->first]);
      }
      else
      {
        s += fmt::format("{:8.2f}], ", *q);
      }
    }
    ++idx;
  }
  // width = 3 * width + 14 = 38;
  s += fmt::format("{:<38s} {:>8d}\n", "Total", total);
  std::cout << s << std::endl;
}

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->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