CurvatureBandsWithGlyphs

#include <vtkActor.h>
#include <vtkArrowSource.h>
#include <vtkAxesActor.h>
#include <vtkBandedPolyDataContourFilter.h>
#include <vtkCamera.h>
#include <vtkCameraOrientationWidget.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 <vtkGenerateIds.h>
#include <vtkGlyph3D.h>
#include <vtkIdList.h>
#include <vtkInteractorStyleTrackballCamera.h>
#include <vtkLookupTable.h>
#include <vtkMaskPoints.h>
#include <vtkNamedColors.h>
#include <vtkNew.h>
#include <vtkOrientationMarkerWidget.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 <vtkTransformFilter.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

namespace fs = std::filesystem;

namespace {

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<vtkLookupTable> ReverseLUT(vtkLookupTable* lut);

struct BandedGlyphs
{
  vtkSmartPointer<vtkGlyph3D> glyphs;
  vtkSmartPointer<vtkBandedPolyDataContourFilter> bcf;
  vtkSmartPointer<vtkLookupTable> lut;
  vtkSmartPointer<vtkLookupTable> lutr;
  vtkSmartPointer<vtkLookupTable> lut1;
  vtkSmartPointer<vtkLookupTable> lut1r;
  std::array<double, 2> scalarRange;
  std::vector<std::string> labels;
};

/** Get curvature glyphs and the corresponding banded polydata filter for the
 * surface.
 *
 * @param surfaceName - The surface name.
 * @param source - A vtkPolyData object corresponding to the source.
 * @param curvature - The curvature type: Gauss_Curvature or Mean_Curvature.
 * @param precision - the precision level.
 * @param frequencyTable - If true, display a frequency table corresponding to
 * the bands.
 * @param nearestInteger - If true, use the nearest integer when generating the
 * bands.
 * @return A struct holding glyphs, bcf, lut, lutr, lut1, lut1r scalar_range,
 *  labels.
 */
BandedGlyphs
GetCurvatureGlyphs(std::string const& surfaceName, vtkPolyData* source,
                   std::string const& curvature, unsigned int precision,
                   bool frequencyTable = false, bool nearestInteger = false);

/**
 * Get elevation glyphs and the corresponding banded polydata filter for the
 * surface.
 *
 * @param surfaceName - The surface name.
 * @param source - A vtkPolyData object corresponding to the source.
 * @param needsAdjusting - Surfaces whose curvatures need to be adjusted
 *  along the edges of the surface or constrained.
 * @param frequencyTable - If true, display a frequency table corresponding to
 * the bands.
 * @param nearestInteger - If true, use the nearest integer when generating the
 * bands.
 * @return A struct holding glyphs, bcf, lut, lutr, lut1, lut1r scalar_range,
 *  labels.
 */
BandedGlyphs GetElevationGlyphs(std::string const& surfaceName,
                                vtkPolyData* source, unsigned int precision,
                                bool frequencyTable = false,
                                bool nearestInteger = false);

/**
 * Generate elevations over the surface.
 *  @param src - the vtkPolyData source.
 *  @return elev - the elevations.
 */
vtkSmartPointer<vtkPolyData> GenerateElevations(vtkPolyData* src);

/**
 * Adjust curvatures along the edges of the surface.
 *
 * This function adjusts curvatures along the edges of the surface by replacing
 *  the value with the average value of the curvatures of points in the
 *  neighborhood.
 *
 * Remember to update the vtkCurvatures object before calling this.
 *
 * @param source - A vtkPolyData object corresponding to the vtkCurvatures
 * object.
 * @param curvatureName - The name of the curvature, "Gauss_Curvature" or
 * "Mean_Curvature".
 * @param epsilon - Curvature values less than this will be set to zero.
 * @return
 */
void AdjustEdgeCurvatures(vtkPolyData* source, std::string const& curvatureName,
                          double const& epsilon = 1.0e-08);

/**
 * Constrain curvatures to the range [lower_bound ... upper_bound].
 *
 * Remember to update the vtkCurvatures object before calling this.
 *
 * @param source - A vtkPolyData object corresponding to the vtkCurvatures
 * object.
 * @param curvatureName - The name of the curvature, "Gauss_Curvature" or
 * "Mean_Curvature".
 * @param lowerBound - The lower bound.
 * @param upperBound - The upper bound.
 * @return
 */
void ConstrainCurvatures(vtkPolyData* source, std::string const& curvatureName,
                         double const& lowerBound = 0.0,
                         double const& upperBound = 0.0);

/**
 *  Glyph the normals on the surface.
 *
 * @param src - The vtkPolyData source.
 * @param scaleFactor - The scale factor for the glyphs.
 * @param reverseNormals - If True the normals on the surface are reversed.
 * @return The glyphs.
 */
vtkNew<vtkGlyph3D> GetGlyphs(vtkPolyData* src, double const& scaleFactor = 1.0,
                             bool const& reverseNormals = false);

/**
 *  Divide a range into bands.
 *
 * @param scalar_range - [min, max] the range that is to be covered by the
 * bands.
 * @param numberOfBands - The number of bands, a positive integer.
 * @param precision - The decimal precision of the bounds..
 * @param nearestInteger - If True then [floor(min), ceil(max)] is used.
 * @return  A map consisting of the band index and [min, midpoint, max] for each
 * band.
 */
std::map<int, std::vector<double>>
GetBands(std::array<double, 2> const& scalarRange, 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 scalar_range - [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 index and [min, midpoint, max] for each
 * band.
 */
std::map<int, std::vector<double>>
GetCustomBands(std::array<double, 2> const& scalarRange,
               int const& numberOfBands,
               std::vector<std::array<double, 2>> const& myBands);

/**
 * Divide a range into integral bands.
 *
 * @param scalar_range - [min, max] the range that is to be covered by the
 * bands.
 * @return A map consisting of the band index and [min, midpoint, max] for each
 * band.
 */
std::map<int, std::vector<double>>
GetIntegralBands(std::array<double, 2> const& scalarRange);

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

/**
 * Print each band and the number of scalars in each band.
 *
 * @param bands: The bands.
 * @param freq: The frequencies.
 * @param precision: The precision for the ranges in each band.
 */
void PrintBandsFrequencies(std::map<int, std::vector<double>> const& bands,
                           std::map<int, int>& freq, int const& precision = 2);

typedef std::map<std::string, std::array<double, 2>> TTextPosition;
typedef std::map<std::string, TTextPosition> TTextPositions;

struct ScalarBarProperties
{
  // 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 numberOfLabels{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);

/**
 * Position the horizontal scalar bar widget.
 *
 *  @param numBands - The number of bands in the scalar bar.
 *  @param maxBands - The maximum number of bands.
 *  @return The scalar bar position.
 */
TTextPosition PositionSBWH(int numBands, int maxBands);

/**
 * Position the vertical scalar bar widget.
 *
 *  @param numBands - The number of bands in the scalar bar.
 *  @param maxBands - The maximum number of bands.
 *  @return The scalar bar position.
 */
TTextPosition PositionSBWV(int numBands, int maxBands);

/**
 * 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.
 *
 * Note: Only works for single-byte characters.
 *
 * @param s: The string to be converted.
 * @param freq: The locale.
 * @return The converted string.
 */
std::string Title(std::string s, const std::locale& loc = std::locale());

/**
 * Adjust the camera parameters.
 *
 * @param surfaceName: The name of the surface.
 * @param ren: The surface renderer.
 */
void AdjustCameraParameters(std::string const& surfaceName, vtkRenderer* ren);

} // 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.");
  auto useOmw{false};
  app.add_flag(
      "-o, --omw", useOmw,
      "Use an OrientationMarkerWidget instead of a CameraOrientationWidget.");

  CLI11_PARSE(app, argc, argv);

  if (surfaceName.empty())
  {
    surfaceName = "parametric 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 == "parametrichills" || surfaceName == "random hills" ||
      surfaceName == "randomhills")
  {
    surfaceName = "parametric hills";
  }
  if (surfaceName == "parametrictorus")
  {
    surfaceName = "parametric torus";
  }

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

  auto it = std::find(availableSurfaces.begin(), availableSurfaces.end(),
                      surfaceName);
  if (it == availableSurfaces.end())
  {
    std::string res = fmt::format("Nonexistent surface: {:s}\n"
                                  "Available 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. "
           "\"parametric hills\"";
    std::cout << res << std::endl;

    return EXIT_FAILURE;
  }

  auto source = GetSource(surfaceName);

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

  std::string curvature{""};

  if (meanCurvature)
  {
    cc->SetCurvatureTypeToMean();
    cc->Update();
    curvature = cc->GetOutput()->GetPointData()->GetScalars()->GetName();
    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")
    {
      // The sphere radius is 1.
      // Mean curvature is 1/r.
      ConstrainCurvatures(cc->GetOutput(), curvature, 1.0, 1.0);
    }
  }
  else
  {
    cc->SetCurvatureTypeToGaussian();
    cc->Update();
    curvature = cc->GetOutput()->GetPointData()->GetScalars()->GetName();
    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")
    {
      // The sphere radius is 1.
      // Gaussian curvature is 1/r^2.
      ConstrainCurvatures(cc->GetOutput(), curvature, 1.0, 1.0);
    }
  }

  BandedGlyphs curvBG = GetCurvatureGlyphs(
      surfaceName, cc->GetOutput(), curvature, 10, frequencyTable, false);

  BandedGlyphs elevBG =
      GetElevationGlyphs(surfaceName, source, 10, frequencyTable, false);

  // ------------------------------------------------------------
  // Create the mappers and actors
  // ------------------------------------------------------------
  vtkNew<vtkNamedColors> colors;

  vtkNew<vtkPolyDataMapper> srcMapper;
  srcMapper->SetInputConnection(curvBG.bcf->GetOutputPort());
  srcMapper->SetScalarRange(curvBG.scalarRange.data());
  srcMapper->SetLookupTable(curvBG.lut);
  srcMapper->SetScalarModeToUseCellData();

  vtkNew<vtkActor> srcActor;
  srcActor->SetMapper(srcMapper);

  // Create contour edges
  vtkNew<vtkPolyDataMapper> edgeMapper;
  edgeMapper->SetInputData(curvBG.bcf->GetContourEdgesOutput());
  edgeMapper->SetResolveCoincidentTopologyToPolygonOffset();

  vtkNew<vtkActor> edgeActor;
  edgeActor->SetMapper(edgeMapper);
  edgeActor->GetProperty()->SetColor(colors->GetColor3d("Black").GetData());

  vtkNew<vtkPolyDataMapper> glyphMapper;
  glyphMapper->SetInputConnection(elevBG.glyphs->GetOutputPort());
  glyphMapper->SetScalarModeToUsePointFieldData();
  glyphMapper->SetColorModeToMapScalars();
  glyphMapper->ScalarVisibilityOn();
  glyphMapper->SelectColorArray("Elevation");
  // Colour by scalars.
  glyphMapper->SetLookupTable(elevBG.lut1);
  glyphMapper->SetScalarRange(elevBG.scalarRange.data());

  vtkNew<vtkActor> glyphActor;
  glyphActor->SetMapper(glyphMapper);

  auto windowWidth = 800;
  auto windowHeight = 800;

  // ------------------------------------------------------------
  // Create the RenderWindow, Renderer and Interactor
  // ------------------------------------------------------------
  vtkNew<vtkRenderer> ren;
  ren->SetBackground(colors->GetColor3d("ParaViewBlueGrayBkg").GetData());
  vtkNew<vtkRenderWindow> renWin;
  renWin->SetSize(windowWidth, windowHeight);
  auto appFn = fs::path((app.get_name())).stem().string();
  renWin->SetWindowName(appFn.c_str());
  renWin->AddRenderer(ren);
  vtkNew<vtkRenderWindowInteractor> iren;
  iren->SetRenderWindow(renWin);

  vtkNew<vtkInteractorStyleTrackballCamera> style;
  iren->SetInteractorStyle(style);

  // Add actors.
  ren->AddViewProp(srcActor);
  ren->AddViewProp(edgeActor);
  ren->AddViewProp(glyphActor);

  // Position the source name according to its length.
  auto textPositions =
      GetTextPositions(availableSurfaces, VTK_TEXT_LEFT, VTK_TEXT_TOP, 0.25);

  vtkNew<vtkTextProperty> titleTextProp;
  // titleTextProp->SetFontFamilyToCourier();
  titleTextProp->SetFontSize(12);
  titleTextProp->SetVerticalJustificationToTop();
  titleTextProp->SetJustificationToLeft();
  titleTextProp->BoldOn();
  titleTextProp->ItalicOn();
  titleTextProp->ShadowOn();
  titleTextProp->SetColor(colors->GetColor3d("AliceBlue").GetData());

  vtkNew<vtkTextProperty> labelTextProp;
  //     labelTextProp->SetFontFamilyToCourier();
  labelTextProp->SetFontSize(12);
  labelTextProp->SetVerticalJustificationToTop();
  labelTextProp->SetJustificationToLeft();
  labelTextProp->BoldOff();
  labelTextProp->ItalicOff();
  labelTextProp->ShadowOn();
  labelTextProp->SetColor(colors->GetColor3d("AliceBlue").GetData());

  auto title = surfaceName;
  std::replace(title.begin(), title.end(), '_', ' ');
  vtkNew<vtkTextActor> textActor;
  textActor->SetInput(Title(title).c_str());
  textActor->SetTextScaleModeToNone();
  textActor->SetTextProperty(titleTextProp);

  // Create the text representation. Used for positioning the text actor.
  vtkNew<vtkTextRepresentation> textRep;
  textRep->EnforceNormalizedViewportBoundsOn();
  textRep->SetPosition(textPositions[surfaceName]["p"][0],
                       textPositions[surfaceName]["p"][1]);
  textRep->SetPosition2(textPositions[surfaceName]["p2"][0],
                        textPositions[surfaceName]["p2"][1]);

  vtkNew<vtkTextWidget> textWidget;
  textWidget->SetRepresentation(textRep);
  textWidget->SetTextActor(textActor);
  textWidget->SetInteractor(iren);
  textWidget->SelectableOff();
  textWidget->ResizableOn();
  textWidget->On();

  // Add scalar bars.
  auto curvSBTtle = curvature;
  std::replace(curvSBTtle.begin(), curvSBTtle.end(), '_', ' ');

  auto curvSBP = ScalarBarProperties();
  curvSBP.titleText = curvSBTtle + "\n";
  curvSBP.numberOfLabels = static_cast<int>(curvBG.labels.size());
  // lut puts the lowest value at the top of the scalar bar.
  // lutr puts the highest value at the top of the scalar bar.
  curvSBP.lut = curvBG.lut;
  curvSBP.orientation = false;
  auto maxBands = 9;
  if (surfaceName == "hills")
  {
    curvSBP.positionH = std::move(PositionSBWH(9, maxBands));
  }
  else if (surfaceName == "parametric hills")
  {
    curvSBP.positionH = std::move(PositionSBWH(7, maxBands));
  }
  else if (surfaceName == "plane" || surfaceName == "sphere")
  {
    curvSBP.positionH = std::move(PositionSBWH(1, maxBands));
  }
  else
  {
    curvSBP.positionH = std::move(PositionSBWH(5, maxBands));
  }

  auto cbw =
      MakeScalarBarWidget(curvSBP, titleTextProp, labelTextProp, ren, iren);

  auto elevSBP = ScalarBarProperties();
  elevSBP.titleText = "Elevation\n";
  elevSBP.numberOfLabels = static_cast<int>(elevBG.labels.size());
  // lut puts the lowest value at the top of the vertical scalar bar.
  // lutr puts the highest value at the top of the vertical scalar bar.
  elevSBP.lut = elevBG.lutr;
  elevSBP.orientation = true;
  maxBands = 8;
  if (surfaceName == "hills" || surfaceName == "parametric hills")
  {
    elevSBP.positionV = std::move(PositionSBWV(8, maxBands));
  }
  else if (surfaceName == "plane")
  {
    elevSBP.positionV = std::move(PositionSBWV(1, maxBands));
  }
  else
  {
    elevSBP.positionV = std::move(PositionSBWV(5, maxBands));
  }

  auto sbw =
      MakeScalarBarWidget(elevSBP, titleTextProp, labelTextProp, ren, iren);

  // Important: The interactor must be set prior to enabling the widget.
  vtkNew<vtkCameraOrientationWidget> cow;
  vtkNew<vtkOrientationMarkerWidget> omw;
  if (useOmw)
  {
    vtkNew<vtkAxesActor> axes;
    double rgba[4]{0.0, 0.0, 0.0, 0.0};
    colors->GetColor("Carrot", rgba);

    omw->SetOrientationMarker(axes);
    omw->SetInteractor(iren);
    omw->SetDefaultRenderer(ren);
    omw->SetOutlineColor(rgba[0], rgba[1], rgba[2]);
    omw->SetViewport(0.8, 0.8, 1.0, 1.0);
    omw->EnabledOn();
    omw->InteractiveOn();
  }
  else
  {
    cow->SetParentRenderer(ren);
    // Enable the widget.
    cow->On();
  }

  ren->ResetCamera();

  AdjustCameraParameters(surfaceName, ren);

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

  return EXIT_SUCCESS;
}

namespace {
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";
  vtkNew<vtkGenerateIds> idFilter;
  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());
}

BandedGlyphs GetCurvatureGlyphs(std::string const& surfaceName,
                                vtkPolyData* source,
                                std::string const& curvature,
                                unsigned int precision, bool frequencyTable,
                                bool nearestInteger)
{
  // The length of the normal arrow glyphs.
  auto scaleFactor = 1.0;
  if (surfaceName == "hills")
  {
    scaleFactor = 0.5;
  }
  if (surfaceName == "sphere")
  {
    scaleFactor = 0.2;
  }

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

  vtkNew<vtkColorSeries> colorSeries;
  if (surfaceName == "hills")
  {
    colorSeries->SetColorScheme(colorSeries->BREWER_DIVERGING_SPECTRAL_7);
  }
  else if (surfaceName == "parametric hills")
  {
    colorSeries->SetColorScheme(colorSeries->BREWER_DIVERGING_SPECTRAL_9);
  }
  else if (surfaceName == "plane" || surfaceName == "sphere")
  {
    colorSeries->SetColorScheme(colorSeries->CITRUS);
  }
  else
  {
    colorSeries->SetColorScheme(colorSeries->BREWER_DIVERGING_SPECTRAL_5);
  }

  vtkNew<vtkLookupTable> lut;
  colorSeries->BuildLookupTable(lut, colorSeries->CATEGORICAL);
  lut->SetNanColor(0, 0, 0, 1);
  lut->SetTableRange(scalarRange.data());

  vtkNew<vtkLookupTable> lut1;
  colorSeries->BuildLookupTable(lut1, colorSeries->ORDINAL);
  lut1->SetNanColor(0, 0, 0, 1);
  lut1->SetTableRange(scalarRange.data());

  auto numberOfBands = lut->GetNumberOfTableValues();
  lut1->SetNumberOfTableValues(numberOfBands);

  auto bands = GetBands(scalarRange, numberOfBands, precision, nearestInteger);

  if (curvature == "Gaussian_Curvature")
  {
    if (surfaceName == "random hills")
    {
      // These are my custom bands.
      // Generated by first running:
      // bands = GetBands(curvScalarRanges, curvNumBands, precision, false);
      // then:
      //  std::vector<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.
      bands = GetCustomBands(scalarRange, numberOfBands, myBands);
    }
    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.
      bands = GetCustomBands(scalarRange, numberOfBands, myBands);
    }
  }

  // Adjust the number of table values and scalar range.
  scalarRange[0] = bands.begin()->second[0];
  scalarRange[1] = std::prev(bands.end())->second[2];
  lut->SetTableRange(scalarRange.data());
  lut->SetNumberOfTableValues(bands.size());
  lut1->SetTableRange(scalarRange.data());
  lut1->SetNumberOfTableValues(bands.size());

  if (frequencyTable)
  {
    auto freq = GetFrequencies(bands, source);
    AdjustFrequencyRanges(bands, freq);
    fmt::println("{:s} {:s}", Title(surfaceName), curvature);
    PrintBandsFrequencies(bands, freq);
  }

  // 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]));
  }
  for (vtkIdType i = 0; i < values->GetNumberOfTuples(); ++i)
  {
    lut->SetAnnotation(i, values->GetValue(i).ToString());
  }

