/***********************************************************************
pftrail.cpp v1.01 (31 March 2020)

Copyright 2020 Herman Johannes Haverkort

Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at

  http://www.apache.org/licenses/LICENSE-2.0

Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.

This file is part of the pftrail package. 
The package includes the following files:

* pftrail.cpp
* ifs-classics
* ifs-inventions
* generate-preamble.cpp
* colours
* postamble
* manual.pdf

New versions of the package may be available from 
http://spacefillingcurves.net or http://herman.haverkort.net.

***********************************************************************/

const char* pftrailInfo = 
"pftrail v1.01 by Herman Haverkort.";
 
#include <string>
using std::string;
#include <string.h>
#include <sstream>
using std::stringstream;
#include <iostream>
using std::istream;
using std::ostream;
using std::cout;
using std::cerr;
using std::endl;
using std::fixed;
#include <fstream>
using std::ifstream;
#include <iomanip>
using std::setprecision;
#include <cstdlib>
#include <utility>
#include <stack>
#include <vector>
#include <algorithm>
using std::min;
using std::max;
using std::swap;
#include <cmath>
using std::abs;

/*********************************************************
* Tunable constants etc.                                 *
*********************************************************/

// for polynomial zoom, the precision of doubles may be limiting:
typedef long double DOUBLE;

// precision of the coordinates in the output:
const int outputPrecision = 5;

// constant used to:
// - not choke on computation-induced rounding errors in computeMinDisk;
// - not choke on computation- or IFS-specification-induced rounding 
//   errors in the IFS parser, more specifically the part that checks
//   whether the IFS indeed generates a curve of fractal dimension
//   two, and which part of the curve fills how much of plane;
// - round coordinates that are probably zero to actually zero when
//   the IFS parser converts the input to Cartesian coordinates:
const DOUBLE practicallyZero = 0.0000001;

// maxMeshFill says with how many vertices (modulo the addition of
// another polygon) a mesh is considered to be full:
const int maxMeshFill = 65000;

// cornerDepth determines how tall/flat the model is;
// its value is half of the ratio between the edge length 
// (or equivalently: the short diagonal) and the height (1)
// of the rendering grid:
const DOUBLE cornerDepth = 1.06066; 

// where to put the base of the model (which is normally 
// invisible, because the lower ref. plane at elevation 0 
// hides the base):
const DOUBLE foundation = -0.5;

// the minimum elevation difference for which parapets
// are created:
const DOUBLE noParapetDrop = 0.001;

// the minimum elevation difference for which parapets
// have their maximum height: 
const DOUBLE maxParapetDrop = 0.02;

// the next constant controls the maximum horizontal thickness of
// parapets; the higher the number, the thinner the parapets.
// setting it too large could create artefacts from rounding errors
// in the output
const DOUBLE minSurfaceToParapetRatio = 10;

// the zoom exponent for vertical distance as a function
// of the zoom exponent for horizontal distance:
DOUBLE elevationZoomExponent(DOUBLE zoomExponent)
{
  return max((DOUBLE) 1.0, (zoomExponent * 1.5 - 0.5));
} 

/*********************************************************
* Input/output headers                                   *
*********************************************************/

enum ErrorLevel 
{ 
  warning, 
  userError, 
  ifsParsingError, 
  usageError, 
  internalError 
};

class Complaint 
// an object to which an error message can be witten with <<. 
// when the object is destroyed, the error message is produced
// and, if it is not a mere warning, the program halts.
{
    private:
  static const char* errorLabel[];
  stringstream complaint;
  ErrorLevel errorLevel;
    public:
  Complaint(ErrorLevel _errorLevel = userError);
  Complaint(ErrorLevel _errorLevel, const string& expected, const string& read);
  ~Complaint();
  template <typename T> Complaint& operator<<(const T& t);
};

class CommentedInput 
// an object that reads a file and while skipping comments.
// Comments are from // to the end of the line, or 
// between /* and */; in the latter case, /* has to be at
// the beginning of the line.
{
    private:
  ifstream source;
  stringstream buffer;
  void fillBuffer();

    public:
  CommentedInput(string& fileName);
  bool fail() const;
  template < class T > CommentedInput& operator>>(T& t);
  void getLine(string& s);
  void getLowerCaseWord(string& s);

  // The following method scans the input and stops after 
  // reading two consecutive words keyword and name 
  // (case-insensitive). Returns true iff the combination
  // was found. 
  // More precisely: after an occurrence of keyword 
  // the next word is compared to name; if it is not a 
  // match, the word is ignored. So, keyword keyword name
  // is not a match, as the first occurrence of keyword
  // "consumes" the second.
  bool find(const char* keyword, const string& name);
};

void makeLowerCase(string& s);

// convert a string to a number, verifying that the string 
// does not contain any further content:
DOUBLE strtofstop(const char* str);
int strtoistop(const char* str);


/*********************************************************
* Geometry headers                                       *
*********************************************************/

struct Point
{
  DOUBLE x; DOUBLE y;
  Point(DOUBLE _x = 0, DOUBLE _y = 0);
  Point operator+(const Point& p) const;
  Point& operator+=(const Point& p);
  Point operator-(const Point& p) const;
  Point& operator-=(const Point& p);
  Point operator-() const;
  DOUBLE operator*(const Point& p) const; // inner product
  Point operator*(const DOUBLE d) const;
  Point operator/(const DOUBLE d) const;
  DOUBLE sqLength() const;
  DOUBLE length() const;
};
ostream& operator<<(ostream& s, const Point& p);

struct ElevatedPoint: public Point
{
  DOUBLE h;
  ElevatedPoint();
  ElevatedPoint(const Point& _p, DOUBLE _h);
  ElevatedPoint(DOUBLE _x, DOUBLE _y, DOUBLE _h);
  ElevatedPoint operator+(const ElevatedPoint& e) const;
  ElevatedPoint& operator+=(const ElevatedPoint& e);
  ElevatedPoint& operator-=(const Point& p);
  ElevatedPoint& operator-=(const ElevatedPoint& e);
};
ostream& operator<<(ostream& s, const ElevatedPoint& p);

enum SampleCategory { sampleFill, sampleJump, sampleJumpCore };

struct AnnotatedPoint: public ElevatedPoint
{
  SampleCategory a;
  std::pair< DOUBLE, DOUBLE > range;
  AnnotatedPoint();
  AnnotatedPoint(const ElevatedPoint& _e, SampleCategory _a, DOUBLE r = 0);
  AnnotatedPoint(const Point& _p, DOUBLE _h, SampleCategory _a, DOUBLE r = 0);
};

// to compute the radius of a smallest enclosing circle:
DOUBLE radius(const std::vector< AnnotatedPoint >& points);

struct Rotation
// a rotation with scaling in 2D, defaults to identity transformation
{
  Point x; Point y; // corresponding rows of the rotation matrix
  Rotation(const Point& _x = Point(1,0), const Point& _y = Point(0,1));
  Rotation transpose() const;

  // composition *this after r:
  Rotation operator*(const Rotation& r) const; 
  Rotation operator*=(const Rotation& r);

  // apply to point:
  Point operator*(const Point& p) const;

  // (squared) scaling factor, assuming the transformation is indeed
  // a similarity transformation (not just any linear transformation):       
  DOUBLE sqScale() const; 
  DOUBLE scale() const;
};

struct Similarity
// similarity transformations that consist of rotation, scaling, 
// translation and/or reversal (changing the sign of the 3rd coordinate);
// these can be understood as a 4x4 matrix of which:
// - the upper left quadrant encodes rotation and scaling of x, y;
// - the last elements of the first rows encode translation of x, y, t
// - the third element of the third row encodes scaling/reversal of t
// - the fourth element of the fourth row is 1
// - all other elements are zero.
{
  Rotation rotation;         // upper left quadrant of the 4x4 matrix
  ElevatedPoint translation; // last elements of first three rows
  DOUBLE scale;              // third element of third row
  Similarity(
    const Rotation& _rotation = Rotation(), 
    const ElevatedPoint& _translation = ElevatedPoint(), 
    DOUBLE _scale = 1
  );

  // composition *this after s:
  Similarity operator*(const Similarity& s) const;
  Similarity& operator*=(const Similarity& s);

  // composition s after *this (store result in *this):
  Similarity& apply(const Similarity& s);
  
  // apply to point:
  Point operator()(const Point& p) const;
  ElevatedPoint operator()(const ElevatedPoint& e) const;

  // print with labels left and right of the first row:
  void print( 
    std::ostream& o, 
    const string& leftLabel, 
    const string& rightLabel = string()
  ) const; 
};

// useful constants
const DOUBLE SQRTTHIRD = sqrt((DOUBLE) 1.0/ (DOUBLE) 3.0);
const DOUBLE SQRTTHREE = sqrt((DOUBLE) 3.0);

template < class T >
class TrianglePoint
// points with coordinates on three axes at 120 degree angles 
// (or equivalent: 3D points that represent their orthogonal 
// projections on the plane x + y + z = 0);
// with conversions from/to 2D Cartesian points 
// note: in the three coordinates' system, lines x = a and x = b
// are a distance b-a apart when measured along the y- or z-axis;
// and a distance sqrt(3/4) (b-a) apart when measured orthogonal
// to the x-axis. 
{
    private:
  T x[3];
  static const Rotation TriToCar;

    public:
  // constructor sets x[s mod 3] = a, x[s+1 mod 3] = b, x[s+2 mod 3] = -(a+b):
  TrianglePoint(int s = 0, T a = 0, T b = 0); 
  TrianglePoint(const Point& p);
  T operator[](int i) const;
  operator Point() const;
};

// to round a TrianglePoint to the nearest integer TrianglePoint:
TrianglePoint< int > round(const TrianglePoint< DOUBLE >& t);

// to compute the minimum axis-aligned bounding diamond (that is, 
// roughly, what will be the rendering grid) for a set of points,
// the result is stored in the 2nd and 3rd arguments:
void findMinimumBoundingDiamond(
  const std::vector< AnnotatedPoint >& pp,
  Point& centre,
  DOUBLE& edgeLength
);


/*********************************************************
* IFS headers (parsing, analysis and sampling)           *
*********************************************************/

class Options;

struct Segment 
{
  Similarity transformation;
  int generator;
  Segment(
    Similarity _transformation = Similarity(),
    int _generator = 0
  );
};

struct Generator 
{
  std::vector< Segment > segments;
  DOUBLE maxRadius;
  DOUBLE maxSampleLength;
  DOUBLE maxSqSampleLength;
};

enum BridgeStyle 
{ 
  noBridges, 
  pointedBridges, 
  rectangularBridges 
};

class IFS
{
    private:
  std::vector< Generator > generators;
  bool continuous;
  DOUBLE diameter; // lower bound on ratio image diameter : generator length
  DOUBLE curveSampleRadius; // radius of region that is covered by a sample
  BridgeStyle bridgeStyle;
  DOUBLE bridgeSampleDistance;
  DOUBLE bridgeWidthFactor;
  DOUBLE bridgeCoreWidth;
  ElevatedPoint centre; // centre for polynomial zooming
  DOUBLE zoomExponent;
  DOUBLE invZoomExponent;
  DOUBLE invElevationZoomExponent;
  DOUBLE minScaledElev;
  DOUBLE invScaledElevRange;
  stringstream definition;

  // some helper functions for the parser
  void skip(CommentedInput& text, const char* s) const;
  void skip(CommentedInput& text, const char c) const;
  string generatorName(int g) const;
  typedef std::vector< DOUBLE > Equation;
  static void writeAreaEquation(ostream& s, const Equation& e);
  static void writeAreaEquations(ostream& s, const std::vector< Equation >& ee);

    public:
  void parse(CommentedInput& text);
  bool isContinuous() const;
  DOUBLE diameterLowerBound() const;  

  // calculate coordinates of point at start of subsegment indexed by string:
  ElevatedPoint operator()(std::string index) const; 

    private:
  // some helper functions for the sample point generator:
  ElevatedPoint zoom(const ElevatedPoint& p);
  void buildBridge(
    std::vector< AnnotatedPoint >& pp, 
    Point a, Point b, DOUBLE h
  );
  void generatePoints(
    std::vector< AnnotatedPoint >& pp, 
    int generator,
    const Similarity& s,
    bool& jump // must be set to true if jumping in; this includes
               // the initial call; the method will set the variable
               // to true (or false) if (not) jumping on the way out.
  );

    public:
  void applySketchSettings(DOUBLE sampleLength);
  void applyUserSettings(const Options& options);
  void generatePoints(
    std::vector< AnnotatedPoint >& pp, // vector to store the output
    const Similarity& s = Similarity() // transformation to apply
  );
  void report(ostream& o);
};

enum CoordinateSystem 
{
  planarCoordinates, // Cartesian
  cubicCoordinates   // 3D points representing projection on plane x+y+z = 0
};

void readSegmentParameters(
  // input:
  stringstream& text,
  // output:
  DOUBLE& x,
  DOUBLE& y,
  int& g, // generator, must be preset to number of generators minus 1,
          // and will only be read if larger than 0
  int& d,
  int& r,
  // context:
  CoordinateSystem c = planarCoordinates
);


/*********************************************************
* Options headers                                        *
*********************************************************/

struct OptionSyntax 
{
  int rank;
  string id;
  bool requiresArgument;
  OptionSyntax(int _rank, string _id, bool _requiresValue = false);
};

class View
{
    public:
  Similarity orientation;
  bool ridge;
  DOUBLE ridgeLine;
  bool defaultCentre;
  ElevatedPoint centre;
  DOUBLE altitude;
  DOUBLE azimuth;
  DOUBLE distance;

  View();
  void reset();
  void parse(CommentedInput& text, const IFS& ifs);
  void report(ostream& s) const;

    private:
  static const OptionSyntax syntax[];
};

struct Options
{
    public:
  string ifsFileName; string ifsName; IFS ifs; 
  string cfgFileName; string cfgName; View view;
  DOUBLE accuracy;
  BridgeStyle bridgeStyle;
  DOUBLE bridgeWidth; 
  DOUBLE bridgeCoreWidth;
  DOUBLE clearance;
  DOUBLE strength;
  int nrColourClasses;
  int erosionRadius;
  DOUBLE resolution;
  DOUBLE parapet;
  bool supported; // supported by rock, as opposed to floating
  DOUBLE zoomExponent;
  static int verbosity;

  Options();
  void reset();
  void parse(const int argc, const char* argv[]);
  void report(ostream& s) const;
  static void usage();

