ElevationBandsWithGlyphs

#include <vtkActor.h>
#include <vtkArrowSource.h>
#include <vtkAxesActor.h>
#include <vtkBandedPolyDataContourFilter.h>
#include <vtkCamera.h>
#include <vtkCameraOrientationWidget.h>
#include <vtkCleanPolyData.h>
#include <vtkColorSeries.h>
#include <vtkColorTransferFunction.h>
#include <vtkDelaunay2D.h>
#include <vtkDoubleArray.h>
#include <vtkElevationFilter.h>
#include <vtkFloatArray.h>
#include <vtkGlyph3D.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 <vtkTextActor.h>
#include <vtkTextProperty.h>
#include <vtkTextRepresentation.h>
#include <vtkTextWidget.h>
#include <vtkTransform.h>
#include <vtkTransformFilter.h>
#include <vtkTriangleFilter.h>
#include <vtkVariant.h>
#include <vtkVariantArray.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 {
/**
 *  Generate elevations over the surface.
 *
 * @param src - The vtkPolyData source.
 * @return elev - The elevations.
 */
vtkSmartPointer<vtkPolyData> GenerateElevations(vtkPolyData* src);

vtkSmartPointer<vtkPolyData> GetHills();
vtkSmartPointer<vtkPolyData> GetParametricHills();
vtkSmartPointer<vtkPolyData> GetParametricTorus();
vtkSmartPointer<vtkPolyData> GetPlane();
vtkSmartPointer<vtkPolyData> GetSphere();
vtkSmartPointer<vtkPolyData> GetTorus();
vtkSmartPointer<vtkPolyData> GetSource(std::string const& source);

vtkSmartPointer<vtkLookupTable> ReverseLUT(vtkLookupTable* lut);

/**
 *  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 AdjustRanges(std::map<int, std::vector<double>>& bands,
                  std::map<int, int>& freq);

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

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

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

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

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

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

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

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

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

  // Create contour edges
  vtkNew<vtkPolyDataMapper> edgeMapper;
  edgeMapper->SetInputData(elevBG.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->SetLookupTable(elevBG.lut1);
  glyphMapper->SetScalarRange(elevBG.scalarRange.data());
  glyphMapper->SetColorModeToMapScalars();
  glyphMapper->ScalarVisibilityOn();
  glyphMapper->SetScalarModeToUsePointFieldData();
  glyphMapper->SelectColorArray("Elevation");

  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 namePositions =
      GetTextPositions(availableSurfaces, VTK_TEXT_LEFT, VTK_TEXT_TOP, 0.25);

  vtkNew<vtkTextProperty> titleTextProperty;
  titleTextProperty->SetColor(colors->GetColor3d("AliceBlue").GetData());
  titleTextProperty->BoldOn();
  titleTextProperty->ItalicOn();
  titleTextProperty->ShadowOn();
  titleTextProperty->SetFontSize(12);
  titleTextProperty->SetJustification(VTK_TEXT_LEFT);

  vtkNew<vtkTextProperty> labelTextProperty;
  labelTextProperty->SetColor(colors->GetColor3d("AliceBlue").GetData());
  labelTextProperty->BoldOff();
  labelTextProperty->ItalicOff();
  labelTextProperty->ShadowOn();
  labelTextProperty->SetFontSize(12);
  labelTextProperty->SetJustification(VTK_TEXT_LEFT);

  vtkNew<vtkTextActor> textActor;
  textActor->SetInput(Title(surfaceName).c_str());
  textActor->SetTextScaleMode(vtkTextActor::TEXT_SCALE_MODE_NONE);

  // Create the text representation. Used for positioning the text actor.
  vtkNew<vtkTextRepresentation> textRepresentation;
  textRepresentation->EnforceNormalizedViewportBoundsOn();
  textRepresentation->GetPositionCoordinate()->SetValue(
      namePositions[surfaceName]["p"].data());
  textRepresentation->GetPosition2Coordinate()->SetValue(
      namePositions[surfaceName]["p2"].data());

  vtkNew<vtkTextWidget> textWidget;
  textWidget->SetRepresentation(textRepresentation);
  textWidget->SetTextActor(textActor);
  textWidget->SetDefaultRenderer(ren);
  textWidget->SetInteractor(iren);
  textWidget->SelectableOff();
  textWidget->EnabledOn();

  auto elevSBP = ScalarBarProperties();
  elevSBP.titleText = elevBG.glyphType + "\n";
  std::replace(elevSBP.titleText.begin(), elevSBP.titleText.end(), '_', '\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;
  elevSBP.positionV = PositionSBWV(elevSBP.numberOfLabels, elevBG.lut1MaxBands);

  auto sbw = MakeScalarBarWidget(elevSBP, titleTextProperty, labelTextProperty,
                                 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 {
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->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();
}

vtkSmartPointer<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] += 1;
        break;
      }
    }
  }
  return freq;
}

void AdjustRanges(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;
}

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

  int maxBands = 0;

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

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

  if (surfaceName == "plane")
  {
    lut->SetNumberOfTableValues(1);
  }
  auto numberOfBands = lut->GetNumberOfTableValues();

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

  auto bands = GetBands(scalarRange, numberOfBands, precision, nearestInteger);
  if (surfaceName == "parametric hills")
  {
    // These are my custom bands.
    // Generated by first running:
    // bands = GetBands(scalarRange, numberOfBands, precision, false);
    // then:
    //  std::vector<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);
    AdjustRanges(bands, freq);
    fmt::println("{:s} Elevation", Title(surfaceName));
    PrintBandsFrequencies(bands, freq);
  }

  // We will use the midpoint of the band as the label.
  std::vector<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.glyphType = "Elevation";
  bg.glyphs = glyphs;
  bg.lutMaxBands = maxBands;
  bg.bcf = bcf;
  bg.lut = lut;
  bg.lutr = ReverseLUT(lut);
  bg.lut1MaxBands = maxBands;
  bg.lut1 = lut1;
  bg.lut1r = ReverseLUT(lut1);
  bg.scalarRange = scalarRange;
  bg.labels = labels;

  return bg;
}

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

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