  // Create the contour bands.
  vtkNew<vtkBandedPolyDataContourFilter> bcf;
  bcf->SetInputData(source);
  // Use either the minimum or maximum value for each band.
  int i = 0;
  for (std::map<int, std::vector<double>>::const_iterator p = bands.begin();
       p != bands.end(); ++p)
  {
    bcf->SetValue(i, p->second[2]);
    ++i;
  }
  // We will use an indexed lookup table.
  bcf->SetScalarModeToIndex();
  bcf->GenerateContourEdgesOn();

  // Generate the glyphs on the original surface.
  auto glyphs = GetGlyphs(source, scaleFactor, false);

  BandedGlyphs bg;
  bg.glyphs = glyphs;
  bg.bcf = bcf;
  bg.lut = lut;
  bg.lutr = ReverseLUT(lut);
  bg.lut1 = lut1;
  bg.lut1r = ReverseLUT(lut1);
  bg.scalarRange = scalarRange;
  bg.labels = labels;

  return bg;
}

BandedGlyphs GetElevationGlyphs(std::string const& surfaceName,
                                vtkPolyData* source, unsigned int precision,
                                bool frequencyTable, bool nearestInteger)
{
  // The length of the normal arrow glyphs.
  auto scaleFactor = 1.0;
  if (surfaceName == "hills")
  {
    scaleFactor = 0.5;
  }
  if (surfaceName == "sphere")
  {
    scaleFactor = 0.2;
  }

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

  vtkNew<vtkColorSeries> colorSeries;
  if (surfaceName == "hills" || surfaceName == "parametric hills")
  {
    colorSeries->SetColorScheme(
        colorSeries->BREWER_DIVERGING_BROWN_BLUE_GREEN_8);
  }
  else
  {
    colorSeries->SetColorScheme(
        colorSeries->BREWER_DIVERGING_BROWN_BLUE_GREEN_5);
  }

  vtkNew<vtkLookupTable> lut;
  colorSeries->BuildLookupTable(lut, colorSeries->CATEGORICAL);
  lut->SetNanColor(0, 0, 0, 1);
  lut->SetTableRange(scalarRange.data());

  vtkNew<vtkLookupTable> lut1;
  colorSeries->BuildLookupTable(lut1, colorSeries->ORDINAL);
  lut1->SetNanColor(0, 0, 0, 1);
  lut1->SetTableRange(scalarRange.data());

  auto numberOfBands = lut->GetNumberOfTableValues();
  lut1->SetNumberOfTableValues(numberOfBands);

  auto bands = GetBands(scalarRange, numberOfBands, precision, nearestInteger);

  if (surfaceName == "parametric hills")
  {
    // These are my custom bands.
    // Generated by first running:
    // bands = GetBands(scalarRange, numberOfBands, precision, false);
    // then:
    //  std::vector<int> freq = Frequencies(bands, source);
    //  PrintBandsFrequencies(bands, freq);
    // Finally using the output to create this table:
    std::vector<std::array<double, 2>> myBands = {
        {0, 1.0},   {1.0, 2.0}, {2.0, 3.0}, {3.0, 4.0},
        {4.0, 5.0}, {5.0, 6.0}, {6.0, 7.0}, {7.0, 8.0}};
    // Comment this out if you want to see how allocating
    // equally spaced bands works.
    bands = GetCustomBands(scalarRange, numberOfBands, myBands);
  }

  // Adjust the number of table values and scalar range.
  scalarRange[0] = bands.begin()->second[0];
  scalarRange[1] = std::prev(bands.end())->second[2];
  lut->SetTableRange(scalarRange.data());
  lut->SetNumberOfTableValues(bands.size());
  lut1->SetTableRange(scalarRange.data());
  lut1->SetNumberOfTableValues(bands.size());

  if (frequencyTable)
  {
    auto freq = GetFrequencies(bands, source);
    AdjustFrequencyRanges(bands, freq);
    fmt::println("{:s} Elevation", Title(surfaceName));
    PrintBandsFrequencies(bands, freq);
  }

  // We will use the midpoint of the band as the label.
  std::vector<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]));
  }
  for (vtkIdType i = 0; i < values->GetNumberOfTuples(); ++i)
  {
    lut->SetAnnotation(i, values->GetValue(i).ToString());
  }

  // Create the contour bands.
  vtkNew<vtkBandedPolyDataContourFilter> bcf;
  bcf->SetInputData(source);
  // Use either the minimum or maximum value for each band.
  int i = 0;
  for (std::map<int, std::vector<double>>::const_iterator p = bands.begin();
       p != bands.end(); ++p)
  {
    bcf->SetValue(i, p->second[2]);
    ++i;
  }
  // We will use an indexed lookup table.
  bcf->SetScalarModeToIndex();
  bcf->GenerateContourEdgesOn();

  // Generate the glyphs on the original surface.
  auto glyphs = GetGlyphs(source, scaleFactor, false);

  BandedGlyphs bg;
  bg.glyphs = glyphs;
  bg.bcf = bcf;
  bg.lut = lut;
  bg.lutr = ReverseLUT(lut);
  bg.lut1 = lut1;
  bg.lut1r = ReverseLUT(lut1);
  bg.scalarRange = scalarRange;
  bg.labels = labels;

  return bg;
}