    private:
  static const OptionSyntax syntax[];
  static void cutId(const string id, string& fileName, string& name);
};


/*********************************************************
* Hex Grid headers (model creation)                      *
*********************************************************/

// in rendering the curve, four coordinate systems are involved,
// with coinciding origins:
// - sample space (x,y,h) (equal to model space, see below,
//   up to a scaling factor such that it fits on the grid)
// - grid space (row,col,h), where row and col are in 
//   [-gir, +gir]. (gir = gridIndexRadius)
// - model space (x,y,h): (c = cornerDepth)
//   grid (-(gir-1),-(gir-1),h) -> model ( 0,  sqrt(3) c, h)
//   grid (-(gir-1),+(gir-1),h) -> model (-c,  0,         h)
//   grid (+(gir-1),+(gir-1),h) -> model ( 0, -sqrt(3) c, h)
//   grid (+(gir-1),-(gir-1),h) -> model ( c,  0,         h)
// - output space: model (x,y,h) -> output (-y,h,-x)

// neighbours of hex cells are numbered 0 to 5, in 
// counterclockwise order, starting with the one in 
// direction (+1,0) in grid space, which is direction 
// (1, -sqrt(3)) in sample/model space

class HexCell;

struct Slab
// Models a solid object that covers an entire grid cell between
// specified bottom and top elevations
{
  DOUBLE bottom, top;
  Slab();
  Slab(DOUBLE _bottom, DOUBLE _top);

  // decide what walls (from top to bottom) and what parapet
  // (i.e. its height) needs to be created for this slab on the
  // boundary with the specified neighbouring grid cell:
  void decideWallAndParapet(
    const HexCell& neighbour, 
    std::vector< DOUBLE >& wallStartStop, 
    DOUBLE& parapet
  ) const;
};

struct Sample
// An elevation sample. Collected by the constructor of HexGrid
// (see below) and passed to HexCell to calculate slabs
{
  DOUBLE h;
  SampleCategory a;
  std::pair< DOUBLE, DOUBLE > range;
  Sample();
  Sample(AnnotatedPoint p); 
  Sample(DOUBLE _h);
  DOUBLE weight() const;   // size of range
  DOUBLE priority() const; // weight of jump core samples, 0 otherwise 
  bool operator<(const Sample& s) const; // sorts top-down
};

class HexGrid;

class HexCell
// a cell in the hex grid, storing a sequence of slabs, order
// from top to bottom
{
  std::vector< Slab > slabs;
  friend class HexGrid;

  public:
  HexCell();
  typedef std::vector< Slab >::const_iterator SlabIterator;
  SlabIterator begin() const;
  SlabIterator end() const;

  // the following method turns the requested floor levels into
  // a sequence of slabs from top to bottom for the solid and
  // floating drawing styles:

  void calculateSlabs(
    std::vector< Sample >& samples, // must be pre-sorted top-down
    DOUBLE strength,  // minimum thickness of a slab
    DOUBLE clearance, // minimum distance between slabs
    bool supported = true // true = solid, false = floating
  );

  // the following methods are used for the eroded style:
  
  // create slab with top at mean of highest-weight floors:
  void calculateMeanFloor(const std::vector< Sample >& samples);
  // add a copy of the first slab of h:
  void copySlab(const HexCell& h); 
  // replace all slabs by a single slab:
  void calculateMeanFloor();
};

class MeshVertex;

class Mesh
// A class to generate COLLADA geometries (meshes) of a specified
// material. To add a polygon, add vertices with addVertex, then
// call endPolygon(). Meshes are output and cleared when clear()
// is called, or when the destructor is called. Blender refuses to
// load meshes with huge coordinate arrays, so before adding a 
// polygon, use full() to check whether the mesh is getting too 
// large (and if so, call clear()). (This approach is sufficient 
// as long as polygons are small). The material names for all meshes
// that are output are collected in an array which can be accessed
// with materials(); the COLLADA nodes library should create a node
// for each mesh and bind the material "SOLID" to the material name
// in materials(); (the main function does this.)
{
    private:
  string material;
  std::vector< ElevatedPoint > coordinates;
  std::vector< int > polygonSizes;
  std::vector< int > polygons;
  int nrPending; // number of vertices in unfinished polygon
  static std::vector< string > theMaterials;
    public:
  Mesh(const string& _material = "");
  void setMaterial(const string& _material);
  static const std::vector< string >& materials();

  // for each vertex, a normal index should be specified:
  // 0-5: pointing to neighbours in directions 0-5
  // 6:   pointing upwards
  // 7:   pointing downwards 
  void addVertex(MeshVertex& v, int normalIndex);
  void addVertex(const Point& p, DOUBLE h, int normalIndex);
  void addVertex(DOUBLE x, DOUBLE y, DOUBLE h, int normalIndex);

  void endPolygon();
  bool full() const;
  void clear();
  ~Mesh();
};

class MeshVertex
// A class for re-usable vertices in a mesh. The class keeps 
// track of an index in the coordinates array of the mesh. Only 
// the first time the vertex is used, the coordinates are added 
// to the coordinates array. Note: MeshVertex objects become 
// invalid if passed to different Mesh objects, or whenever the 
// Mesh is cleared.
{
    private:
  ElevatedPoint point;
  int index;
    public:
  MeshVertex(Point p = Point(), DOUBLE h = 0);
  friend class Mesh;
};

class HexGrid
{
    private:
  ////////////////////////////
  // dimensions in grid space:

  int gridWidth; 
  int gridIndexRadius; // (gridWidth-1)/2
  int gridSize;        // gridWidth * gridWidth
  int indexDistanceToNeighbour[6]; 

  // the upper reference plane lies in the grid in the halfplane 
  // row - column >= discreteRidge. To omit the upper plane,  
  // discreteRidge is set to -3 * gridIndexRadius:
  int discreteRidge; 

  /////////////////////////////
  // dimensions in model space:

  DOUBLE cellSpacing;   // width of a cell, between parallel sides
  Point cellShape[6];   // corners of a grid cell, relative to the centre

  // camera position *before* rotation according to azimuth and 
  // altitude settings (rotation is done in the COLLADA output, so
  // that the viewing direction rotates with it):
  ElevatedPoint camera; 

  //////////////////////////////
  // dimensions in sample space:

  // number of cells in a grid row/column, per unit length in sample space
  int cellsPerUnit; 

  /////////////////
  // other members:

  // grid with slabs generated from points given to the constructor
  HexCell* grid;

  // calculate index in grid from row and column coordinates:
  int gridIndex(const int r, const int c) const;

  // calculate index in grid from coordinates in sample sapce:
  int gridIndex(const Point& p) const;

  // calculate coordinates in model space from index in grid
  // (so: *not* the inverse of the above!)
  Point location(int gridIndex) const;

    public:
  HexGrid( // all coordinates/dimensions given in sample space:
    const Options& options,
    const DOUBLE edgeLength, // edge length of bounding diamond
    const DOUBLE ridgeLine, // upper ref. plane at x > ridgeLine
    const std::vector< AnnotatedPoint >& pp  // sample points
  );

  // create geometries (meshes) for the COLLADA geometry library:
  void draw(int nrColours, DOUBLE parapetHeight = 0);

  // get the camera position (before azimuth/altitude rotation):
  ElevatedPoint cameraPosition() const;

  ~HexGrid();
};


/*********************************************************
* main header                                            *
*********************************************************/

int main(const int argc, const char* argv[]);


/*********************************************************
**********************************************************
**                                                      **
** IMPLEMENTATION                                       **
**                                                      **
**********************************************************
*********************************************************/


/*********************************************************
* Input/Output                                           *
*********************************************************/

const char* Complaint::errorLabel[5] = 
{ 
  "WARNING", 
  "Error",
  "Error while parsing IFS", 
  "Error", 
  "Internal error" 
};

Complaint::Complaint(ErrorLevel _errorLevel): 
  errorLevel(_errorLevel) 
{
};
  
Complaint::Complaint(
  ErrorLevel _errorLevel, 
  const string& expected, 
  const string& read
):
  errorLevel(_errorLevel)
{
  complaint << "expected \"" << expected << "\", but read \"" << read << "\".";
}

Complaint::~Complaint()
{
  cerr << errorLabel[(int) errorLevel] << ": " << complaint.str() << endl;
  if (errorLevel == usageError) 
    cerr << "Run pftrail without arguments for usage." << endl;
  if (errorLevel > warning) exit((int) errorLevel); 
}

template <typename T> Complaint& Complaint::operator<<(const T& t)
{
  complaint << t; return *this;
}

void CommentedInput::fillBuffer()
{
  string text;
  while (text.empty())
  {
    std::getline(source, text);
    if (source.fail()) return;
    // remove long comments
    if (text.substr(0,2) == "/*")
    {
      text.erase(0,2);
      while (text.find("*/") > text.length())
      {
        std::getline(source, text);
        if (source.fail()) 
        {
          Complaint(warning) << "unclosed comment in file";
          return;
        }
      }
      text.erase(0, text.find("*/") + 2);
    }
    // remove short comments:
    text.erase(std::min(text.find("//"), text.length()));
  }
  buffer.clear();
  buffer.str(text);
}

CommentedInput::CommentedInput(string& fileName):
  source(fileName.c_str())
{};

bool CommentedInput::fail() const
{
  return buffer.fail() || source.fail();
}

template < class T > CommentedInput& CommentedInput::operator>>(T& t)
{
  if (source.fail()) return *this;
  buffer >> t;
  while (buffer.fail() && buffer.eof() && !source.fail()) 
  {
    fillBuffer();
    buffer >> t;
  }
  return *this;
}

void CommentedInput::getLine(string& s)
{
  if (source.fail()) return;
  std::getline(buffer >> std::ws, s);
  if (s.empty()) 
  {
    fillBuffer();
    std::getline(buffer >> std::ws, s);
  }
}

void CommentedInput::getLowerCaseWord(string& s)
{
  *this >> s; makeLowerCase(s);
}

bool CommentedInput::find(const char* keyword, const string& name)
{
  string lowerCaseKeyword(keyword); makeLowerCase(lowerCaseKeyword);
  string lowerCaseName(name); makeLowerCase(lowerCaseName);
  string inputWord;
  while (!fail() && inputWord != lowerCaseName)
  {
    while (!fail() && inputWord != lowerCaseKeyword) 
      getLowerCaseWord(inputWord);
    if (fail()) return false;
    getLowerCaseWord(inputWord);
  }
  return (!fail() && inputWord == lowerCaseName);
}

void makeLowerCase(string& s)
{
  for (int i = 0; i < s.length(); ++i) s[i] = tolower(s[i]);
}

DOUBLE strtofstop(const string& str)
{
  stringstream s(str);
  DOUBLE x;
  s >> x >> std::ws;
  if (!s.eof()) Complaint(userError)
    << "expected a floating-point number but read \"" << str << "\".";
  return x; 
}

int strtoistop(const string& str)
{
  stringstream s(str);
  int n;
  s >> n >> std::ws;
  if (!s.eof()) Complaint(userError)
    << "expected an integer but read \"" << str << "\".";
  return n; 
}


/*********************************************************
* Geometry                                               *
*********************************************************/

Point::Point(DOUBLE _x, DOUBLE _y): 
  x(_x), y(_y) 
{};

Point Point::operator+(const Point& p) const 
{
  return Point(x + p.x, y + p.y);
}

Point& Point::operator+=(const Point& p) 
{
  x += p.x; y += p.y; return *this;
}

Point Point::operator-(const Point& p) const 
{
  return Point(x - p.x, y - p.y);
}

Point& Point::operator-=(const Point& p) 
{
  x -= p.x; y -= p.y; return *this;
}

Point Point::operator-() const 
{
  return Point(-x, -y);
}
  
DOUBLE Point::operator*(const Point& p) const 
{
  return x * p.x + y * p.y;
}

Point Point::operator*(const DOUBLE d) const 
{
  return Point(x * d, y * d);
}
 
Point Point::operator/(const DOUBLE d) const 
{
  return Point(x / d, y / d);
}
 
DOUBLE Point::sqLength() const 
{
  return x * x + y * y;
}

DOUBLE Point::length() const 
{
  return sqrt(sqLength());
}

ostream& operator<<(ostream& s, const Point& p)
{
  return s << '(' << p.x << ',' << p.y << ')';
}

ElevatedPoint::ElevatedPoint(): 
  Point(), h(0)  
{};

ElevatedPoint::ElevatedPoint(const Point& _p, DOUBLE _h): 
  Point(_p), h(_h)
{};

ElevatedPoint::ElevatedPoint(DOUBLE _x, DOUBLE _y, DOUBLE _h):
  Point(_x, _y), h(_h)
{};

ElevatedPoint ElevatedPoint::operator+(const ElevatedPoint& e) const 
{
  return ElevatedPoint(x + e.x, y + e.y, h + e.h);
}

ElevatedPoint& ElevatedPoint::operator+=(const ElevatedPoint& e) 
{
  x += e.x; y += e.y; h += e.h; return *this;
}

ElevatedPoint& ElevatedPoint::operator-=(const Point& p) 
{
  x -= p.x; y -= p.y; return *this;
}

ElevatedPoint& ElevatedPoint::operator-=(const ElevatedPoint& e) 
{
  x -= e.x; y -= e.y; h -= e.h; return *this;
}

ostream& operator<<(ostream& s, const ElevatedPoint& p)
{
  return s << '(' << p.x << ',' << p.y << ',' << p.h << ')';
}

AnnotatedPoint::AnnotatedPoint():
  ElevatedPoint(), a(sampleFill), range(-1,-1)
{};

AnnotatedPoint::AnnotatedPoint(
  const ElevatedPoint& _e, SampleCategory _a, DOUBLE r
):
  ElevatedPoint(_e), a(_a), range(_e.h - r, _e.h + r)
{};

AnnotatedPoint::AnnotatedPoint(
  const Point& _p, DOUBLE _h, SampleCategory _a, DOUBLE r
):
  ElevatedPoint(_p, _h), a(_a), range(h - r, h + r)
{};

// Below are helper functions for the computation of smallest 
// enclosing circles, computing the centre and squared radius 
// of a minimum enclosing circle for a set of points with 
// 1, 2 or 3 points on the boundary given.