vtkSmartPointer<vtkPolyData> GenerateElevations(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[2] = 0.0;
    bounds[3] = 1.0e-8;
  }
  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.
  vtkNew<vtkPoints> points;
  for (auto i = 0; i < xRes; ++i)
  {
    auto x = xMin + i * dx;
    for (auto j = 0; j < yRes; ++j)
    {
      auto y = yMin + j * dy;
      points->InsertNextPoint(x, y, 0);
    }
  }

  // Add the grid points to a polydata object.
  vtkNew<vtkPolyData> plane;
  plane->SetPoints(points);

  // Triangulate the grid.
  vtkNew<vtkDelaunay2D> delaunay;
  delaunay->SetInputData(plane);
  delaunay->Update();

  auto polydata = delaunay->GetOutput();

  vtkNew<vtkDoubleArray> elevation;
  elevation->SetNumberOfTuples(points->GetNumberOfPoints());

  //  We define the parameters for the hills here.
  // [[0: x0, 1: y0, 2: x variance, 3: y variance, 4: amplitude]...]
  std::vector<std::array<double, 5>> hd{{-2.5, -2.5, 2.5, 6.5, 3.5},
                                        {2.5, 2.5, 2.5, 2.5, 2},
                                        {5.0, -2.5, 1.5, 1.5, 2.5},
                                        {-5.0, 5, 2.5, 3.0, 3}};
  // 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};
  for (auto i = 0; i < points->GetNumberOfPoints(); ++i)
  {
    auto x = polydata->GetPoint(i);
    for (size_t j = 0; j < hd.size(); ++j)
    {
      xx[0] = std::pow(x[0] - hd[j][0] / hd[j][2], 2.0);
      xx[1] = std::pow(x[1] - hd[j][1] / hd[j][3], 2.0);
      x[2] += hd[j][4] * std::exp(-(xx[0] + xx[1]) / 2.0);
    }
    polydata->GetPoints()->SetPoint(i, x);
    elevation->SetValue(i, x[2]);
  }

  vtkNew<vtkFloatArray> textures;
  textures->SetNumberOfComponents(2);
  textures->SetNumberOfTuples(2 * polydata->GetNumberOfPoints());
  textures->SetName("Textures");

  for (auto i = 0; i < xRes; ++i)
  {
    float tc[2];
    tc[0] = i / (xRes - 1.0);
    for (auto j = 0; j < yRes; ++j)
    {
      // tc[1] = 1.0 - j / (yRes - 1.0);
      tc[1] = j / (yRes - 1.0);
      textures->SetTuple(static_cast<vtkIdType>(i) * yRes + j, tc);
    }
  }

  polydata->GetPointData()->SetScalars(elevation);
  polydata->GetPointData()->GetScalars()->SetName("Elevation");
  polydata->GetPointData()->SetTCoords(textures);

  vtkNew<vtkPolyDataNormals> normals;
  normals->SetInputData(polydata);
  normals->SetInputData(polydata);
  normals->SetFeatureAngle(30);
  normals->SplittingOff();

  vtkNew<vtkTransform> tr1;
  tr1->RotateX(-90);

  vtkNew<vtkTransformFilter> tf1;
  tf1->SetInputConnection(normals->GetOutputPort());
  tf1->SetTransform(tr1);
  tf1->Update();

  return tf1->GetPolyDataOutput();
}

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<vtkTransformFilter> transformFilter;
  transformFilter->SetInputConnection(source->GetOutputPort());
  transformFilter->SetTransform(transform);
  transformFilter->Update();

  return transformFilter->GetPolyDataOutput();
}

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<vtkTransformFilter> transformFilter;
  transformFilter->SetInputConnection(source->GetOutputPort());
  transformFilter->SetTransform(transform);
  transformFilter->Update();

  return transformFilter->GetPolyDataOutput();
}