// three boundary points given:
void computeMinDisk(
  Point p1, Point p2, Point p3, 
  Point& centre, DOUBLE& sqRadius
)
{
  // check if a pair determines the diameter:
  for (int i = 0; i < 3; ++i)
  {
    centre = (p1+p2)/2;
    sqRadius = (p2-p1).sqLength()/4 + practicallyZero;
    if ((p3-centre).sqLength() <= sqRadius) return;
    std::swap(p1,p2); std::swap(p2,p3);
  }
  // centre of bisectors satisfies:
  // (p2-p1)*2*centre = ||p2||^2 - ||p1||^2
  DOUBLE c[3][3]; // coefficients of the three equations
  c[0][0] = (p2.x-p1.x)*2; c[0][1] = (p2.y-p1.y)*2; c[0][2] = p2.sqLength()-p1.sqLength(); 
  c[1][0] = (p3.x-p2.x)*2; c[1][1] = (p3.y-p2.y)*2; c[1][2] = p3.sqLength()-p2.sqLength(); 
  c[2][0] = (p1.x-p3.x)*2; c[2][1] = (p1.y-p3.y)*2; c[2][2] = p1.sqLength()-p3.sqLength(); 
  // solve:
  if (abs(c[2][0]) > abs(c[1][0])) std::swap(c[1], c[2]);
  if (abs(c[1][0]) > abs(c[0][0])) std::swap(c[0], c[1]);
  c[0][1] /= c[0][0]; c[0][2] /= c[0][0]; c[0][0] = 1;
  for (int i = 1; i < 3; ++i)
    for (int j = 1; j < 3; ++j)
      c[i][j] -= c[i][0] * c[0][j];
  // 1st column is henceforth to be ignored, should be regarded as 1,0,0
  if (abs(c[2][1]) > abs(c[1][1])) std::swap(c[1], c[2]);
  c[1][2] /= c[1][1];
  c[0][2] -= c[0][1] * c[1][2];
  // 2nd column is henceforth to be ignored, should be regarded as 0,1,0
  centre = Point(c[0][2], c[1][2]);
  sqRadius = (p1-centre).sqLength();
}

// two boundary points given:
void computeMinDisk(std::vector< Point >& pp, 
  const Point& p1, const Point& p2,
  Point& centre, DOUBLE& sqRadius
)
{
  computeMinDisk(pp[0], p1, p2, centre, sqRadius);
  for (int i = 1; i < pp.size(); ++i)
    if ((pp[i]-centre).sqLength() > sqRadius) 
      computeMinDisk(p1, p2, pp[i], centre, sqRadius);
}

// one boundary point given:
void computeMinDisk(std::vector< Point >& pp, 
  const Point& p1,
  Point& centre, DOUBLE& sqRadius
)
{
  computeMinDisk(pp[0], pp[1], p1, centre, sqRadius);
  for (int i = 2; i < pp.size(); ++i)
    if ((pp[i]-centre).sqLength() > sqRadius) 
      computeMinDisk(pp, p1, pp[i], centre, sqRadius);
}

DOUBLE radius(const std::vector< AnnotatedPoint >& points)
{
  if (points.empty()) return 0;
  std::vector < Point > pp;
  for (
    std::vector< AnnotatedPoint >::const_iterator p = points.begin(); 
    p != points.end(); 
    ++p
  )
    pp.push_back(*p);

  // ensure minimum input size of three points:
  while (pp.size() < 3) pp.push_back(points[0]); 

  // permute randomly:
  for (int i = 0; i < pp.size(); ++i) 
    std::swap(pp[i], pp[i + (rand() % (pp.size()-i))]);

  // solve using incremental construction:
  Point centre; DOUBLE sqRadius;
  computeMinDisk(pp[0], pp[1], pp[2], centre, sqRadius);
  for (int i = 3; i < pp.size(); ++i)
    if ((pp[i]-centre).sqLength() > sqRadius) 
      computeMinDisk(pp, pp[i], centre, sqRadius);
  return sqrt(sqRadius);
}

Rotation::Rotation(const Point& _x, const Point& _y): 
  x(_x), y(_y) 
{};

Rotation Rotation::transpose() const 
{
  return Rotation(Point(x.x, y.x), Point(x.y, y.y));
}

Rotation Rotation::operator*(const Rotation& r) const 
{
  Rotation t = r.transpose();
  return Rotation(
    Point( x * t.x, x * t.y ),
    Point( y * t.x, y * t.y ) 
  );
}

Rotation Rotation::operator*=(const Rotation& r) 
{
  Rotation t = r.transpose();
  x = Point( x * t.x, x * t.y );
  y = Point( y * t.x, y * t.y );
  return *this;
}

Point Rotation::operator*(const Point& p) const 
{
  return Point(x * p, y * p);
}

DOUBLE Rotation::sqScale() const
{
  return x.x * x.x + y.x * y.x;
}

DOUBLE Rotation::scale() const
{
  return sqrt(sqScale()); 
}

Similarity::Similarity(
  const Rotation& _rotation, 
  const ElevatedPoint& _translation, 
  DOUBLE _scale
):
  rotation(_rotation), translation(_translation), scale(_scale) 
{};

Similarity Similarity::operator*(const Similarity& s) const 
{
  return Similarity(
    rotation * s.rotation,
    operator()(s.translation),
    scale * s.scale
  );
}

Similarity& Similarity::operator*=(const Similarity& s) 
{
  translation = operator()(s.translation);
  rotation *= s.rotation;
  scale *= s.scale;
  return *this;
}

Similarity& Similarity::apply(const Similarity& s) 
{
  translation = s(translation);
  rotation = s.rotation * rotation;
  scale *= s.scale;
  return *this;
}

Point Similarity::operator()(const Point& p) const 
{
  return rotation * p + (Point) translation;
}

ElevatedPoint Similarity::operator()(const ElevatedPoint& e) const 
{
  return ElevatedPoint(rotation * (Point) e, scale * e.h) + translation;
}

void Similarity::print(
  std::ostream& o, const string& leftLabel, const string& rightLabel
) const
{
  o << leftLabel << std::fixed << std::setprecision(3)
    << "| " << std::setw(6) << rotation.x.x << ' '
            << std::setw(6) << rotation.x.y << ' '
            << "  0   " << ' '
            << std::setw(6) << translation.x << " |" << rightLabel << endl;
  o << std::setw(leftLabel.length() + 2)
    << "| " << std::setw(6) << rotation.y.x << ' '
            << std::setw(6) << rotation.y.y << ' '
            << "  0   " << ' '
            << std::setw(6) << translation.y << " |" << endl;
  o << std::setw(leftLabel.length() + 2) 
    << "| " << "  0   " << ' '
            << "  0   " << ' '
            << std::setw(6) << scale << ' '
            << std::setw(6) << translation.h << " |" << endl;
  o << std::setw(leftLabel.length() + 2) 
    << "| " << "  0   " << ' '
            << "  0   " << ' '
            << "  0   " << ' '
            << "  1   " << " |" << endl;
  o.unsetf(std::ios_base::floatfield);
}

template < class T > 
const Rotation TrianglePoint< T >::TriToCar = 
  Rotation(Point(1.0, 0.5), Point(0, sqrt((DOUBLE) 0.75)));

template < class T >
TrianglePoint< T >::TrianglePoint(int s, T a, T b) 
{
  x[s%3] = a; x[(s+1)%3] = b; x[(s+2)%3] = -(a+b);
}

template < class T >
TrianglePoint< T >::TrianglePoint(const Point& p) 
{
  x[0] =  p.x - SQRTTHIRD * p.y;
  x[1] = -p.x - SQRTTHIRD * p.y;
  x[2] = -(x[0] + x[1]); 
};

template < class T >
T TrianglePoint< T >::operator[](int i) const 
{ 
  return x[i % 3]; 
}

template < class T >
TrianglePoint<T>::operator Point() const 
{
  return TriToCar * Point(x[0], x[2]);
}

TrianglePoint< int > round(const TrianglePoint< DOUBLE >& t) 
{
  DOUBLE roundingError[3];
  for (int i = 0; i < 3; ++i) roundingError[i] = round(t[i]) - t[i];

  // first find coordinate with smallest difference to rounded number
  int closest = 0;
  for (int i = 1; i < 3; ++i) {
    if (abs(roundingError[i]) >= abs(roundingError[closest])) continue;
    closest = i;
  }
  // round this coordinate; distribute the rounding error over the other two
  // and round them accordingly

  TrianglePoint< int > r(closest,
    (int) round(t[closest]),
    (int) round(t[closest+1] - 0.5*roundingError[closest])
  );

  return r;
}

void findMinimumBoundingDiamond(
  const std::vector< AnnotatedPoint >& pp,
  Point& centre,
  DOUBLE& edgeLength
)
{
  if (pp.size() == 0) { edgeLength = 0; return; }
  DOUBLE minCoord[2];
  DOUBLE maxCoord[2];
  TrianglePoint< DOUBLE > t(pp[0]);
  minCoord[0] = t[0]; maxCoord[0] = t[0]; 
  minCoord[1] = t[1]; maxCoord[1] = t[1]; 
  for (int i = 1; i < pp.size(); ++i) 
  {
    TrianglePoint< DOUBLE > t(pp[i]);
    minCoord[0] = min(minCoord[0], t[0]); maxCoord[0] = max(maxCoord[0], t[0]); 
    minCoord[1] = min(minCoord[1], t[1]); maxCoord[1] = max(maxCoord[1], t[1]); 
  }
  centre = (Point) TrianglePoint< DOUBLE >(0,
    0.5*(minCoord[0] + maxCoord[0]),
    0.5*(minCoord[1] + maxCoord[1])
  );
  edgeLength = max(
    maxCoord[0] - minCoord[0],
    maxCoord[1] - minCoord[1]
  );
}


/*********************************************************
* IFS (parsing, analysis and sampling)                   *
*********************************************************/

Segment::Segment(Similarity _transformation, int _generator):
  transformation(_transformation),
  generator(_generator)
{};

void IFS::skip(CommentedInput& text, const char* s) const
{
  std::string aword;
  text.getLowerCaseWord(aword);
  if (text.fail() || aword != s) Complaint(ifsParsingError, s, aword);
}

void IFS::skip(CommentedInput& text, const char c) const
{
  char acharacter;
  text >> acharacter; acharacter = tolower(acharacter);
  if (text.fail() || acharacter != c) 
    Complaint(ifsParsingError, string() + c, string() + acharacter);
}

string IFS::generatorName(int g) const
{
  if (g == -1) return string("zero/bridge generator");
  if (g == 0)  return string("principal generator");
  else         return string("generator ") + (char) ('A' + g);
}

void IFS::writeAreaEquation(ostream& s, const Equation& e)
{
  bool zero = true;
  for (int i = 0; i < e.size()-1; ++i)
  {
    if (abs(e[i]) < practicallyZero) continue;
    if (zero)
    {
      if (e[i] < 0) s << '-';
    }
    else s << (e[i] < 0 ? " - " : " + ");
    zero = false;
    DOUBLE coefficient = abs(e[i]);
    if (coefficient != 1) s << coefficient << ' '; 
    s << "f(" << (char) ('A' + i) << ')';
  }
  if (zero) s << '0';
  s << " = " << std::setprecision(10) << e.back();
}

void IFS::writeAreaEquations(ostream& s, const std::vector< Equation >& ee)
{
  for (std::vector< Equation >::const_iterator e = ee.begin(); e != ee.end(); ++e)
  {
    writeAreaEquation(s, *e); 
    s << endl;
  }
}

void IFS::parse(CommentedInput& text)
{
  definition.clear();
  std::string aword;

  // skip any (alternative) names of the IFS
  while (!text.fail()) 
  {
    text.getLowerCaseWord(aword);
    if (aword == "ifs") text >> aword; else break;
  }

  // determine coordinate system
  Rotation grid;
  CoordinateSystem coordinateSystem = planarCoordinates;
  if (aword == "square")
    grid = Rotation(Point(1, 0), Point(0, 1));
  else if (aword == "triangular")
    grid = Rotation(Point(1, 0.5), Point(0, SQRTTHREE * 0.5));
  else if (aword == "cubic" || aword == "barycentric")
  {
    grid = Rotation(Point(1.5, 0), Point(-SQRTTHREE * 0.5, -SQRTTHREE));
    coordinateSystem = cubicCoordinates;
  }    
  else 
    Complaint(ifsParsingError, "square\", \"triangular\", or \"cubic", aword);
  skip(text, "grid");

  // the definition is recorded in Cartesian coordinates, 
  // regardless of the input format
  definition << "square grid" << endl;

  // determine number of generators
  int n; text >> n;
  if (text.fail()) Complaint(ifsParsingError) 
    << "could not read number of generators or segments.";
  int nrGenerators;
  bool multipleGeneratorSyntax;
  text.getLowerCaseWord(aword);
  if (aword == "segments") {
    multipleGeneratorSyntax = false;
    nrGenerators = 1;
    definition << n << ' ' << "segments" << endl;
  }
  else if (aword == "generators") {
    multipleGeneratorSyntax = true;
    nrGenerators = n;
    definition << n << ' ' << "generators" << endl;
  }
  else
    Complaint(ifsParsingError, "segments\" or \"generators", aword);
  
  // prepare for collecting area equations and continuity status 
  continuous = true;
  std::vector< Equation > areaEquations;
  if (Options::verbosity >= 2)
  {
    cerr << "f(x):= area filled divided by squared distance between endpoints of generator x." << endl
         << "The following equations are derived from the generators:" << endl;
  }

  // read generators one by one
  for (int i = 0; i < nrGenerators; ++i)
  {
    bool areaGiven = false;
    DOUBLE givenArea;
    if (multipleGeneratorSyntax)
    {
      skip(text, "generator");
      char generatorLabel = (char) ('a' + i);
      skip(text, generatorLabel); 
      skip(text, ':');
      definition << endl << "generator " << generatorLabel << ':' << endl;

      text.getLowerCaseWord(aword);
      if (aword == "fills")
      {
        text >> givenArea;
        if (text.fail()) Complaint(ifsParsingError) 
          << "could not read area filled by " << generatorName(i) << ".";
        if (givenArea < 0) Complaint(ifsParsingError)
          << "negative area specified for " << generatorName(i) << ".";
        areaGiven = true; 
        definition << "fills " << givenArea << endl;
        text.getLowerCaseWord(aword);
      }

      n = strtoistop(aword);
      if (text.fail() || n == 0) Complaint(ifsParsingError)
        << "could not read number of segments of " << generatorName(i) << ".";
      skip(text, "segments");
      definition << n << " segments" << endl;
    }

    if (n < 2) Complaint(ifsParsingError)
      << "less than two segments specified for " << generatorName(i) << ".";
    skip(text, "segment"); skip(text, "values:");
    definition << endl << "segment values:" << endl;

    // initialize the generator
    generators.push_back(Generator());
    Point location(0,0);
    areaEquations.push_back(Equation(nrGenerators+1));

    // read the segments one by one
    for (int j = 1; j <= n; ++j)
    {
      int segmentLabel; text >> segmentLabel;
      if (text.fail()) Complaint(ifsParsingError) 
        << "expected segment index " << j << ".";
      if (segmentLabel != j) Complaint(ifsParsingError) 
        << "expected segment index " << j << ", but read " << segmentLabel << ".";
      skip(text, ':');

      DOUBLE x, y; int d, r; int g = nrGenerators-1;
      string segmentParameters;
      text.getLine(segmentParameters);
      stringstream segmentParameterStream(segmentParameters);
      readSegmentParameters(segmentParameterStream, x, y, g, d, r, coordinateSystem);
      string aword;
      segmentParameterStream >> aword;
      if (!aword.empty()) Complaint(userError) << "unexpected text \"" << aword << "\" on segment line.";
      Point edge = grid * Point(x, y);
      if (abs(edge.x) < practicallyZero) edge.x = 0;
      if (abs(edge.y) < practicallyZero) edge.y = 0;
 
      definition << j << ":  " << edge.x << ", " << edge.y << ", ";
      if (nrGenerators > 1) 
      {
        if (g == -1) definition << 'X' << endl;
        else definition << (char) ('A' + g) << ", " << d << ", " << r << endl;
      }
      else definition << d << ", " << r << endl;

      if (d == 0) continuous = false;
      Point anchor = location;
      location += edge;
      if (d == -1) { edge = -edge; anchor = location; }
      generators.back().segments.push_back(
        Segment(
          Similarity(
            Rotation(
              Point(edge.x, -r * d * edge.y), 
              Point(edge.y,  r * d * edge.x)
            ),
            ElevatedPoint(anchor, 0),
            d
          ),
          g
        )
      );
      if (d != 0) areaEquations.back()[g] -= edge.sqLength();
    }

    // verify that the generator has non-zero length
    DOUBLE generatorSize = location.sqLength();
    if (generatorSize == 0) Complaint(userError) 
      << "distance between end points of " << generatorName(i) << " is zero"; 

    // finish the area equations
    areaEquations.back()[i] += location.sqLength();
    if (Options::verbosity >= 2)
    { 
      cerr << "- from " << generatorName(i) << ": ";
      writeAreaEquation(cerr, areaEquations.back());
      cerr << endl;
    }
    if (areaGiven) 
    {
      areaEquations.push_back(Equation(nrGenerators+1));
      areaEquations.back()[i] = generatorSize;
      areaEquations.back()[nrGenerators] = givenArea;
      if (Options::verbosity >= 2) 
      {
        cerr << "- specified area filled by " << generatorName(i) << ": ";
        writeAreaEquation(cerr, areaEquations.back());
        cerr << endl;
      }
    }

    // rotate and scale such that endpoints are (0,0) and (1,0)
    Similarity normalization(
      Rotation(
        Point( location.x/generatorSize, location.y/generatorSize),
        Point(-location.y/generatorSize, location.x/generatorSize)
      ),
      ElevatedPoint(0, 0, 0),
      1
    );
    for (
      std::vector< Segment >::iterator s = generators[i].segments.begin(); 
      s != generators[i].segments.end(); 
      ++s
    )
      s->transformation.apply(normalization);
  }

  if (areaEquations.size() == nrGenerators) 
  {
    // no area was explicitly specified, w.l.o.g. assume f(A) = 1
    areaEquations.push_back(Equation(nrGenerators+1));
    areaEquations.back()[0] = 1;
    areaEquations.back()[nrGenerators] = 1;
    if (Options::verbosity >= 2)
    { 
      cerr << "- for a suitable choice of area unit we may assume: ";
      writeAreaEquation(cerr, areaEquations.back());
      cerr << endl;
    }
  }

  // solve the area equations; implementation is suboptimal but "n" is really small
  for (int j = nrGenerators-1; j >= 0; --j)
  {
    // find an equation i with coefficients j+1...n-1 all zero, 
    // and among those, the one with largest absolute value for coefficient j
    int bestEquation = -1; DOUBLE divisor = 0; 
    for (int i = 0; i < areaEquations.size(); ++i)
    {
      if (abs(areaEquations[i][j]) <= abs(divisor)) continue;
      bool eligible = true; 
      for (int k = j+1; eligible && k < nrGenerators; ++k)
        eligible = (areaEquations[i][k] == 0);
      if (!eligible) continue;
      bestEquation = i;
      divisor = areaEquations[i][j]; 
    }
    if (abs(divisor) <= practicallyZero)
    {
      if (Options::verbosity > 0) 
      {
        cerr << "Area equations reduce to:" << endl;
        writeAreaEquations(cerr, areaEquations);
      }
      Complaint(userError) << "cannot determine area filled by " << generatorName(j) << ".";
    }
    if (bestEquation != j) swap(areaEquations[bestEquation], areaEquations[j]);

    // normalize:
    for (int k = 0; k <= nrGenerators; ++k) areaEquations[j][k] /= divisor;

    // eliminate generator j from all other rows
    for (int i = 0; i < areaEquations.size(); ++i)
    { 
      if (i == j) continue;
      DOUBLE multiplier = areaEquations[i][j];
      for (int k = 0; k <= nrGenerators; ++k) 
        areaEquations[i][k] -= multiplier * areaEquations[j][k];
    }
  }
  // now, the coefficients matrix should be in the form
  // |   1   0   ...   0   f(A)   |
  // |   0   1   ...   0   f(B)   |
  // |   :   :    :    :    :     |
  // |   0   0   ...   1   f(N)   |
  // |   0   0   ...   0    0     |
  // |   :   :    :    :    :     |

  // check for inconsistencies:
  for (int i = nrGenerators; i < areaEquations.size(); ++i)
    if (abs(areaEquations[i].back()) > practicallyZero) 
    {
      if (Options::verbosity > 0) 
      {
        cerr << "Area equations reduce to:" << endl;
        writeAreaEquations(cerr, areaEquations);
      }
      Complaint(userError) << " specified areas are inconsistent.";
    }
  
  for (int i = 0; i < nrGenerators; ++i)
    if (abs(areaEquations[i].back()) <= practicallyZero) 
    {
      if (Options::verbosity > 0) 
      {
        cerr << "Area equations reduce to:" << endl;
        writeAreaEquations(cerr, areaEquations);
      }
      Complaint(userError) 
        << generatorName(i) << " does not generate a plane-filling curve.";
    }

  if (Options::verbosity >= 2 && nrGenerators >= 2)
  {
    cerr << "(f(A), ..., f(" << (char) ('A' + nrGenerators-1) << ")) = (";
    cerr << std::setprecision(3) << areaEquations[0].back();
    for (int i = 1; i < nrGenerators; ++i)
      cerr << ", " << std::setprecision(3) << areaEquations[i].back();
    cerr << ')' << endl;
  }

  // fill in the third rows of the transformation matrices
  for (int i = 0; i < nrGenerators; ++i)
  {
    DOUBLE t = 0;
    for (
      std::vector< Segment >::iterator s = generators[i].segments.begin(); 
      s != generators[i].segments.end(); 
      ++s
    )
    {
      if (s->transformation.scale != 0)
        s->transformation.scale *= 
          s->transformation.rotation.sqScale() * 
          areaEquations[s->generator].back() / areaEquations[i].back();
      if (s->transformation.scale > 0)
      {
        s->transformation.translation.h = t;
        t += s->transformation.scale;
      }
      else
      {
        t -= s->transformation.scale;
        s->transformation.translation.h = t;
      }
    }
  }

  if (Options::verbosity >= 2) 
  {
    cerr << "The following transformation matrices have been constructed:" << endl;
    for (int i = 0; i < nrGenerators; ++i) {
      if (nrGenerators > 1) cerr << "Generator " << (char) ('A' + i) << ':' << endl;
      for (
        std::vector< Segment >::const_iterator s = generators[i].segments.begin(); 
        s != generators[i].segments.end(); 
        ++s
      )
      {
        string rightLabel;
        if (nrGenerators > 1 && s->generator >= 0) 
          rightLabel = string(" ") + (char) ('A' + s->generator); 
        s->transformation.print(cerr, "+ ", rightLabel);
      }
    }
  }

  // calculate lower bound on diameter and, for each generator, upper bound 
  // on "radius" (max. distance from curve point to closest endpoint) 
  // (these numbers are needed to calculate the sampling density) 
  if (Options::verbosity >= 2)
    cerr << "Generating sketches for all generators to estimate shape parameters..." << endl;

  // apply settings for rough sketches  
  const DOUBLE maxEdgeLength = 0.05;
  applySketchSettings(maxEdgeLength); 

  DOUBLE overallMaxRadius = 0;
  for (int i = 0; i < nrGenerators; ++i) {
    Generator& g = generators[i];
    std::vector< AnnotatedPoint > pp;
    bool skipJump = true;
    generatePoints(pp, i, Similarity(), skipJump);
    if (i == 0) diameter = radius(pp) * 2.0;
    // find maximum squared distance from any sample point to closest end point
    g.maxRadius = 0; 
    for (std::vector< AnnotatedPoint >::const_iterator p = pp.begin(); p != pp.end(); ++p) 
      g.maxRadius = std::max(g.maxRadius,
        std::min(
          ((Point) *p - Point(0,0)).sqLength(), 
          ((Point) *p - Point(1,0)).sqLength())
      );
    // calculate actual distance
    g.maxRadius = sqrt(g.maxRadius);
    overallMaxRadius = std::max(overallMaxRadius, g.maxRadius);
  }
  // calculate maximum further expansion in recursion
  DOUBLE expansion = overallMaxRadius * maxEdgeLength / (1.0 - maxEdgeLength);
  for (int i = 0; i < nrGenerators; ++i)
    generators[i].maxRadius += expansion;

  if (Options::verbosity >= 2)
  {
    cerr << "Diameter divided by the distance between endpoints is at least " << diameter << "." << endl;
    cerr << "For any point on the curve";
    if (nrGenerators > 1) cerr << "s generated by generators A..." << (char) ('A' + nrGenerators-1);
    cerr << ", the distance to the closest endpoint divided by the distance between endpoints is at most";
    if (nrGenerators == 1)
      cerr << ' ' << generators[0].maxRadius << '.' << endl;
    else 
    { 
      cerr << ", respectively" << endl << generators[0].maxRadius;
      for (int i = 1; i < nrGenerators; ++i) cerr << ", " << generators[i].maxRadius;
      cerr << '.' << endl;
    }
  }
}

bool IFS::isContinuous() const
{
  return continuous;
}

DOUBLE IFS::diameterLowerBound() const
{
  return diameter;
}

ElevatedPoint IFS::operator()(std::string index) const
{
  if (index[0] - '0' == generators[0].segments.size()) return ElevatedPoint(Point(1,0), 1);
  Similarity s;
  int whichEnd = 0;
  int g = 0;
  for (int i = 0; i < index.length() && index[i] != ' ' && index[i] != '\t'; ++i) 
  {
    if (g == -1 || s.scale == 0) Complaint(userError) 
      << "IFS::operator() called on \"" << index.substr(0,i) << "...\", but \""
      << index.substr(0,i) << "\" denotes curve section with zero two-dimensional content.";
    const std::vector< Segment >& segments = generators[g].segments;
 
    int j = index[i] - '0';
    if (j < 0 || j >= segments.size()) Complaint(userError) 
      << "IFS::operator() called on \"" << index.substr(0,i+1) << "...\", but \""
      << index[i] << "\" does not denote a subsegment of \"" << index.substr(0,i) << "\".";

    if (s.scale < 0) j = segments.size() - 1 - j;
    if (segments[j].transformation.scale < 0) whichEnd = 1 - whichEnd;
    s *= segments[j].transformation;
    g = segments[j].generator;
  }
  return s(ElevatedPoint(Point(whichEnd,0), whichEnd));
}

ElevatedPoint IFS::zoom(const ElevatedPoint& p)
{
  Point pz; DOUBLE hz;

  // scale distance x to x^{1/a}
  DOUBLE oldDistance = p.length();
  DOUBLE newDistance = pow(oldDistance, 1.0/zoomExponent);
  if (newDistance == 0) pz = Point(0,0);
  else pz = ((Point) centre) + ((Point) p) * newDistance / oldDistance;

  // same for elevation, but with additional scaling (to stay within 0-1 range)
  // and with exponent 1/(2a-1) instead of exponent a (much stronger zoom),
  // to compensate for the fact that small parts of the curve are "flatter"
  DOUBLE oldVerticalDistance = abs(p.h);
  DOUBLE newVerticalDistance = pow(oldVerticalDistance, invElevationZoomExponent);
  if (newVerticalDistance == 0) hz = 0;
  else hz = invScaledElevRange * 
    ((p.h < 0 ? -newVerticalDistance : newVerticalDistance) - minScaledElev);

  return ElevatedPoint(pz, hz);
}
 
void IFS::buildBridge(
  std::vector< AnnotatedPoint >& pp, 
  Point a, Point b, const DOUBLE h
)
{
  // calculate length and width in steps of size (roughly) cell diameter * sqrt(3)/4 
  Point d = b - a;
  DOUBLE length = d.length();
  int nrSamples = (int) (length / bridgeSampleDistance) + 1;
  DOUBLE halfwidth = bridgeWidthFactor * pow(d.sqLength(), 1.0/6.0);
  int corew = (bridgeCoreWidth == 0 ? -1 : 
    bridgeCoreWidth / (bridgeSampleDistance * 2.0));

  // calculate single step size:
  Point lngStep = d / nrSamples;
  Point latStep = Point(lngStep.y, -lngStep.x); // side step 

  for (int i = 0; i <= nrSamples; ++i)
  {
    // calculate local width w (steps left/right of centre line)
    int w;
    switch ((int) bridgeStyle)
    {
      case pointedBridges: 
        w = (int) (
          pow(pow((DOUBLE) nrSamples, 4.0) - pow((DOUBLE) 2*i-nrSamples, 4.0), 0.25) 
          * halfwidth / nrSamples + 0.5
        ); 
        break;
      case rectangularBridges: 
        w = (int) (halfwidth + 0.5);
        break;
      default: 
        Complaint(internalError) << "unrecognized bridge style";
        w = -1;
    }
    Point c = a; c -= (latStep * w);
    for (int j = -w; j <= w; ++j)
    {
      pp.push_back(AnnotatedPoint(
        c, h, 
        ( -corew <= j && j <= corew ? sampleJumpCore : sampleJump ), 
        length
      )); 
      c += latStep;
    }
    a += lngStep;
  }
}