vtkSmartPointer<vtkPolyData> GetPlane()
{
  vtkNew<vtkPlaneSource> source;
  source->SetOrigin(-10.0, -10.0, 0.0);
  source->SetPoint2(-10.0, 10.0, 0.0);
  source->SetPoint1(10.0, -10.0, 0.0);
  source->SetXResolution(10);
  source->SetYResolution(10);
  source->Update();

  vtkNew<vtkTransform> transform;
  transform->Translate(0.0, 0.0, 0.0);
  transform->RotateX(-90.0);
  vtkNew<vtkTransformFilter> 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(1.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->SetPhiResolution(64);
  source->SetThetaResolution(64);
  source->SetThetaRoundness(1);
  source->SetThickness(0.5);
  source->SetSize(10);
  source->SetToroidal(1);
  source->Update();

  // The quadric is made of strips, so pass it through a triangle filter as
  // the curvature filter only operates on polys
  vtkNew<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); });
  if (surface == "hills")
  {
    return GetHills();
  }
  else if (surface == "parametric hills")
  {
    return GetParametricHills();
  }
  else if (surface == "parametric torus")
  {
    return GetParametricTorus();
  }
  else if (surface == "plane")
  {
    return GenerateElevations(GetPlane());
  }
  else if (surface == "sphere")
  {
    return GenerateElevations(GetSphere());
  }
  else if (surface == "torus")
  {
    return GenerateElevations(GetTorus());
  }
  std::cout << "The surface is not available." << std::endl;
  std::cout << "Using parametric hills instead." << std::endl;
  return GetParametricHills();
}

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(vtkPolyData* src, double const& scaleFactor,
                             bool const& reverseNormals)
{
  // Sometimes the contouring algorithm can create a volume whose gradient
  // vector and ordering of polygon(using the right hand rule) are
  // inconsistent. vtkReverseSense cures this problem.
  vtkNew<vtkReverseSense> reverse;
  vtkNew<vtkMaskPoints> maskPts;
  maskPts->SetOnRatio(5);
  maskPts->RandomModeOn();
  if (reverseNormals)
  {
    reverse->SetInputData(src);
    reverse->ReverseCellsOn();
    reverse->ReverseNormalsOn();
    maskPts->SetInputConnection(reverse->GetOutputPort());
  }
  else
  {
    maskPts->SetInputData(src);
  }

  // Source for the glyph filter
  vtkNew<vtkArrowSource> arrow;
  arrow->SetTipResolution(16);
  arrow->SetTipLength(0.3);
  arrow->SetTipRadius(0.1);

  vtkNew<vtkGlyph3D> glyph;
  glyph->SetSourceConnection(arrow->GetOutputPort());
  glyph->SetInputConnection(maskPts->GetOutputPort());
  glyph->SetVectorModeToUseNormal();
  glyph->SetScaleFactor(scaleFactor);
  glyph->SetColorModeToColorByVector();
  glyph->SetScaleModeToScaleByVector();
  glyph->OrientOn();
  glyph->Update();

  return glyph;
}