void IFS::generatePoints(
  std::vector< AnnotatedPoint >& pp, 
  int generator, 
  const Similarity& s,
  bool& jump
) 
{
  if (generator == -1) 
  {
    jump = true;
    return;
  }
  Generator& g = generators[generator];
  bool mustRefine = false;
  ElevatedPoint p[3];
  p[s.scale >= 0 ? 0 : 1] = s(ElevatedPoint(Point(0,0), 0));
  p[s.scale >= 0 ? 1 : 0] = s(ElevatedPoint(Point(1,0), 1));

  if (zoomExponent == 0)
    mustRefine = s.rotation.sqScale() > g.maxSqSampleLength;
  else 
  {
    // Points at distance x are moved to x^{1/z}. This appears to lead to the following:
    // Let X be the origin (centre of zoom), A be any point, and r a fixed radius.
    // Points at distance r from A that are furthest from A after zooming are found:
    // -  if r <= 2 ||XA||: on the circle around X with radius ||XA||;
    //    The distance after zooming is:  r ||XA||^{1/z-1}
    // -  if r >= 2 ||XA||: on the ray from A through X.
    //    The distance after zooming is:  ||XA||^{1/z} + (r-||XA||)^{1/z}
    // The effect of zooming is stronger near the origin, so what matters is the case:
    // A is the segment endpoint closest to X; r is the segment's maxRadius.
 
    Point A = p[0]; DOUBLE sqDistA = A.sqLength();
    Point B = p[1]; DOUBLE sqDistB = B.sqLength();
    if (sqDistA < sqDistB) { A = B; sqDistA = sqDistB; }
    // A is now the segment endpoint furthest from X. 
    DOUBLE distA = sqrt(sqDistA);

    DOUBLE r = g.maxRadius * s.rotation.scale();
    if (r <= distA * 2.0) 
    {
      // first a test that is sufficient if distA is really small, to avoid 
      // division by zero or overflow in the real test
      if (pow(distA, invZoomExponent) <= curveSampleRadius)
        mustRefine = false;
      else 
        mustRefine = (r * pow(distA, invZoomExponent-1.0) > curveSampleRadius);
    }
    else 
      mustRefine = 
        (pow(distA, invZoomExponent) + pow(r - distA, invZoomExponent) > curveSampleRadius);
  }
  if (!mustRefine)   
  {
    if (zoomExponent != 0) p[1] = zoom(p[1]);
    if (jump)
    {
      if (zoomExponent != 0) p[0] = zoom(p[0]);
      if (bridgeStyle != noBridges && !pp.empty()) buildBridge(pp, pp.back(), p[0], p[0].h);
      pp.push_back(AnnotatedPoint(p[0], sampleFill));
    }
    pp.push_back(AnnotatedPoint(p[1], sampleFill));
    jump = false;
  }
  else if (s.scale > 0)
    for (int i = 0; i < g.segments.size(); ++i)
      generatePoints(pp, g.segments[i].generator, s * g.segments[i].transformation, jump);
  else // (s.scale < 0)
    for (int i = g.segments.size()-1; i >= 0; --i)
      generatePoints(pp, g.segments[i].generator, s * g.segments[i].transformation, jump);
}

void IFS::applySketchSettings(DOUBLE sampleLength)
{
  centre = ElevatedPoint(0,0,0);
  zoomExponent = 0;
  invZoomExponent = 1;
  invElevationZoomExponent = 1;
  curveSampleRadius = sampleLength; // should not have any effect
  for (std::vector< Generator >::iterator g = generators.begin(); g != generators.end(); ++g)
  {
    g->maxSampleLength = sampleLength;
    g->maxSqSampleLength = pow(g->maxSampleLength, 2.0);
  }
  bridgeStyle = noBridges;
}

void IFS::applyUserSettings(const Options& options)
{
  zoomExponent = options.zoomExponent;
  centre = options.view.orientation(options.view.centre);
  invZoomExponent = (DOUBLE) 1.0 / max((DOUBLE) 1.0, zoomExponent);
  invElevationZoomExponent = (DOUBLE) 1.0 / elevationZoomExponent(zoomExponent);
  minScaledElev = -pow(centre.h, invElevationZoomExponent);
  DOUBLE maxScaledElev =  pow((DOUBLE) 1.0-centre.h, invElevationZoomExponent);
  invScaledElevRange = (DOUBLE) 1.0 / (maxScaledElev - minScaledElev);
  curveSampleRadius = 0.25 * pow(diameter, invZoomExponent) / options.resolution / options.accuracy;
  for (std::vector< Generator >::iterator g = generators.begin(); g != generators.end(); ++g)
  {
    g->maxSampleLength = curveSampleRadius / g->maxRadius;
    g->maxSqSampleLength = pow(g->maxSampleLength, 2.0);
  }
  bridgeStyle = options.bridgeStyle;
  bridgeSampleDistance = SQRTTHREE * 0.25 * diameter / options.resolution;
  bridgeWidthFactor = 0.5 * options.bridgeWidth * pow(diameter, 1.0/3.0) / bridgeSampleDistance;
  bridgeCoreWidth = options.bridgeCoreWidth;
}

void IFS::generatePoints(std::vector< AnnotatedPoint >& pp, const Similarity& s)
{
  bool jump = true;

  // when zooming in, to make the precision of DOUBLE last longest (especially w.r.t.
  // elevation), we shift the traversal so that the centre of zoom lies in the origin;
  Similarity shifted = s; 
  if (zoomExponent != 0) shifted.translation -= centre;

  generatePoints(pp, 0, shifted, jump);
  // the zoom function that is used by generatePoints undoes the shift so that the 
  // output is where it should be
}

void IFS::report(ostream& o)
{
  o << "IFS definition:" << endl << endl << definition.str() << endl;
}

void skipChar(istream& text, const char c)
// only used by readSegmentParameters
{
  char acharacter;
  text >> acharacter; acharacter = tolower(acharacter);
  if (text.fail() || acharacter != c) 
    Complaint(ifsParsingError, string() + c, string() + acharacter);
}

void readSegmentParameters(
  stringstream& text,
  DOUBLE& x,
  DOUBLE& y,
  int& g, // generator, must be preset to number of generators minus 1,
          // and will only be read if larger than 0
          // is set to -1 if the segment is a bridge
  int& d,
  int& r,
  CoordinateSystem c)
{
  bool readGenerator = (g > 0);

  text >> x;
  if (text.fail()) Complaint(ifsParsingError) << "could not read segment's x-coordinate.";

  skipChar(text, ',');
  text >> y;
  if (text.fail()) Complaint(ifsParsingError) << "could not read segment's y-coordinate.";

  if (c == cubicCoordinates)
  {
    DOUBLE z;
    skipChar(text, ',');
    text >> z;
    if (text.fail()) Complaint(ifsParsingError) << "could not read segment's z-coordinate.";
    z += x; z += y; z /= 3.0;
    x -= z; y -= z;
  }

  if (readGenerator)
  {
    skipChar(text, ',');
    char gchar;
    text >> gchar;
    if (text.fail()) Complaint(ifsParsingError) << "could not read segment's generator.";
    if (gchar == 'X') { 
      g = -1; d = 0; r = 0; return;
    }
    if (gchar < 'A' || gchar > 'A' + g) 
      Complaint(ifsParsingError) << "invalid generator identifier \"" << gchar << "\".";
    g = gchar - 'A';
  }
  else g = 0;

  skipChar(text, ',');
  text >> d;
  if (text.fail()) Complaint(ifsParsingError) << "could not read segment's direction indicator.";
  if (d*d > 1 || (readGenerator && d == 0)) 
    Complaint(ifsParsingError) << "invalid direction indicator (" << d << ").";

  skipChar(text, ',');
  text >> r;
  if (text.fail()) Complaint(ifsParsingError) << "could not read segment's reflection indicator.";
  if (r*r > d*d) 
    Complaint(ifsParsingError) << "invalid reflection indicator (" << d << ").";

  if (d == 0) g = -1;
}


/*********************************************************
* Options                                                *
*********************************************************/

// first declare some helpers:

const OptionSyntax* findOption(const OptionSyntax options[], string argv);

struct Argument 
// class to store options with arguments;
// different options must have unique ranks and ids
{
  int rank;
  string id;
  string value;
  Argument(int _rank, string _id, const string& _value);

  // functions to extract a number from the string value, and 
  // verify that the number is within given bounds:
  int intValue(const int min, const int max) const;
  DOUBLE floatValue(const DOUBLE min, const DOUBLE max) const;

  // comparison operator to sort by rank:
  bool operator<(const Argument& a) const;
};


// and here is the implementation:

OptionSyntax::OptionSyntax(int _rank, string _id, bool _requiresValue):
    rank(_rank), id(_id), requiresArgument(_requiresValue)
{
};

View::View() 
{
  reset();
}

void View::reset()
{
  orientation = Similarity();
  ridge = true;
  ridgeLine = 0.5;
  centre = ElevatedPoint();
  defaultCentre = true;
  altitude = 18.5;
  azimuth = 0;
  distance = 1.1;
}

void View::parse(CommentedInput& text, const IFS& ifs)
{
  reset();

  // collect view/rendering instructions
  std::vector< Argument > arguments;
  string inputWord;
  while (true)
  {
    text >> inputWord;
    if (text.fail()) break;
    const OptionSyntax* o = findOption(syntax, inputWord);
    if (!o) 
      Complaint(userError) << "unrecognized keyword \"" << inputWord << "\" in rendering/view configuration.";
    if (o->id == "ifs") break;
    if (o->id == "view" && !arguments.empty()) break;
    if (o->requiresArgument)
    {     
      text.getLine(inputWord);
      if (inputWord.empty()) 
        Complaint(userError) << "missing value for \"" << o->id << "\" in rendering/view configuration.";
    }
    if (o->id != "view") arguments.push_back(Argument(o->rank, o->id, inputWord));           
  }

  // sort and check for doubles
  std::sort(arguments.begin(), arguments.end());
  for (int i = 2; i < arguments.size(); ++i)
    if (arguments[i].rank == arguments[i-1].rank && arguments[i].value != arguments[i-1].value)
      Complaint(userError) << 
      "multiple values for \"" << arguments[i].id + "\" (\"" << 
      arguments[i-1].value << "\" and \"" <<
      arguments[i].value << "\").";

  // process values
  for (std::vector< Argument >::const_iterator a = arguments.begin(); a != arguments.end(); ++a)
  {
    if (a->id == "render") 
    {
      DOUBLE x, y; int g, d, r; g = 0;
      stringstream value(a->value);
      readSegmentParameters(value, x, y, g, d, r);
      if (g == -1) Complaint(userError) 
        << "generating line segment on \"render\" line must have non-zero orientation.";
      inputWord.clear(); value >> inputWord;
      if (!inputWord.empty()) Complaint(userError) 
        << "unexpected text \"" << inputWord << "\" on \"render\" line.";
      DOUBLE l = Point(x,y).length();
      if (l == 0) Complaint(userError) << "generating line segment on \"render\" line must have non-zero length.";
      x /= l; y /= l; 
      orientation = Similarity(
        Rotation(
          Point(d * x, -r * y),
          Point(d * y,  r * x)
        ),
        d == 1 ? ElevatedPoint(Point(0, 0), 0)
               : ElevatedPoint(Point(x, y), 1),       
        d
      );
      ridgeLine = 0.5 * x; // note: cannot undo ridge command, because cmnds are sorted, render first.
    }
    else if (a->id == "ridge")    ridgeLine = a->floatValue(-1,1000);
    else if (a->id == "centre")
    {
      defaultCentre = false;
      centre = ifs(a->value);
    }
    else if (a->id == "azimuth")  azimuth = a->floatValue(-360,360);
    else if (a->id == "altitude") altitude = a->floatValue(-90,90);
    else if (a->id == "distance") distance = a->floatValue(-2,2);
  }
  if (ridgeLine < 0) ridge = false;
}

void View::report(ostream& s) const 
{
  s.unsetf(std::ios_base::floatfield);
  s << "View settings:" << endl;
  orientation.print(s, "transformation       "); 
  s.unsetf(std::ios_base::floatfield);
  s << "upper ref. plane     ";
  if (ridge)         s << "starts at line x = " << ridgeLine << endl;
  else               s << "none" << endl;
  if (defaultCentre) s << "camera centre        automatic" << endl;
  else               s << "camera/zoom centre   " << centre << endl;
  s << "camera azimuth       " << azimuth << endl;
  s << "camera altitude      " << altitude << endl;
  s << "camera distance      " << distance << " (1 = where it sees all)" << endl; 
  s << endl;
}

const OptionSyntax View::syntax[9] =
{
  OptionSyntax(0, "render", true),
  OptionSyntax(1, "ridge", true),
  OptionSyntax(2, "centre", true),
  OptionSyntax(3, "altitude", true),
  OptionSyntax(4, "azimuth", true),
  OptionSyntax(5, "distance", true),
  OptionSyntax(6, "ifs", false),
  OptionSyntax(7, "view", true),
  OptionSyntax(99, "")
};

int Options::verbosity = 1;
 
Options::Options()
{
  reset();
}

void Options::reset() 
{
  view.reset();
  accuracy = 1;
  bridgeStyle = pointedBridges;
  bridgeWidth = 0.015; 
  bridgeCoreWidth = 0.010;
  clearance = 0.010;
  strength = 0.005;
  nrColourClasses = 1;
  erosionRadius = -1;
  resolution = 200;
  parapet = 0.002;
  verbosity = 1;
  supported = true;
  zoomExponent = 0; // 0 means 1, but with more efficient non-zooming algo
}

void Options::parse(const int argc, const char* argv[])
{
  // restore defaults
  reset();

  // collect command line arguments, check syntax
  std::vector< Argument > arguments;
  for (int i = 1; i < argc; ++i)
  {
    if (argv[i][0] != '-') 
    {
      arguments.push_back(Argument(0, "-i", argv[i]));
      continue;
    }
    const OptionSyntax* o = findOption(syntax, string(argv[i]));
    if (!o) Complaint(usageError) << "unrecognized option \"" << argv[i] << "\".";
    if (o->requiresArgument)
    {     
      ++i; 
      if (i == argc) Complaint(usageError) << "missing value for \"" << o->id << "\" option.";
    }
    arguments.push_back(Argument(o->rank, o->id, argv[i]));           
  }

  // sort command line arguments, check for missing or double arguments
  std::sort(arguments.begin(), arguments.end());
  if (arguments.size() < 1 || arguments[0].rank != 0) 
    Complaint(usageError) << "missing input file name.";
  if (arguments.size() > 1 && arguments[1].rank == 0)
    Complaint(usageError) << 
      "multiple inputs specified (\"" << 
      arguments[1].value << "\" and \"" <<
      arguments[2].value << "\").";
  for (int i = 2; i < arguments.size(); ++i)
    if (arguments[i].rank == arguments[i-1].rank && arguments[i].value != arguments[i-1].value)
      Complaint(usageError) << 
      "multiple values for \"" << arguments[i].id + "\" option (\"" << 
      arguments[i-1].value << "\" and \"" <<
      arguments[i].value << "\").";
 
  // process command line arguments
  for (std::vector< Argument >::const_iterator a = arguments.begin(); a != arguments.end(); ++a)
  {
    switch(a->id[1]) 
    {
      case 'i': cutId(a->value, ifsFileName, ifsName); break;
      case 's':
        if (a->value == "eroded")
        { 
          erosionRadius = 0;
          supported = true;
          parapet = 0;
        }
        else if (a->value == "floating") 
          supported = false;
        else if (a->value != "solid")
          Complaint(usageError) << "unrecognized style \"" << a->value << "\".";
        break;
      case 'r': resolution = a->floatValue(5, 10000); break;
      case 'a': accuracy = a->floatValue(0.25, 20.0); break;
      case 'b': bridgeWidth = a->floatValue(0, 0.1); break;
      case 'c': cutId(a->value, cfgFileName, cfgName); break;
      case 'f': bridgeStyle = rectangularBridges; break;
      case 'g': nrColourClasses = a->intValue(1, 1000); break;
      case 'h': clearance = a->floatValue(0, 1); break;
      case 'k': 
        if (erosionRadius < 0) Complaint(usageError) << 
          "\"-k\" option invalid for \"solid\" or \"floating\" style.";
        erosionRadius = a->intValue(0, resolution);
        break;
      case 'm': bridgeCoreWidth = a->floatValue(0, 0.1); break;
      case 'p':
        if (erosionRadius >= 0) Complaint(usageError) <<
          "\"-p\" option invalid for \"eroded\" style.";
        parapet = a->floatValue(0, 1);
        break;
      case 'q': verbosity = 0; break;
      case 'v': 
        if (verbosity != 1) Complaint(usageError) << "conflicting options \"-q\" and \"-v\"."; 
        verbosity = 2;
        break;
      case 'w': strength = a->floatValue(0, 1); break;
      case 'z': zoomExponent = a->floatValue(1, 10); break;
      default: Complaint(internalError) << "unrecognized option " << a->id << ".";
    } 
  }

  // read IFS and view configuration
  CommentedInput ifsFile(ifsFileName);
  if (ifsFile.fail()) Complaint(userError) << "cannot read file \"" << ifsFileName << "\".";
  if (!ifsName.empty() && !ifsFile.find("IFS", ifsName))
    Complaint(userError) << "IFS \"" << ifsName << "\" not found in file \"" << ifsFileName << "\".";
  if (verbosity >= 2) 
  {
    cerr << "Reading IFS";
    if (!ifsName.empty()) cerr << " \"" << ifsName << "\"";
    cerr << " from file \"" << ifsFileName << "\"..." << endl;
  } 
  ifs.parse(ifsFile);
  view.parse(ifsFile, ifs);

  if (!cfgFileName.empty() || !cfgName.empty()) 
  {
    if (cfgFileName.empty()) cfgFileName = ifsFileName;
    CommentedInput cfgFile(cfgFileName);
    if (cfgFile.fail()) Complaint(userError) << "cannot read file \"" << cfgFileName << "\".";
    if (!cfgName.empty() && !cfgFile.find("view", cfgName))
      Complaint(userError) << "View \"" << cfgName << "\" not found in file \"" << cfgFileName << "\".";
    if (verbosity >= 2) 
    {
      cerr << "Reading view";
      if (!cfgName.empty()) cerr << " \"" << cfgName << "\"";
      cerr << " from file \"" << cfgFileName << "\"..." << endl;
    } 
    view.parse(cfgFile, ifs);
  }

  if (parapet > clearance * 0.5) 
    Complaint(warning) << "Parapet height is set to more than half of the minimum vertical distance between layers (minimum tunnel height). Beware of unexpected artefacts.";
 
  if (!supported && view.ridge)
  {
    Complaint(warning) << "no upper reference plane and ridge are drawn with \"floating\" style.";
    view.ridge = false;
  }

  if (zoomExponent != 0) 
  {
    if (view.defaultCentre)
      Complaint(userError) << "rendering/view configuration does not specify centre for polynomial zoom."; 
    if (view.ridge) 
    {
      Complaint(warning) << "no upper reference plane and ridge are drawn in polynomial zoom.";
      view.ridge = false;
    }
    if (zoomExponent > 1 && !ifs.isContinuous())
      Complaint(warning) << "jumps are drawn as straight line segments despite polynomial zoom.";
  }

  ifs.applyUserSettings(*this);

  if (verbosity >= 2) 
  {
    cerr << endl;
    report(cerr);
    view.report(cerr);
  }
}

void Options::report(ostream& s) const
{
  s.unsetf(std::ios_base::floatfield);
  s << "IFS file/name        " << ifsFileName; 
  if (!ifsName.empty()) s << ':' << ifsName;
  s << endl;
  s << "View file/name       ";
  if (cfgFileName.empty()) s << "IFS default";
  else { s << cfgFileName; if (!cfgName.empty()) s << ':' << cfgName; }
  s << endl;

  s << endl;
  s << "General settings:" << endl;
  s << "zoom exponent        "; 
  if (zoomExponent == 0) s << "no zoom"; else s << zoomExponent; 
  s << endl;

  s << "style                ";
  s << (supported ? (erosionRadius < 0 ? "solid" : "eroded") : "floating");
  s << endl; 
 
  s << "resolution           " << (int) resolution      << endl;
  s << "accuracy             " << accuracy        << endl;
  s << "number of colours    " << nrColourClasses << endl;

  if (erosionRadius >= 0) 
    s << "erosion radius       " << erosionRadius   << endl;
      
  if (erosionRadius < 0) 
  {
    s << "parapet height       " << std::setprecision(4) << parapet         << endl;
    s << "min. tunnel height   " << std::setprecision(4) << clearance       << endl;
    s << "min. thickness       " << std::setprecision(4) << strength        << endl;
  }
  s << "bridge/tunnel style  ";
  if (ifs.isContinuous())
    s << "not applicable (traversal is continuous)" << endl;
  else
  {
    switch (bridgeStyle) 
    { 
      case noBridges: s << "none"; break;
      case pointedBridges: s << "pointed"; break;
      case rectangularBridges: s << "rectangular"; break;
      default: Complaint(internalError) << "bridge style undefined";
    }
    s << endl;
    s << "bridge/tunnel width  " << bridgeWidth     << endl;
    s << "brdg/tnl core width  " << bridgeCoreWidth << endl;
  }
  
  s << endl;
  s.unsetf(std::ios_base::floatfield);
} 

void Options::usage() 
{
  cerr << pftrailInfo << endl
       << "Usage: pftrail <ifsfile>[:<ifsid>] [<options>]" << endl
       << "Output: geometries and nodes libraries. Concatenate preamble, output, and postamble"
       << " to produce a COLLADA file that can be imported in Blender." << endl
       << endl
       << "<ifsfile>   name of a file with an iterated function system in extended" << endl
       << "            Ventrella format (see documentation)" << endl
       << "<ifsid>     name of the IFS, if the file contains more than one" << endl
       << "<options>   the following options are supported:" << endl
       << "   -a <number>   accuracy, use values under 1 for undersampling (default: 1)" << endl
       << "   -b <number>   maximum bridge/tunnel width (default: 0.015)" << endl
       << "   -c [<ifsfile>][:<viewid>]   name of view configuration in ifs-file" << endl
       << "                 at least the file name or the view id must be specified" << endl
       << "                 (default: the configuration that follows the curve definition)" << endl
       << "   -f            use fixed-width instead of pointed bridges" << endl
       << "   -g <number>   number of colour classes in colour gradient (default: 0)" << endl
       << "   -h <number>   minimum tunnel height (default: 0.010)" << endl
       << "   -k <number>   size of the kernel in eroded style (default: 0)" << endl
       << "   -m <number>   width of bridge/tunnel core/median strip (default: 0.010)" << endl
       << "   -r <number>   resolution: image diameter / grid cell diameter (default: 200)" << endl
       << "   -p <number>   height of parapets in solid/floating style (default: 0.002)" << endl
       << "   -q            quiet mode: less info on error output" << endl
       << "   -s <style>    \"solid\", \"floating\", or \"eroded\" (default: solid; see manual)" << endl
       << "   -v            verbose mode: more info on error output" << endl
       << "   -w <number>   minimum weight (thickness) of surfaces (default: 0.005)" << endl
       << "   -z <number>   exponent for polynomial zooming (default: 1)" << endl
       << endl;
  exit(-1);
}

const OptionSyntax Options::syntax[16] = 
{
  OptionSyntax( 1, "-s", true),
  OptionSyntax( 2, "-r", true),
  OptionSyntax( 3, "-a", true),
  OptionSyntax( 4, "-b", true),
  OptionSyntax( 5 ,"-c", true),
  OptionSyntax( 6, "-f"),
  OptionSyntax( 7, "-g", true),
  OptionSyntax( 8, "-h", true),
  OptionSyntax( 9, "-k", true),
  OptionSyntax(10, "-m", true),
  OptionSyntax(11, "-p", true),
  OptionSyntax(12, "-q"),
  OptionSyntax(13, "-v"),
  OptionSyntax(14, "-w", true),
  OptionSyntax(15, "-z", true),
  OptionSyntax(99, "") // sentinel end
}; 

void Options::cutId(const string id, string& fileName, string& name)
{
  fileName = id;
  name.erase();
  int cut = fileName.find(':');
  if (cut >= fileName.length()) return;
  name = fileName.substr(cut+1);
  fileName.erase(cut);
}

const OptionSyntax* findOption(const OptionSyntax options[], string argv)
{
  string lowerCaseArg(argv.substr(0,4)); makeLowerCase(lowerCaseArg);
  for (const OptionSyntax* o = &options[0]; ; ++o) 
  {
    if (o->id.empty()) return 0;
    if (lowerCaseArg == o->id.substr(0,4)) return o;
  }
}

Argument::Argument(int _rank, string _id, const string& _value):
  rank(_rank), id(_id), value(_value) 
{
};
 
int Argument::intValue(const int min, const int max) const
{
  int x = strtoistop(value);
  if (x < min || x > max) Complaint(usageError) 
    << "unexpected value \"" + value + "\" for \"" << id << "\" option.\n"
    << "Expected value between " << min << " and " << max << ".";
  return x;
}

DOUBLE Argument::floatValue(const DOUBLE min, const DOUBLE max) const
{
  DOUBLE x = strtofstop(value);
  if (x < min || x > max) Complaint(usageError) 
    << "unexpected value \"" + value + "\" for \"" << id << "\" option.\n"
    << "Expected value between " << min << " and " << max << ".";
  return x;
}

bool Argument::operator<(const Argument& a) const 
{ 
  return rank < a.rank; 
}


/*********************************************************
* Hex Grid (model creation)                              *
*********************************************************/

// the vertices of a hexagon of unit edge length in sample 
// or model space, vertex i is the vertex shared with 
// neighbours i-1 and i (see above for indexing of neighbours):

const Point hexagon[6] = {
  Point( 0  , (DOUBLE) -1.0 * SQRTTHIRD),
  Point( 0.5, (DOUBLE) -0.5 * SQRTTHIRD),
  Point( 0.5, (DOUBLE)  0.5 * SQRTTHIRD),
  Point( 0  , (DOUBLE)  1.0 * SQRTTHIRD),
  Point(-0.5, (DOUBLE)  0.5 * SQRTTHIRD),
  Point(-0.5, (DOUBLE) -0.5 * SQRTTHIRD)
};

Slab::Slab(): 
  bottom(0), top(0) 
{};

Slab::Slab(DOUBLE _bottom, DOUBLE _top): 
  bottom(_bottom), top(_top) 
{};

void Slab::decideWallAndParapet(
  const HexCell& neighbour, 
  std::vector< DOUBLE >& wallStartStop, DOUBLE& parapet
) const
{
  // first find highest neighbour ceiling not above current floor
  HexCell::SlabIterator s = neighbour.begin(); // int neighbourIndex = 1;
  while (s != neighbour.end() && s->bottom > top) ++s;

  DOUBLE drop = top;
  if (s != neighbour.end()) drop -= s->top;

  if (top == 1.0 || drop < noParapetDrop) parapet = 0;
  else parapet = std::min((DOUBLE) 1.0, drop / maxParapetDrop);
  
  // former check for top > 0 removed (in preparation for sound 3D models)

  if (s == neighbour.end())
  {
    wallStartStop.push_back(top);
    wallStartStop.push_back(bottom);
    return;
  }
  if (s->top < top)
  {
    wallStartStop.push_back(top);
    wallStartStop.push_back(max(s->top, bottom));
  }
  while (s->bottom > bottom)
  {
    wallStartStop.push_back(s->bottom);
    ++s;
    if (s == neighbour.end() || s->top <= bottom)
    {
      wallStartStop.push_back(bottom);
      return;
    }
    wallStartStop.push_back(s->top);
  }
}

Sample::Sample(): 
  h(0), a(sampleFill), range(-1,-1) 
{};

Sample::Sample(AnnotatedPoint p): 
  h(p.h), a(p.a), range(p.range) 
{}; 

Sample::Sample(DOUBLE _h): 
  h(_h), a(sampleFill), range(_h,_h) 
{};
  
DOUBLE Sample::weight() const
{
  return range.second - range.first;
}