std::map<int, std::vector<double>>
GetBands(std::array<double, 2> const& scalarRange, 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 ((scalarRange[1] < scalarRange[0]) || (numberOfBands <= 0))
  {
    return bands;
  }
  double x[2];
  for (int i = 0; i < 2; ++i)
  {
    x[i] = scalarRange[i];
  }
  if (nearestInteger)
  {
    x[0] = std::floor(x[0]);
    x[1] = std::ceil(x[1]);
  }
  double dx = (x[1] - x[0]) / static_cast<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(std::array<double, 2> const& scalarRange,
               int const& numberOfBands,
               std::vector<std::array<double, 2>> const& myBands)
{
  std::map<int, std::vector<double>> bands;
  if ((scalarRange[1] < scalarRange[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 (scalarRange[0] < myBands[idx][1] && scalarRange[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 (scalarRange[1] < myBands[idx][1] && scalarRange[1] >= myBands[idx][0])
    {
      idxMax = static_cast<int>(idx);
      break;
    }
  }

  // Set the minimum to match the range minimum.
  x[idxMin][0] = scalarRange[0];
  x[idxMax][1] = scalarRange[1];
  for (int i = idxMin; i < idxMax + 1; ++i)
  {
    std::vector<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(std::array<double, 2> const& scalarRange)
{
  if (scalarRange[1] < scalarRange[0])
  {
    std::map<int, std::vector<double>> bands;
    return bands;
  }
  std::array<double, 2> x{0, 0};
  for (int i = 0; i < 2; ++i)
  {
    x[i] = scalarRange[i];
  }
  x[0] = std::floor(x[0]);
  x[1] = std::ceil(x[1]);
  int numberOfBands = static_cast<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;
  std::string res{""};
  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())
      {
        res += fmt::format("{:4d} [", idx);
      }
      if (q == std::prev(p->second.end()))
      {
        res += fmt::format("{:{}.{}f}]: {:{}d}\n", *q, width, precision,
                           freq[p->first], width);
      }
      else
      {
        res += fmt::format("{:{}.{}f}, ", *q, width, precision);
      }
    }
    ++idx;
  }
  auto width1 = 3 * width + 13;
  fmt::println("{:s}\n{:<{}s}{:>{}d}", res, "Total", width1, total, width);
  std::cout << std::endl;
}

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.numberOfLabels);
  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;
}

TTextPosition PositionSBWH(int numBands, int maxBands)
{
  maxBands = std::abs(maxBands);
  numBands = std::abs(numBands);
  if (numBands > maxBands)
  {
    numBands = maxBands;
  }
  if (numBands == 0)
  {
    numBands = 1;
  }

  // Origin of the scalar bar.
  std::array<double, 2> xy0{0.125, 0.05};
  // Width and height of the scalar bar.
  std::array<double, 2> dxy{0.75, 0.1};

  if (numBands >= maxBands)
    return {{"p", xy0}, {"p2", dxy}};

  auto dx = dxy[0] - xy0[0] * numBands / maxBands;
  dxy[0] = dxy[0] * numBands / maxBands;
  if (numBands == 1)
  {
    xy0[0] = 0.5 - dx * numBands / (maxBands * 2);
  }
  else
  {
    xy0[0] = 0.5 - dx * (numBands + 1) / (maxBands * 2);
  }
  return {{"p", xy0}, {"p2", dxy}};
}

TTextPosition PositionSBWV(int numBands, int maxBands)
{
  maxBands = std::abs(maxBands);
  numBands = std::abs(numBands);
  if (numBands > maxBands)
  {
    numBands = maxBands;
  }
  if (numBands == 0)
  {
    numBands = 1;
  }

  // Origin of the scalar bar.
  std::array<double, 2> xy0{0.9, 0.25};
  // Width and height of the scalar bar.
  std::array<double, 2> dxy{0.08, 0.5};

  if (numBands >= maxBands)
    return {{"p", xy0}, {"p2", dxy}};

  auto dx = dxy[1] - xy0[1] * numBands / maxBands;
  dxy[1] = dxy[1] * numBands / maxBands;
  if (numBands == 1)
  {
    xy0[1] = 0.5 - dx * numBands / (maxBands * 2);
  }
  else
  {
    xy0[1] = 0.5 - dx * (numBands + 1) / (maxBands * 2);
  }
  return {{"p", xy0}, {"p2", dxy}};
}

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(std::string s, const std::locale& loc)
{
  auto newWord = true;
  for (char& c : s)
  {
    c = static_cast<char>(newWord ? std::toupper(c, loc)
                                  : std::tolower(c, loc));
    newWord = std::isspace(c, loc);
  }
  return s;
}

void AdjustCameraParameters(std::string const& surfaceName, vtkRenderer* ren)
{
  if (surfaceName == "hills")
  {
    auto camera = ren->GetActiveCamera();
    camera->SetPosition(16.3424, 19.8311, 0.46492);
    camera->SetFocalPoint(0.209609, 0.432443, -1.18699);
    camera->SetViewUp(-0.755535, 0.640179, -0.13906);
    camera->SetDistance(25.2845);
    camera->SetClippingRange(13.1133, 37.6179);
  }
  else if (surfaceName == "parametric hills")
  {
    auto camera = ren->GetActiveCamera();
    camera->SetPosition(10.9299, 59.1505, 24.9823);
    camera->SetFocalPoint(2.21692, 7.97545, 7.75135);
    camera->SetViewUp(-0.230136, 0.345504, -0.909761);
    camera->SetDistance(54.6966);
    camera->SetClippingRange(36.3006, 77.9852);
  }
  else if (surfaceName == "parametric torus")
  {
    auto camera = ren->GetActiveCamera();
    camera->SetPosition(-1.38419, 24.2883, 34.9246);
    camera->SetFocalPoint(-2.07248e-07, 3.63658e-06, 0.016056);
    camera->SetViewUp(0.010284, 0.821007, -0.570825);
    camera->SetDistance(42.5493);
    camera->SetClippingRange(25.2917, 64.5115);
  }
  else if (surfaceName == "plane")
  {
    auto camera = ren->GetActiveCamera();
    camera->SetPosition(-0.516003, 22.5763, 51.9171);
    camera->SetFocalPoint(-5.77108e-08, 0.500002, 5.80651e-06);
    camera->SetViewUp(-0.000956134, 0.920254, -0.391321);
    camera->SetDistance(56.4182);
    camera->SetClippingRange(36.7854, 81.268);
  }
  else if (surfaceName == "torus")
  {
    auto camera = ren->GetActiveCamera();
    camera->SetPosition(-2.02659, 35.5605, 51.1256);
    camera->SetFocalPoint(0, 0, 0.0160508);
    camera->SetViewUp(0.010284, 0.821007, -0.570825);
    camera->SetDistance(62.2964);
    camera->SetClippingRange(38.14, 92.8545);
  }
}

} // namespace