DOUBLE Sample::priority() const
{
  return (a == sampleJumpCore ? weight() : 0); 
}

bool Sample::operator<(const Sample& s) const // sorts top-down
{ 
  return (h > s.h);
}

HexCell::HexCell() 
{};
  
HexCell::SlabIterator HexCell::begin() const
{
  return slabs.begin();
}

HexCell::SlabIterator HexCell::end() const
{
  return slabs.end();
}

void HexCell::calculateSlabs(
  std::vector< Sample >& samples, 
  DOUBLE strength, 
  DOUBLE clearance, 
  bool supported
) 
{
  slabs.clear();
  if (samples.empty()) return;
  bool fill = false;
    
  // maintain a stack of requested floors that have not been turned
  // into slabs yet, in order of decreasing elevation and priority
  std::stack< Sample > pending;
  for (std::vector< Sample >::const_iterator s = samples.begin(); ; ++s)
  {
    // remove any lower-priority conflicting floors from the stack
    if (s != samples.end())
      while (
        !pending.empty() && 
        pending.top().h < s->h + clearance + strength &&
        pending.top().priority() < s->priority()
      ) 
        pending.pop();
 
    // if the lowest floor on the stack (top of the stack) does not
    // conflict with the new floor, turn the stack into a slab and
    // empty the stack    
    if (!pending.empty())
      if (s == samples.end() || pending.top().h >= s->h + clearance + strength)
      {
        DOUBLE bottom = pending.top().h - strength;
        DOUBLE top;
        while (!pending.empty()) { top = pending.top().h; pending.pop(); }
        slabs.push_back(Slab(bottom, top));
      }

    if (s == samples.end()) break;

    pending.push(*s);
    fill = (fill || s->a == sampleFill);
  }
  if (supported && fill) slabs.back().bottom = foundation;
}

void HexCell::calculateMeanFloor(const std::vector< Sample >& samples)
{
  static const DOUBLE practicallyZero = 0.0000001;
  slabs.clear();
  DOUBLE sum = 0; int count = 0; DOUBLE weight = 0;
  for (std::vector< Sample >::const_iterator s = samples.begin(); s != samples.end(); ++s)
  {
    if (s->weight() < weight - practicallyZero) continue;
    if (s->weight() > weight + practicallyZero) 
    { 
      sum = 0; count = 0; weight = s->weight(); 
    }
    sum += s->h; ++count;      
  }
  if (count == 0) slabs.push_back(Slab(foundation, -0.001)); // just below the ground plane
  else slabs.push_back(Slab(foundation, sum/count));
}

void HexCell::copySlab(const HexCell& h)
{  
  slabs.push_back(h.slabs[0]);
}
 
void HexCell::calculateMeanFloor()
{
  if (slabs.empty()) return;
  DOUBLE sum = 0;
  for (SlabIterator s = begin(); s != end(); ++s) sum += s->top;
  sum /= slabs.size();
  slabs.clear();
  slabs.push_back(Slab(foundation, sum));
}

Mesh::Mesh(const string& _material): 
  material(_material), nrPending(0) 
{};
  
void Mesh::setMaterial(const string& _material)
{
  material = _material;
}

const std::vector< string >& Mesh::materials()
{
  return theMaterials;
}

void Mesh::addVertex(MeshVertex& v, int normalIndex)
{
  if (v.index < 0) 
  {
    v.index = coordinates.size();
    coordinates.push_back(v.point);
  }
  polygons.push_back(v.index);
  polygons.push_back(normalIndex);
  ++nrPending;
}

void Mesh::addVertex(const Point& p, DOUBLE h, int normalIndex)
{
  polygons.push_back(coordinates.size());
  polygons.push_back(normalIndex);
  coordinates.push_back(ElevatedPoint(p, h));
  ++nrPending;
}

void Mesh::addVertex(DOUBLE x, DOUBLE y, DOUBLE h, int normalIndex)
{
  addVertex(Point(x, y), h, normalIndex);
}

void Mesh::endPolygon()
{
  polygonSizes.push_back(nrPending);
  nrPending = 0;
}

bool Mesh::full() const 
{
  return (coordinates.size() > maxMeshFill);
}

void Mesh::clear() 
{ 
  if (nrPending > 0) 
    Complaint(internalError) << "clear() called on a mesh with an unfinished polygon.";
  cout << "  <geometry id=\"landscape" << theMaterials.size() << "\">" << endl
       << "    <mesh>" << endl
       << "      <source id=\"Coordinates" << theMaterials.size() << "\">" << endl
       << "        <float_array id=\"CoordinatesArray" << theMaterials.size() 
       << "\" count=\"" << coordinates.size() * 3 << "\">" << endl;
  for (int i = 0; i < coordinates.size(); ++i)
  {
    cout << "          " << fixed << setprecision(outputPrecision) << -coordinates[i].y
         << ' '          << fixed << setprecision(outputPrecision) <<  coordinates[i].h
         << ' '          << fixed << setprecision(outputPrecision) << -coordinates[i].x
         << endl;
  }
  cout << "        </float_array>" << endl
       << "        <technique_common>" << endl
       << "          <accessor source=\"#CoordinatesArray" << theMaterials.size() 
       << "\" count=\"" << coordinates.size() << "\" stride=\"3\">" << endl
       << "            <param name=\"X\" type=\"float\" />" << endl
       << "            <param name=\"Y\" type=\"float\" />" << endl
       << "            <param name=\"Z\" type=\"float\" />" << endl
       << "          </accessor>" << endl
       << "        </technique_common>" << endl
       << "      </source>" << endl;
  cout << "      <source id=\"Normals" << theMaterials.size() << "\">" << endl
       << "        <float_array id=\"NormalsArray" << theMaterials.size() 
       << "\" count=\"24\">" << endl
       << "           0.86602540  0 -0.5" << endl
       << "           0           0 -1  " << endl
       << "          -0.86602540  0 -0.5" << endl
       << "          -0.86602540  0  0.5" << endl
       << "           0           0  1  " << endl
       << "           0.86602540  0  0.5" << endl
       << "           0           1  0  " << endl
       << "           0          -1  0  " << endl 
       << "        </float_array>" << endl
       << "        <technique_common>" << endl
       << "          <accessor source=\"#NormalsArray" << theMaterials.size() 
       << "\" count=\"8\" stride=\"3\">" << endl
       << "            <param name=\"X\" type=\"float\" />" << endl
       << "            <param name=\"Y\" type=\"float\" />" << endl
       << "            <param name=\"Z\" type=\"float\" />" << endl
       << "          </accessor>" << endl
       << "        </technique_common>" << endl
       << "      </source>" << endl;
  cout << "      <vertices id=\"Vertices" << theMaterials.size() << "\">" << endl
       << "        <input semantic=\"POSITION\" source=\"#Coordinates" << theMaterials.size() << "\"/>" << endl
       << "      </vertices>" << endl
       << "      <polylist count=\"" << polygonSizes.size() << "\" material=\"SOLID" << theMaterials.size() << "\">" << endl
       << "        <input semantic=\"VERTEX\" source=\"#Vertices" << theMaterials.size() << "\" offset=\"0\"/>" << endl
       << "        <input semantic=\"NORMAL\" source=\"#Normals" << theMaterials.size() << "\" offset=\"1\"/>" << endl
       << "        <vcount>" << endl;
  for (int i = 0; i < polygonSizes.size(); ++i) cout << ' ' << polygonSizes[i];
  cout << endl 
       << "        </vcount>" << endl;
  int offset = 0;
  for (int i = 0; i < polygonSizes.size(); ++i)
  {
    cout << "        <p>";
    for (int j = 0; j < 2*polygonSizes[i]; ++j) cout << ' ' << polygons[offset+j];
    cout << " </p>" << endl;
    offset += 2*polygonSizes[i];
  }
  cout << "      </polylist>" << endl
       << "    </mesh>" << endl
       << "  </geometry>" << endl;
  theMaterials.push_back(material);
  coordinates.clear();
  polygonSizes.clear();
  polygons.clear();
}

Mesh::~Mesh() 
{
  if (!coordinates.empty()) clear();
}

std::vector< string > Mesh::theMaterials;

MeshVertex::MeshVertex(Point p, DOUBLE h):
  point(p, h), index(-1) 
{};

int HexGrid::gridIndex(const int r, const int c) const {
  return (r + gridIndexRadius) * gridWidth + (c + gridIndexRadius);
}

int HexGrid::gridIndex(const Point& p) const {
  TrianglePoint< int > t = 
    round(TrianglePoint< DOUBLE >(p * cellsPerUnit));
  for (int j = 0; j < 2; ++j)
    if (t[j] > gridIndexRadius || t[j] < -gridIndexRadius) 
      Complaint(internalError) << 
        "grid index (" << t[0] << ',' << t[1] << ") for point " << p << " out of bounds.";
  return gridIndex(t[0], t[1]);
}

Point HexGrid::location(int gridIndex) const
{
  TrianglePoint< int > t(0,
    gridIndex / gridWidth - gridIndexRadius,
    gridIndex % gridWidth - gridIndexRadius
  );
  return ((Point) t) * cellSpacing;
}

HexGrid::HexGrid(
  const Options& options,
  const DOUBLE edgeLength,
  const DOUBLE ridgeLine, 
  const std::vector< AnnotatedPoint >& pp
)
{
  // ifs.generatePoints guarantees that sample density suffices for grid
  // cells with diameter options.ifs.diameterLowerBound() / resolution, 
  // which means 
  // 2.0 * SQRTTHIRD * resolution / options.ifs.diameterLowerBound()
  // cells per unit along an edge/row/column of the grid
  cellsPerUnit = SQRTTHIRD * 2.0 * options.resolution / options.ifs.diameterLowerBound();

  gridIndexRadius = round(edgeLength * cellsPerUnit * 0.5) 
    + max(options.erosionRadius, 0) 
    + 2  // for rounding errors
    + 1; // for sentinel boundary (not rendered)
  gridWidth = 2 * gridIndexRadius + 1;
  gridSize = gridWidth * gridWidth;
  if (Options::verbosity >= 2) 
    cerr << "Rendering grid size: " << gridWidth << " x " << gridWidth << endl;

  cellSpacing = cornerDepth / (gridIndexRadius - 1); 
  for (int i = 0; i < 6; ++i) cellShape[i] = hexagon[i] * cellSpacing;

  indexDistanceToNeighbour[0] = gridWidth;
  indexDistanceToNeighbour[1] = gridWidth-1;
  indexDistanceToNeighbour[2] = -1;
  indexDistanceToNeighbour[3] = -gridWidth;
  indexDistanceToNeighbour[4] = 1-gridWidth;
  indexDistanceToNeighbour[5] = 1;

  if (Options::verbosity >= 1) cerr << "Sorting sample points into grid cells..." << endl;
  std::vector< Sample >* samples = new std::vector< Sample >[gridSize];
  for (int i = 0; i < pp.size(); ++i)
  {
    AnnotatedPoint p = pp[i];
    // limit elevation to very slightly below 1, to be able to tell
    // the difference with the upper reference plane:
    p.h = std::min((DOUBLE) 0.9999, p.h);
    samples[gridIndex(p)].push_back(Sample(p));
  }

  if (Options::verbosity >= 1) cerr << "Calculating camera position..." << endl;
  // first figure out where the camera needs to be to see everything:
  const DOUBLE t = 0.33; // tan of (field of view / 2)
  const DOUBLE degtorad = 3.14159265358979323846 / 180.0;
  const DOUBLE sinazm = sin(options.view.azimuth * degtorad);
  const DOUBLE cosazm = cos(options.view.azimuth * degtorad);
  const DOUBLE sinalt = sin(options.view.altitude * degtorad);
  const DOUBLE cosalt = cos(options.view.altitude * degtorad);
  DOUBLE topEdgeConstraint = 100;    // t*camera.x - camera.h < topEdgeConstraint
  DOUBLE bottomEdgeConstraint = 100; // t*camera.x + camera.h < bottomEdgeConstraint
  DOUBLE leftEdgeConstraint = 100;   // t*camera.x - camera.y < leftEdgeConstraint
  DOUBLE rightEdgeConstraint = 100;  // t*camera.x + camera.y < rightEdgeConstraint
  for (int r = -(gridIndexRadius-1); r <= (gridIndexRadius-1); ++r)
  {
    for (int c = -(gridIndexRadius-1); c <= (gridIndexRadius-1); ++c)
    {
      int i = gridIndex(r, c);
      if (samples[i].empty()) continue;
      std::sort(samples[i].begin(), samples[i].end());
      Point p = location(i);
      ElevatedPoint e[2];
      e[0] = ElevatedPoint(p, samples[i].front().h);
      e[1] = ElevatedPoint(p, samples[i].back().h);
      for (int j = 0; j < 2; ++j)
      {
        // rotate scene so that camera looks into direction of positive x-axis
        ElevatedPoint tmp(
          e[j].x*cosazm + e[j].y*sinazm, -e[j].x*sinazm + e[j].y*cosazm, e[j].h);
        e[j] = ElevatedPoint(
          tmp.x*cosalt - tmp.h*sinalt, tmp.y, tmp.x*sinalt + tmp.h*cosalt);
        DOUBLE tqx = t * e[j].x;
        topEdgeConstraint = std::min(topEdgeConstraint, tqx - e[j].h);
        bottomEdgeConstraint = std::min(bottomEdgeConstraint, tqx + e[j].h);
        leftEdgeConstraint = std::min(leftEdgeConstraint, tqx - e[j].y);
        rightEdgeConstraint = std::min(rightEdgeConstraint, tqx + e[j].y);
      }
    }
  }
  camera.y = (options.view.defaultCentre ? 
      (rightEdgeConstraint - leftEdgeConstraint) / 2 : 0);
  camera.h = (options.view.defaultCentre ? 
      (bottomEdgeConstraint - topEdgeConstraint) / 2 : options.view.centre.h);
  camera.x = options.view.distance * std::min(
    std::min(topEdgeConstraint + camera.h,
    bottomEdgeConstraint - camera.h),
    std::min(leftEdgeConstraint + camera.y,
    rightEdgeConstraint - camera.y)
  ) / t;

  discreteRidge = -3 * gridIndexRadius;
  if (options.view.ridge) {
    if (Options::verbosity >= 1) cerr << "Constructing upper reference plane..." << endl;
    discreteRidge = 2 * ridgeLine * cellsPerUnit + 1;
    for (int r = -gridIndexRadius+1; r <= gridIndexRadius; ++r)
      for (int c = min(gridIndexRadius-1, r-discreteRidge); 
        c >= -gridIndexRadius; --c)
      {
        std::vector< Sample >& cell = samples[gridIndex(r, c)];
        if (cell.empty()) 
          cell.push_back(Sample(1.0));
        else // make sure the space under the lowest sample is filled:
          if (cell.back().a != sampleFill) 
            cell.push_back(Sample(cell.back().h));
      }
  }

  if (Options::verbosity >= 1) cerr << "Calculating floor levels..." << endl;
  grid = new HexCell[gridSize];
  if (options.erosionRadius < 0) {
    for (int i = 0; i < gridSize; ++i)
      // the next step requires the samples to be sorted;
      // this has already been done while computing the camera position
      grid[i].calculateSlabs(
        samples[i], 
        options.strength, 
        (options.supported ? options.clearance : 0), 
        options.supported
      );
  }
  else {
    for (int i = 0; i < gridSize; ++i) grid[i].calculateMeanFloor(samples[i]);
    if (options.erosionRadius > 0) 
      if (Options::verbosity >= 2) 
        cerr << "Smoothing " << options.erosionRadius << " times..." << endl;
    for (int k = 0; k < options.erosionRadius; ++k) {
      for (int i = 0; i < gridSize; ++i) {
        for (int j = 0; j < 6; ++j) {
          int nbr = i + indexDistanceToNeighbour[j];
          int r = nbr / gridWidth; if (r <= 1 || r >= gridWidth - 2) continue;
          int c = nbr % gridWidth; if (c <= 1 || c >= gridWidth - 2) continue;
          grid[nbr].copySlab(grid[i]);
        }
      }
      for (int i = 0; i < gridSize; ++i) grid[i].calculateMeanFloor();
    }
  }
  delete [] samples;
}

// array indexed modulo n, used by HexGrid::draw:
template < class T, int n > class WrapArray
{
    private:
  T a[n];
    public:
  WrapArray() {};
  T& operator[](int i)
  {
    return a[(i+n)%n];
  } 
};
 
void HexGrid::draw(int nrColours, DOUBLE parapetHeight) 
{
  Mesh terrain[nrColours];
  string colourNames[nrColours+1];
  for (int i = 0; i < nrColours; ++i)
  {
    stringstream name;
    name << "PathMaterial" << i;
    terrain[i].setMaterial(name.str());
  }
  Mesh parapets("ParapetMaterial");
  Mesh walls("WallMaterial");
  Mesh upperPlane("UpperPlaneMaterial");
  for (int r = -(gridIndexRadius-1); r <= (gridIndexRadius-1); ++r)
  {
    for (int c = -(gridIndexRadius-1); c <= (gridIndexRadius-1); ++c)
    {
      int h = gridIndex(r, c);
      Point p = location(h);
      for (HexCell::SlabIterator s = grid[h].begin(); s != grid[h].end(); ++s)
      { 
        int colour = (int) (s->top * nrColours);
        if (colour < 0) colour = 0; // just to catch rounding errors 
        Mesh& surface = (colour < nrColours ? terrain[colour] : upperPlane);
        if (surface.full()) surface.clear();          
        if (parapets.full()) parapets.clear();
        if (walls.full()) walls.clear();

        // figure out what walls and parapets are needed:
        WrapArray< DOUBLE, 6 > parapetFraction;
        WrapArray< std::vector< DOUBLE >, 6 > wallStartStop;
        for (int i = 0; i < 6; ++i)
        {
          // parapet is constructed where there is a drop > 0.001; 
          // ideal height between 0 and parapetHeight (for drop >= 0.02),
          // fraction of parapetHeight stored in parapetFraction
          s->decideWallAndParapet(
            grid[h + indexDistanceToNeighbour[i]], 
            wallStartStop[i], parapetFraction[i]);
          if (parapetHeight <= 0) parapetFraction[i] = 0;
        }

        // define vertices
        WrapArray< MeshVertex, 6 > surfaceTopOutside;
        WrapArray< MeshVertex, 6 > surfaceTopInside;
        WrapArray< MeshVertex, 6 > parapetTopOutside;
        WrapArray< MeshVertex, 6 > parapetTopInside;
        WrapArray< MeshVertex, 6 > parapetParapet;
        WrapArray< MeshVertex, 6 > wallBottom;
        WrapArray< MeshVertex, 6 > wallTopOutside;
        WrapArray< MeshVertex, 6 > wallParapet;
        for (int i = 0; i < 6; ++i)
        {
          // if parapet top is between two parapets, take max of their ideal heights
          DOUBLE postFraction = max(parapetFraction[i-1], parapetFraction[i]);
          DOUBLE postHeight = postFraction * parapetHeight;
          DOUBLE postDepth = 
            minSurfaceToParapetRatio / (postFraction + minSurfaceToParapetRatio); 
          surfaceTopOutside[i] = MeshVertex(p + cellShape[i], s->top);
          surfaceTopInside[i] = MeshVertex(p + cellShape[i] * postDepth, s->top);
          parapetTopOutside[i] = MeshVertex(p + cellShape[i], s->top);
          parapetTopInside[i] = MeshVertex(p + cellShape[i] * postDepth, s->top);
          parapetParapet[i] = MeshVertex(p + cellShape[i], s->top + postHeight);
          wallBottom[i] = MeshVertex(p + cellShape[i], s->bottom);
          wallTopOutside[i] = MeshVertex(p + cellShape[i], s->top);
          wallParapet[i] = MeshVertex(p + cellShape[i], s->top + postHeight);
        }

        // construct the top
        for (int i = 0; i < 6; ++i)
        {
          if (parapetFraction[i-1] == 0) {
            surface.addVertex(surfaceTopOutside[i], 6);
            if (parapetFraction[i] != 0) surface.addVertex(surfaceTopInside[i], 6);
          } else {
            surface.addVertex(surfaceTopInside[i], 6);
            if (parapetFraction[i] == 0) surface.addVertex(surfaceTopOutside[i], 6);
          }
        }
        surface.endPolygon();

        // construct parapets
        for (int i = 0; i < 6; ++i) {
          if (parapetFraction[i] == 0) continue;            
          if (parapetFraction[i-1] == 0) {
            parapets.addVertex(parapetTopInside[i],  (i+4)%6);
            parapets.addVertex(parapetParapet[i],    (i+4)%6);
            parapets.addVertex(parapetTopOutside[i], (i+4)%6);
            parapets.endPolygon();
          }
          parapets.addVertex(parapetTopInside[i],   (i+3)%6);
          parapets.addVertex(parapetTopInside[i+1], (i+3)%6);
          parapets.addVertex(parapetParapet[i+1],   (i+3)%6);
          parapets.addVertex(parapetParapet[i],     (i+3)%6);
          parapets.endPolygon();
          if (parapetFraction[i+1] == 0) {
            parapets.addVertex(parapetTopInside[i+1],  (i+2)%6);
            parapets.addVertex(parapetParapet[i+1],    (i+2)%6);
            parapets.addVertex(parapetTopOutside[i+1], (i+2)%6);
            parapets.endPolygon();
          }
        }

        // construct the vertical walls
        for (int i = 0; i < 6; ++i)
          for (std::vector< DOUBLE >::const_iterator w = wallStartStop[i].begin();
            w != wallStartStop[i].end(); ++w)
          {
            // top vertices
            if (*w == s->top)
            {
              if (parapetFraction[i] == 0)
              {
                walls.addVertex(wallTopOutside[i+1], i);
                walls.addVertex(wallTopOutside[i],   i);
              } 
              else
              {
                walls.addVertex(wallParapet[i+1], i);
                walls.addVertex(wallParapet[i],   i);
              }
            }
            else
            {
              walls.addVertex(p + cellShape[(i+1)%6], *w, i);
              walls.addVertex(p + cellShape[i],       *w, i);
            }

            // bottom vertices
            ++w;
            if (*w == s->bottom)
            {
              walls.addVertex(wallBottom[i],   i);
              walls.addVertex(wallBottom[i+1], i);
            }
            else
            {
              walls.addVertex(p + cellShape[i],       *w, i);
              walls.addVertex(p + cellShape[(i+1)%6], *w, i);
            }
            walls.endPolygon();
          }

        // construct the bottom (only if above the reference plane)
        if (s->bottom > 0)
        {
          for (int i = 0; i < 6; ++i) walls.addVertex(wallBottom[i], 7);
          walls.endPolygon();
        }
      }   
    }
  }
  Mesh lowerPlane("LowerPlaneMaterial");
  lowerPlane.addVertex(-100,  100, 0, 6);
  lowerPlane.addVertex( 100,  100, 0, 6);
  lowerPlane.addVertex( 100, -100, 0, 6);
  lowerPlane.addVertex(-100, -100, 0, 6);
  lowerPlane.endPolygon();
  if (discreteRidge != -3 * gridIndexRadius)
  {
    DOUBLE ridgeLine = discreteRidge * 0.5 * cellSpacing;
    DOUBLE halfGapWidth = SQRTTHREE * (cornerDepth - abs(ridgeLine));
    walls.addVertex(ridgeLine, 100, 0, 4); 
    walls.addVertex(ridgeLine, 100, 1, 4);
    if (halfGapWidth > 0)
    { 
      walls.addVertex(ridgeLine, halfGapWidth, 1, 4); 
      walls.addVertex(ridgeLine, halfGapWidth, 0, 4); 
      walls.endPolygon();
      walls.addVertex(ridgeLine, -halfGapWidth, 0, 4); 
      walls.addVertex(ridgeLine, -halfGapWidth, 1, 4); 
    }
    walls.addVertex(ridgeLine, -100, 1, 4); 
    walls.addVertex(ridgeLine, -100, 0, 4); 
    walls.endPolygon();
    upperPlane.addVertex(ridgeLine, 100, 1, 6);
    upperPlane.addVertex(ridgeLine, halfGapWidth, 1, 6);
    if (discreteRidge < 0)
      upperPlane.addVertex(0, SQRTTHREE * cornerDepth, 1, 6);
    if (ridgeLine < cornerDepth)
      upperPlane.addVertex(cornerDepth, 0, 1, 6);
    if (discreteRidge < 0)
      upperPlane.addVertex(0, -SQRTTHREE * cornerDepth, 1, 6);
    upperPlane.addVertex(ridgeLine, -halfGapWidth, 1, 6); 
    upperPlane.addVertex(ridgeLine, -100, 1, 6);
    upperPlane.addVertex(100, -100, 1, 6);
    upperPlane.addVertex(100, 100, 1, 6);
    upperPlane.endPolygon();
  }
}

ElevatedPoint HexGrid::cameraPosition() const
{
  return camera;
}

HexGrid::~HexGrid() 
{ 
  delete [] grid; 
}


/*********************************************************
* main                                                   *
*********************************************************/

int main(const int argc, const char* argv[]) {
  if (argc < 2) Options::usage();
  Options options;
  options.parse(argc, argv);

  if (Options::verbosity >= 1) cerr << "Sampling the traversal..." << endl;
  std::vector< AnnotatedPoint > pp; 
  options.ifs.generatePoints(pp, options.view.orientation);

  ElevatedPoint centre = options.view.orientation(options.view.centre);

  if (Options::verbosity >= 1) cerr << "Centering traversal..." << endl;
  Point gridCentre; DOUBLE edgeLength = 0;
  if (options.view.defaultCentre) 
    findMinimumBoundingDiamond(pp, gridCentre, edgeLength);
  else {
    gridCentre = centre;
    for (int i = 0; i < pp.size(); ++i)
    {
      TrianglePoint< DOUBLE > p(pp[i] - gridCentre);
      for (int j = 0; j < 2; ++j)
      edgeLength = std::max(edgeLength, abs(p[j])*2);
    }
  }
  for (int i = 0; i < pp.size(); ++i) pp[i] -= gridCentre;
  DOUBLE ridgeLine = options.view.ridgeLine - gridCentre.x;

  HexGrid h(options, edgeLength, ridgeLine, pp); 

  cout << "<!--" << endl;
  cout << "Path-filling trail model generated by " << pftrailInfo << endl;
  cout << endl;
  options.report(cout);
  options.ifs.report(cout);
  options.view.report(cout);
  cout << "-->" << endl;

  if (Options::verbosity >= 1) cerr << "Computing facets of 3D model..." << endl;
  cout << "<library_geometries>" << endl;
  h.draw(options.nrColourClasses, options.parapet);
  cout << "</library_geometries>" << endl << endl;

  if (Options::verbosity >= 2) cerr << "Writing COLLADA nodes library..." << endl;
  cout << "<library_nodes>" << endl;
  cout << "  <node id=\"Landscape\">" << endl;
  for (int i = 0; i < Mesh::materials().size(); ++i)
    cout << "    <node>" << endl
         << "      <instance_geometry url=\"#landscape" << i << "\">" << endl
         << "        <bind_material>" << endl
         << "                <technique_common>" << endl
         << "                  <instance_material symbol=\"SOLID" << i << "\" target=\"#"
         << Mesh::materials()[i] << "\"/>" << endl
         << "                </technique_common>" << endl
         << "        </bind_material>" << endl
         << "      </instance_geometry>" << endl
         << "    </node>" << endl;
  cout << "  </node>" << endl;
  cout << "  <node id=\"HighCamera\">" << endl
       << "    <rotate> 0 1 0 " << options.view.azimuth << " </rotate>" << endl
       << "    <rotate> 1 0 0 " << -options.view.altitude << " </rotate>" << endl
       << "    <translate> "
       << fixed << setprecision(outputPrecision) << -h.cameraPosition().y << ' ' 
       << fixed << setprecision(outputPrecision) <<  h.cameraPosition().h << ' ' 
       << fixed << setprecision(outputPrecision) << -h.cameraPosition().x
       << " </translate>" << endl
       << "    <instance_camera url=\"#TheCamera\" />" << endl
       << "  </node>" << endl
       << "</library_nodes>" << endl;
  if (Options::verbosity >= 1) cerr << "Done." << endl;
  cout << "<!-- End of pftrail output. Append postamble here to complete the file. -->" << endl;
}


/*********************************************************
**********************************************************
**                                                      **
** HISTORY                                              **
**                                                      **
**********************************************************
**********************************************************

2020-03-31 fixed strtofstop, strtoistop, and IFS::operator() to ignore trailing whitespace

*/

