snake.cpp

#include <algorithm>
#include <iostream>

#include "snake.h"

#include <mapbox/geometry.hpp>
#include <mapbox/polylabel.hpp>

#include <boost/geometry.hpp>
#include <boost/geometry/geometries/adapted/boost_tuple.hpp>
#include <boost/geometry/geometries/box.hpp>
#include <boost/geometry/geometries/polygon.hpp>

#include "clipper/clipper.hpp"
#define CLIPPER_SCALE 1000000

#include "ortools/constraint_solver/routing.h"
#include "ortools/constraint_solver/routing_enums.pb.h"
#include "ortools/constraint_solver/routing_index_manager.h"
#include "ortools/constraint_solver/routing_parameters.h"

using namespace operations_research;

#ifndef NDEBUG
//#define SNAKE_SHOW_TIME
#endif

namespace bg = boost::geometry;
namespace trans = bg::strategy::transform;

BOOST_GEOMETRY_REGISTER_BOOST_TUPLE_CS(bg::cs::cartesian)

namespace snake {
static const IntType stdScale = 1000000;
//=========================================================================
// Geometry stuff.
//=========================================================================

void polygonCenter(const FPolygon &polygon, FPoint &center) {
  using namespace mapbox;
  if (polygon.outer().empty())
    return;
  geometry::polygon<double> p;
  geometry::linear_ring<double> lr1;
  for (size_t i = 0; i < polygon.outer().size(); ++i) {
    geometry::point<double> vertex(polygon.outer()[i].get<0>(),
                                   polygon.outer()[i].get<1>());
    lr1.push_back(vertex);
  }
  p.push_back(lr1);
  geometry::point<double> c = polylabel(p);

  center.set<0>(c.x);
  center.set<1>(c.y);
}

bool minimalBoundingBox(const FPolygon &polygon, BoundingBox &minBBox) {
  /*
  Find the minimum-area bounding box of a set of 2D points

  The input is a 2D convex hull, in an Nx2 numpy array of x-y co-ordinates.
  The first and last points points must be the same, making a closed polygon.
  This program finds the rotation angles of each edge of the convex polygon,
  then tests the area of a bounding box aligned with the unique angles in
  90 degrees of the 1st Quadrant.
  Returns the

  Tested with Python 2.6.5 on Ubuntu 10.04.4 (original version)
  Results verified using Matlab

  Copyright (c) 2013, David Butterworth, University of Queensland
  All rights reserved.

  Redistribution and use in source and binary forms, with or without
  modification, are permitted provided that the following conditions are met:

      * Redistributions of source code must retain the above copyright
        notice, this list of conditions and the following disclaimer.
      * Redistributions in binary form must reproduce the above copyright
        notice, this list of conditions and the following disclaimer in the
        documentation and/or other materials provided with the distribution.
      * Neither the name of the Willow Garage, Inc. nor the names of its
        contributors may be used to endorse or promote products derived from
        this software without specific prior written permission.

  THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
  AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
  IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
  ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
  LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
  CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
  SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
  INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
  CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
  ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
  POSSIBILITY OF SUCH DAMAGE.
  */

  if (polygon.outer().empty() || polygon.outer().size() < 3)
    return false;
  FPolygon convex_hull;
  bg::convex_hull(polygon, convex_hull);

  // cout << "Convex hull: " << bg::wkt<BoostPolygon2D>(convex_hull) << endl;

  //# Compute edges (x2-x1,y2-y1)
  std::vector<FPoint> edges;
  const auto &convex_hull_outer = convex_hull.outer();
  for (long i = 0; i < long(convex_hull_outer.size()) - 1; ++i) {
    FPoint p1 = convex_hull_outer.at(i);
    FPoint p2 = convex_hull_outer.at(i + 1);
    double edge_x = p2.get<0>() - p1.get<0>();
    double edge_y = p2.get<1>() - p1.get<1>();
    edges.push_back(FPoint{edge_x, edge_y});
  }

  //    cout << "Edges: ";
  //    for (auto e : edges)
  //        cout << e.get<0>() << " " << e.get<1>() << ",";
  //    cout << endl;

  // Calculate unique edge angles  atan2(y/x)
  double angle_scale = 1e3;
  std::set<long> angles_long;
  for (auto vertex : edges) {
    double angle = std::fmod(atan2(vertex.get<1>(), vertex.get<0>()), M_PI / 2);
    angle =
        angle < 0 ? angle + M_PI / 2 : angle; // want strictly positive answers
    angles_long.insert(long(round(angle * angle_scale)));
  }
  std::vector<double> edge_angles;
  for (auto a : angles_long)
    edge_angles.push_back(double(a) / angle_scale);

  //    cout << "Unique angles: ";
  //    for (auto e : edge_angles)
  //        cout << e*180/M_PI << ",";
  //    cout << endl;

  double min_area = std::numeric_limits<double>::infinity();
  // Test each angle to find bounding box with smallest area
  // print "Testing", len(edge_angles), "possible rotations for bounding box...
  // \n"
  for (double angle : edge_angles) {

    trans::rotate_transformer<bg::degree, double, 2, 2> rotate(angle * 180 /
                                                               M_PI);
    FPolygon hull_rotated;
    bg::transform(convex_hull, hull_rotated, rotate);
    // cout << "Convex hull rotated: " << bg::wkt<BoostPolygon2D>(hull_rotated)
    // << endl;

    bg::model::box<FPoint> box;
    bg::envelope(hull_rotated, box);
    //        cout << "Bounding box: " <<
    //        bg::wkt<bg::model::box<BoostPoint2D>>(box) << endl;

    //# print "Rotated hull points are \n", rot_points
    FPoint min_corner = box.min_corner();
    FPoint max_corner = box.max_corner();
    double min_x = min_corner.get<0>();
    double max_x = max_corner.get<0>();
    double min_y = min_corner.get<1>();
    double max_y = max_corner.get<1>();
    //        cout << "min_x: " << min_x << endl;
    //        cout << "max_x: " << max_x << endl;
    //        cout << "min_y: " << min_y << endl;
    //        cout << "max_y: " << max_y << endl;

    // Calculate height/width/area of this bounding rectangle
    double width = max_x - min_x;
    double height = max_y - min_y;
    double area = width * height;
    //        cout << "Width: " << width << endl;
    //        cout << "Height: " << height << endl;
    //        cout << "area: " << area << endl;
    //        cout << "angle: " << angle*180/M_PI << endl;

    // Store the smallest rect found first (a simple convex hull might have 2
    // answers with same area)
    if (area < min_area) {
      min_area = area;
      minBBox.angle = angle;
      minBBox.width = width;
      minBBox.height = height;

      minBBox.corners.clear();
      minBBox.corners.outer().push_back(FPoint{min_x, min_y});
      minBBox.corners.outer().push_back(FPoint{min_x, max_y});
      minBBox.corners.outer().push_back(FPoint{max_x, max_y});
      minBBox.corners.outer().push_back(FPoint{max_x, min_y});
      minBBox.corners.outer().push_back(FPoint{min_x, min_y});
    }
    // cout << endl << endl;
  }

  // Transform corners of minimal bounding box.
  trans::rotate_transformer<bg::degree, double, 2, 2> rotate(-minBBox.angle *
                                                             180 / M_PI);
  FPolygon rotated_polygon;
  bg::transform(minBBox.corners, rotated_polygon, rotate);
  minBBox.corners = rotated_polygon;

  return true;
}

void offsetPolygon(const FPolygon &polygon, FPolygon &polygonOffset,
                   double offset) {
  bg::strategy::buffer::distance_symmetric<double> distance_strategy(offset);
  bg::strategy::buffer::join_miter join_strategy(3);
  bg::strategy::buffer::end_flat end_strategy;
  bg::strategy::buffer::point_square point_strategy;
  bg::strategy::buffer::side_straight side_strategy;

  bg::model::multi_polygon<FPolygon> result;

  bg::buffer(polygon, result, distance_strategy, side_strategy, join_strategy,
             end_strategy, point_strategy);

  if (result.size() > 0)
    polygonOffset = result[0];
}

void graphFromPolygon(const FPolygon &polygon, const FLineString &vertices,
                      Matrix<double> &graph) {
  size_t n = graph.n();

  for (size_t i = 0; i < n; ++i) {
    FPoint v1 = vertices[i];
    for (size_t j = i + 1; j < n; ++j) {
      FPoint v2 = vertices[j];
      FLineString path{v1, v2};

      double distance = 0;
      if (!bg::within(path, polygon))
        distance = std::numeric_limits<double>::infinity();
      else
        distance = bg::length(path);
      graph(i, j) = distance;
      graph(j, i) = distance;
    }
  }
}

bool toDistanceMatrix(Matrix<double> &graph) {
  size_t n = graph.n();

  auto distance = [&graph](size_t i, size_t j) -> double {
    return graph(i, j);
  };

  for (size_t i = 0; i < n; ++i) {
    for (size_t j = i + 1; j < n; ++j) {
      double d = graph(i, j);
      if (!std::isinf(d))
        continue;
      std::vector<size_t> path;
      if (!dijkstraAlgorithm(n, i, j, path, d, distance)) {
        return false;
      }
      //            cout << "(" << i << "," << j << ") d: " << d << endl;
      //            cout << "Path size: " << path.size() << endl;
      //            for (auto idx : path)
      //                cout << idx << " ";
      //            cout << endl;
      graph(i, j) = d;
      graph(j, i) = d;
    }
  }
  return true;
}

bool tiles(const FPolygon &area, Length tileHeight, Length tileWidth,
           Area minTileArea, std::vector<FPolygon> &tiles, BoundingBox &bbox,
           string &errorString) {
  if (area.outer().empty() || area.outer().size() < 4) {
    errorString = "Area has to few vertices.";
    return false;
  }

  if (tileWidth <= 0 * bu::si::meter || tileHeight <= 0 * bu::si::meter ||
      minTileArea < 0 * bu::si::meter * bu::si::meter) {
    std::stringstream ss;
    ss << "Parameters tileWidth (" << tileWidth << "), tileHeight ("
       << tileHeight << "), minTileArea (" << minTileArea
       << ") must be positive.";
    errorString = ss.str();
    return false;
  }

  if (bbox.corners.outer().size() != 5) {
    bbox.corners.clear();
    minimalBoundingBox(area, bbox);
  }

  if (bbox.corners.outer().size() < 5)
    return false;
  double bboxWidth = bbox.width;
  double bboxHeight = bbox.height;
  FPoint origin = bbox.corners.outer()[0];

  // cout << "Origin: " << origin[0] << " " << origin[1] << endl;
  // Transform _mArea polygon to bounding box coordinate system.
  trans::rotate_transformer<boost::geometry::degree, double, 2, 2> rotate(
      bbox.angle * 180 / M_PI);
  trans::translate_transformer<double, 2, 2> translate(-origin.get<0>(),
                                                       -origin.get<1>());
  FPolygon translated_polygon;
  FPolygon rotated_polygon;
  boost::geometry::transform(area, translated_polygon, translate);
  boost::geometry::transform(translated_polygon, rotated_polygon, rotate);
  bg::correct(rotated_polygon);
  // cout << bg::wkt<BoostPolygon2D>(rotated_polygon) << endl;

  size_t iMax = ceil(bboxWidth / tileWidth.value());
  size_t jMax = ceil(bboxHeight / tileHeight.value());

  if (iMax < 1 || jMax < 1) {
    std::stringstream ss;
    ss << "Tile width (" << tileWidth << ") or tile height (" << tileHeight
       << ") to large for measurement area.";
    errorString = ss.str();
    return false;
  }

  trans::rotate_transformer<boost::geometry::degree, double, 2, 2> rotate_back(
      -bbox.angle * 180 / M_PI);
  trans::translate_transformer<double, 2, 2> translate_back(origin.get<0>(),
                                                            origin.get<1>());
  for (size_t i = 0; i < iMax; ++i) {
    double x_min = tileWidth.value() * i;
    double x_max = x_min + tileWidth.value();
    for (size_t j = 0; j < jMax; ++j) {
      double y_min = tileHeight.value() * j;
      double y_max = y_min + tileHeight.value();

      FPolygon tile_unclipped;
      tile_unclipped.outer().push_back(FPoint{x_min, y_min});
      tile_unclipped.outer().push_back(FPoint{x_min, y_max});
      tile_unclipped.outer().push_back(FPoint{x_max, y_max});
      tile_unclipped.outer().push_back(FPoint{x_max, y_min});
      tile_unclipped.outer().push_back(FPoint{x_min, y_min});

      std::deque<FPolygon> boost_tiles;
      if (!boost::geometry::intersection(tile_unclipped, rotated_polygon,
                                         boost_tiles))
        continue;

      for (FPolygon t : boost_tiles) {
        if (bg::area(t) > minTileArea.value()) {
          // Transform boost_tile to world coordinate system.
          FPolygon rotated_tile;
          FPolygon translated_tile;
          boost::geometry::transform(t, rotated_tile, rotate_back);
          boost::geometry::transform(rotated_tile, translated_tile,
                                     translate_back);

          // Store tile and calculate center point.
          tiles.push_back(translated_tile);
        }
      }
    }
  }

  if (tiles.size() < 1) {
    std::stringstream ss;
    ss << "No tiles calculated. Is the minTileArea (" << minTileArea
       << ") parameter large enough?";
    errorString = ss.str();
    return false;
  }

  return true;
}

bool joinedArea(const FPolygon &mArea, const FPolygon &sArea,
                const FPolygon &corridor, FPolygon &jArea,
                std::string &errorString) {
  // Measurement area and service area overlapping?
  bool overlapingSerMeas = bg::intersects(mArea, sArea) ? true : false;
  bool corridorValid = corridor.outer().size() > 0 ? true : false;

  // Check if corridor is connecting measurement area and service area.
  bool corridor_is_connection = false;
  if (corridorValid) {
    // Corridor overlaping with measurement area?
    if (bg::intersects(corridor, mArea)) {
      // Corridor overlaping with service area?
      if (bg::intersects(corridor, sArea)) {
        corridor_is_connection = true;
      }
    }
  }

  // Are areas joinable?
  std::deque<FPolygon> sol;
  FPolygon partialArea = mArea;
  if (overlapingSerMeas) {
    if (corridor_is_connection) {
      bg::union_(partialArea, corridor, sol);
    }
  } else if (corridor_is_connection) {
    bg::union_(partialArea, corridor, sol);
  } else {
    std::stringstream ss;
    auto printPoint = [&ss](const FPoint &p) {
      ss << " (" << p.get<0>() << ", " << p.get<1>() << ")";
    };
    ss << "Areas are not overlapping." << std::endl;
    ss << "Measurement area:";
    bg::for_each_point(mArea, printPoint);
    ss << std::endl;
    ss << "Service area:";
    bg::for_each_point(sArea, printPoint);
    ss << std::endl;
    ss << "Corridor:";
    bg::for_each_point(corridor, printPoint);
    ss << std::endl;
    errorString = ss.str();
    return false;
  }

  if (sol.size() > 0) {
    partialArea = sol[0];
    sol.clear();
  }

  // Join areas.
  bg::union_(partialArea, sArea, sol);

  if (sol.size() > 0) {
    jArea = sol[0];
  } else {
    std::stringstream ss;
    auto printPoint = [&ss](const FPoint &p) {
      ss << " (" << p.get<0>() << ", " << p.get<1>() << ")";
    };
    ss << "Areas not joinable." << std::endl;
    ss << "Measurement area:";
    bg::for_each_point(mArea, printPoint);
    ss << std::endl;
    ss << "Service area:";
    bg::for_each_point(sArea, printPoint);
    ss << std::endl;
    ss << "Corridor:";
    bg::for_each_point(corridor, printPoint);
    ss << std::endl;
    errorString = ss.str();
    return false;
  }

  return true;
}

bool joinedArea(const std::vector<FPolygon *> &areas, FPolygon &joinedArea) {
  if (areas.size() < 1)
    return false;
  joinedArea = *areas[0];
  std::deque<std::size_t> idxList;
  for (size_t i = 1; i < areas.size(); ++i)
    idxList.push_back(i);

  std::deque<FPolygon> sol;
  while (idxList.size() > 0) {
    bool success = false;
    for (auto it = idxList.begin(); it != idxList.end(); ++it) {
      bg::union_(joinedArea, *areas[*it], sol);
      if (sol.size() > 0) {
        joinedArea = sol[0];
        sol.clear();
        idxList.erase(it);
        success = true;
        break;
      }
    }
    if (!success)
      return false;
  }

  return true;
}

BoundingBox::BoundingBox() : width(0), height(0), angle(0) {}

void BoundingBox::clear() {
  width = 0;
  height = 0;
  angle = 0;
  corners.clear();
}

bool transectsFromScenario(Length distance, Length minLength, Angle angle,
                           const FPolygon &mArea,
                           const std::vector<FPolygon> &tiles,
                           const Progress &p, Transects &t,
                           string &errorString) {
  // Rotate measurement area by angle and calculate bounding box.
  FPolygon mAreaRotated;
  trans::rotate_transformer<bg::degree, double, 2, 2> rotate(angle.value() *
                                                             180 / M_PI);
  bg::transform(mArea, mAreaRotated, rotate);

  FBox box;
  boost::geometry::envelope(mAreaRotated, box);
  double x0 = box.min_corner().get<0>();
  double y0 = box.min_corner().get<1>();
  double x1 = box.max_corner().get<0>();
  double y1 = box.max_corner().get<1>();

  // Generate transects and convert them to clipper path.
  size_t num_t = int(ceil((y1 - y0) / distance.value())); // number of transects
  vector<ClipperLib::Path> transectsClipper;
  transectsClipper.reserve(num_t);
  for (size_t i = 0; i < num_t; ++i) {
    // calculate transect
    FPoint v1{x0, y0 + i * distance.value()};
    FPoint v2{x1, y0 + i * distance.value()};
    FLineString transect;
    transect.push_back(v1);
    transect.push_back(v2);
    // transform back
    FLineString temp_transect;
    trans::rotate_transformer<bg::degree, double, 2, 2> rotate_back(
        -angle.value() * 180 / M_PI);
    bg::transform(transect, temp_transect, rotate_back);
    // to clipper
    ClipperLib::IntPoint c1{static_cast<ClipperLib::cInt>(
                                temp_transect[0].get<0>() * CLIPPER_SCALE),
                            static_cast<ClipperLib::cInt>(
                                temp_transect[0].get<1>() * CLIPPER_SCALE)};
    ClipperLib::IntPoint c2{static_cast<ClipperLib::cInt>(
                                temp_transect[1].get<0>() * CLIPPER_SCALE),
                            static_cast<ClipperLib::cInt>(
                                temp_transect[1].get<1>() * CLIPPER_SCALE)};
    ClipperLib::Path path{c1, c2};
    transectsClipper.push_back(path);
  }

  if (transectsClipper.size() == 0) {
    std::stringstream ss;
    ss << "Not able to generate transects. Parameter: distance = " << distance
       << std::endl;
    errorString = ss.str();
    return false;
  }

  // Convert measurement area to clipper path.
  ClipperLib::Path mAreaClipper;
  for (auto vertex : mArea.outer()) {
    mAreaClipper.push_back(ClipperLib::IntPoint{
        static_cast<ClipperLib::cInt>(vertex.get<0>() * CLIPPER_SCALE),
        static_cast<ClipperLib::cInt>(vertex.get<1>() * CLIPPER_SCALE)});
  }

  // Perform clipping.
  // Clip transects to measurement area.
  ClipperLib::Clipper clipper;
  clipper.AddPath(mAreaClipper, ClipperLib::ptClip, true);
  clipper.AddPaths(transectsClipper, ClipperLib::ptSubject, false);
  ClipperLib::PolyTree clippedTransecs;
  clipper.Execute(ClipperLib::ctIntersection, clippedTransecs,
                  ClipperLib::pftNonZero, ClipperLib::pftNonZero);
  const auto *transects = &clippedTransecs;

  bool ignoreProgress = p.size() != tiles.size();
  ClipperLib::PolyTree clippedTransecs2;
  if (!ignoreProgress) {
    // Calculate processed tiles (_progress[i] == 100) and subtract them from
    // measurement area.
    size_t numTiles = p.size();
    vector<FPolygon> processedTiles;
    for (size_t i = 0; i < numTiles; ++i) {
      if (p[i] == 100) {
        processedTiles.push_back(tiles[i]);
      }
    }

    if (processedTiles.size() != numTiles) {
      vector<ClipperLib::Path> processedTilesClipper;
      for (const auto &t : processedTiles) {
        ClipperLib::Path path;
        for (const auto &vertex : t.outer()) {
          path.push_back(ClipperLib::IntPoint{
              static_cast<ClipperLib::cInt>(vertex.get<0>() * CLIPPER_SCALE),
              static_cast<ClipperLib::cInt>(vertex.get<1>() * CLIPPER_SCALE)});
        }
        processedTilesClipper.push_back(path);
      }

      // Subtract holes (tiles with measurement_progress == 100) from transects.
      clipper.Clear();
      for (const auto &child : clippedTransecs.Childs) {
        clipper.AddPath(child->Contour, ClipperLib::ptSubject, false);
      }
      clipper.AddPaths(processedTilesClipper, ClipperLib::ptClip, true);
      clipper.Execute(ClipperLib::ctDifference, clippedTransecs2,
                      ClipperLib::pftNonZero, ClipperLib::pftNonZero);
      transects = &clippedTransecs2;
    } else {
      // All tiles processed (t.size() not changed).
      return true;
    }
  }

  // Extract transects from  PolyTree and convert them to BoostLineString
  for (const auto &child : transects->Childs) {
    const auto &clipperTransect = child->Contour;
    FPoint v1{static_cast<double>(clipperTransect[0].X) / CLIPPER_SCALE,
              static_cast<double>(clipperTransect[0].Y) / CLIPPER_SCALE};
    FPoint v2{static_cast<double>(clipperTransect[1].X) / CLIPPER_SCALE,
              static_cast<double>(clipperTransect[1].Y) / CLIPPER_SCALE};

    FLineString transect{v1, v2};
    if (bg::length(transect) >= minLength.value()) {
      t.push_back(transect);
    }
  }

  if (t.size() == 0) {
    std::stringstream ss;
    ss << "Not able to generate transects. Parameter: minLength = " << minLength
       << std::endl;
    errorString = ss.str();
    return false;
  }
  return true;
}

bool route(const FPolygon &area, const Transects &transects,
           std::vector<Solution> &solutionVector, const RouteParameter &par) {

#ifdef SNAKE_SHOW_TIME
  auto start = std::chrono::high_resolution_clock::now();
#endif
  //================================================================
  // Create routing model.
  //================================================================
  // Use integer polygons to increase numerical robustness.
  // Convert area;
  IPolygon intArea;
  for (const auto &v : area.outer()) {
    auto p = float2Int(v);
    intArea.outer().push_back(p);
  }
  for (const auto &ring : area.inners()) {
    IRing intRing;
    for (const auto &v : ring) {
      auto p = float2Int(v);
      intRing.push_back(p);
    }
    intArea.inners().push_back(std::move(intRing));
  }

  // Helper classes.
  struct VirtualNode {
    VirtualNode(std::size_t f, std::size_t t) : fromIndex(f), toIndex(t) {}
    std::size_t fromIndex; // index for leaving node
    std::size_t toIndex;   // index for entering node
  };
  struct NodeToTransect {
    NodeToTransect(std::size_t i, bool r) : transectsIndex(i), reversed(r) {}
    std::size_t transectsIndex; // transects index
    bool reversed;              // transect reversed?
  };
  // Create vertex and node list
  std::vector<IPoint> vertices;
  std::vector<std::pair<std::size_t, std::size_t>> disjointNodes;
  std::vector<VirtualNode> nodeList;
  std::vector<NodeToTransect> nodeToTransectList;
  for (std::size_t i = 0; i < transects.size(); ++i) {
    const auto &t = transects[i];
    // Copy line edges only.
    if (t.size() == 1 || i == 0) {
      auto p = float2Int(t.back());
      vertices.push_back(p);
      nodeToTransectList.emplace_back(i, false);
      auto idx = vertices.size() - 1;
      nodeList.emplace_back(idx, idx);
    } else if (t.size() > 1) {
      auto p1 = float2Int(t.front());
      auto p2 = float2Int(t.back());
      vertices.push_back(p1);
      vertices.push_back(p2);
      nodeToTransectList.emplace_back(i, false);
      nodeToTransectList.emplace_back(i, true);
      auto fromIdx = vertices.size() - 1;
      auto toIdx = fromIdx - 1;
      nodeList.emplace_back(fromIdx, toIdx);
      nodeList.emplace_back(toIdx, fromIdx);
      disjointNodes.emplace_back(toIdx, fromIdx);
    } else { // transect empty
      std::cout << "ignoring empty transect with index " << i << std::endl;
    }
  }
#ifdef SNAKE_DEBUG
  // Print.
  std::cout << "nodeToTransectList:" << std::endl;
  std::cout << "node:transectIndex:reversed" << std::endl;
  std::size_t c = 0;
  for (const auto &n2t : nodeToTransectList) {
    std::cout << c++ << ":" << n2t.transectsIndex << ":" << n2t.reversed
              << std::endl;
  }
  std::cout << "nodeList:" << std::endl;
  std::cout << "node:fromIndex:toIndex" << std::endl;
  c = 0;
  for (const auto &n : nodeList) {
    std::cout << c++ << ":" << n.fromIndex << ":" << n.toIndex << std::endl;
  }
  std::cout << "disjoint nodes:" << std::endl;
  std::cout << "number:nodes" << std::endl;
  c = 0;
  for (const auto &d : disjointNodes) {
    std::cout << c++ << ":" << d.first << "," << d.second << std::endl;
  }
#endif

  // Add polygon vertices.
  for (auto &v : intArea.outer()) {
    vertices.push_back(v);
  }
  for (auto &ring : intArea.inners()) {
    for (auto &v : ring) {
      vertices.push_back(v);
    }
  }

  // Create connection graph (inf == no connection between vertices).
  // Note: graph is not symmetric.
  auto n = vertices.size();
  // Matrix must be double since integers don't have infinity and nan
  Matrix<double> connectionGraph(n, n);
  for (std::size_t i = 0; i < n; ++i) {
    auto &fromVertex = vertices[i];
    for (std::size_t j = 0; j < n; ++j) {
      auto &toVertex = vertices[j];
      ILineString line{fromVertex, toVertex};
      if (bg::covered_by(line, intArea)) {
        connectionGraph(i, j) = bg::length(line);
      } else {
        connectionGraph(i, j) = std::numeric_limits<double>::infinity();
      }
    }
  }
#ifdef SNAKE_DEBUG
  std::cout << "connection grah:" << std::endl;
  std::cout << connectionGraph << std::endl;
#endif

  // Create distance matrix.
  auto distLambda = [&connectionGraph](std::size_t i, std::size_t j) -> double {
    return connectionGraph(i, j);
  };
  auto nNodes = nodeList.size();
  Matrix<IntType> distanceMatrix(nNodes, nNodes);
  for (std::size_t i = 0; i < nNodes; ++i) {
    distanceMatrix(i, i) = 0;
    for (std::size_t j = i + 1; j < nNodes; ++j) {
      auto dist = connectionGraph(i, j);
      if (std::isinf(dist)) {
        std::vector<std::size_t> route;
        if (!dijkstraAlgorithm(n, i, j, route, dist, distLambda)) {
          std::stringstream ss;
          ss << "Distance matrix calculation failed. connection graph: "
             << connectionGraph << std::endl;
          ss << "area: " << bg::wkt(area) << std::endl;
          ss << "transects:" << std::endl;
          for (const auto &t : transects) {

            ss << bg::wkt(t) << std::endl;
          }

          par.errorString = ss.str();
          return false;
        }
        (void)route;
      }
      distanceMatrix(i, j) = dist;
      distanceMatrix(j, i) = dist;
    }
  }
#ifdef SNAKE_DEBUG
  std::cout << "distance matrix:" << std::endl;
  std::cout << distanceMatrix << std::endl;
#endif

  // Create (asymmetric) routing matrix.
  Matrix<IntType> routingMatrix(nNodes, nNodes);
  for (std::size_t i = 0; i < nNodes; ++i) {
    auto fromNode = nodeList[i];
    for (std::size_t j = 0; j < nNodes; ++j) {
      auto toNode = nodeList[j];
      routingMatrix(i, j) = distanceMatrix(fromNode.fromIndex, toNode.toIndex);
    }
  }
  // Insert max for disjoint nodes.
  for (const auto &d : disjointNodes) {
    auto i = d.first;
    auto j = d.second;
    routingMatrix(i, j) = std::numeric_limits<IntType>::max();
    routingMatrix(j, i) = std::numeric_limits<IntType>::max();
  }
#ifdef SNAKE_DEBUG
  std::cout << "routing matrix:" << std::endl;
  std::cout << routingMatrix << std::endl;
#endif

  // Create Routing Index Manager.
  auto minNumTransectsPerRun =
      std::max<std::size_t>(1, par.minNumTransectsPerRun);
  auto maxRuns = std::max<std::size_t>(
      1, std::floor(double(transects.size()) / minNumTransectsPerRun));
  auto numRuns = std::max<std::size_t>(1, par.numRuns);
  numRuns = std::min<std::size_t>(numRuns, maxRuns);

  RoutingIndexManager::NodeIndex depot(0);
  //  std::vector<RoutingIndexManager::NodeIndex> depots(numRuns, depot);
  //  RoutingIndexManager manager(nNodes, numRuns, depots, depots);
  RoutingIndexManager manager(nNodes, numRuns, depot);

  // Create Routing Model.
  RoutingModel routing(manager);
  // Create and register a transit callback.
  const int transitCallbackIndex = routing.RegisterTransitCallback(
      [&routingMatrix, &manager](int64 from_index, int64 to_index) -> int64 {
        // Convert from routing variable Index to distance matrix NodeIndex.
        auto from_node = manager.IndexToNode(from_index).value();
        auto to_node = manager.IndexToNode(to_index).value();
        return routingMatrix(from_node, to_node);
      });
  // Define cost of each arc.
  routing.SetArcCostEvaluatorOfAllVehicles(transitCallbackIndex);
  // Add distance dimension.
  if (numRuns > 1) {
    routing.AddDimension(transitCallbackIndex, 0, 300000000,
                         true, // start cumul to zero
                         "Distance");
    routing.GetMutableDimension("Distance")
        ->SetGlobalSpanCostCoefficient(100000000);
  }

  // Define disjunctions.
#ifdef SNAKE_DEBUG
  std::cout << "disjunctions:" << std::endl;
#endif
  for (const auto &d : disjointNodes) {
    auto i = d.first;
    auto j = d.second;
#ifdef SNAKE_DEBUG
    std::cout << i << "," << j << std::endl;
#endif
    auto idx0 = manager.NodeToIndex(RoutingIndexManager::NodeIndex(i));
    auto idx1 = manager.NodeToIndex(RoutingIndexManager::NodeIndex(j));
    std::vector<int64> disj{idx0, idx1};
    routing.AddDisjunction(disj, -1 /*force cardinality*/, 1 /*cardinality*/);
  }

  // Set first solution heuristic.
  auto searchParameters = DefaultRoutingSearchParameters();
  searchParameters.set_first_solution_strategy(
      FirstSolutionStrategy::PATH_CHEAPEST_ARC);
  // Number of solutions.
  auto numSolutionsPerRun = std::max<std::size_t>(1, par.numSolutionsPerRun);
  searchParameters.set_number_of_solutions_to_collect(numSolutionsPerRun);
  // Set costume limit.
  auto *solver = routing.solver();
  auto *limit = solver->MakeCustomLimit(par.stop);
  routing.AddSearchMonitor(limit);
#ifdef SNAKE_SHOW_TIME
  auto delta = std::chrono::duration_cast<std::chrono::milliseconds>(
      std::chrono::high_resolution_clock::now() - start);
  cout << "create routing model: " << delta.count() << " ms" << endl;
#endif

  //================================================================
  // Solve model.
  //================================================================
#ifdef SNAKE_SHOW_TIME
  start = std::chrono::high_resolution_clock::now();
#endif
  auto pSolutions = std::make_unique<std::vector<const Assignment *>>();
  (void)routing.SolveWithParameters(searchParameters, pSolutions.get());
#ifdef SNAKE_SHOW_TIME
  delta = std::chrono::duration_cast<std::chrono::milliseconds>(
      std::chrono::high_resolution_clock::now() - start);
  cout << "solve routing model: " << delta.count() << " ms" << endl;
#endif
  if (par.stop()) {
    par.errorString = "User terminated.";
    return false;
  }
  if (pSolutions->size() == 0) {
    std::stringstream ss;
    ss << "No solution found." << std::endl;
    par.errorString = ss.str();
    return false;
  }

  //================================================================
  // Construc route.
  //================================================================
#ifdef SNAKE_SHOW_TIME
  start = std::chrono::high_resolution_clock::now();
#endif
  long long counter = -1;
  // Note: route number 0 corresponds to the best route which is the last entry
  // of *pSolutions.
  for (auto solution = pSolutions->end() - 1; solution >= pSolutions->begin();
       --solution) {
    ++counter;
    if (!*solution || (*solution)->Size() <= 1) {
      std::stringstream ss;
      ss << par.errorString << "Solution " << counter << "invalid."
         << std::endl;
      par.errorString = ss.str();
      continue;
    }
    // Iterate over all routes.
    Solution routeVector;
    for (std::size_t vehicle = 0; vehicle < numRuns; ++vehicle) {
      if (!routing.IsVehicleUsed(**solution, vehicle))
        continue;

      // Create index list.
      auto index = routing.Start(vehicle);
      std::vector<size_t> route_idx;
      route_idx.push_back(manager.IndexToNode(index).value());
      while (!routing.IsEnd(index)) {
        index = (*solution)->Value(routing.NextVar(index));
        route_idx.push_back(manager.IndexToNode(index).value());
      }

#ifdef SNAKE_DEBUG
      // Print route.
      std::cout << "route " << counter
                << " route_idx.size() = " << route_idx.size() << std::endl;
      std::cout << "route: ";
      for (const auto &idx : route_idx) {
        std::cout << idx << ", ";
      }
      std::cout << std::endl;
#endif
      if (route_idx.size() < 2) {
        std::stringstream ss;
        ss << par.errorString
           << "Error while assembling route (solution = " << counter
           << ", run = " << vehicle << ")." << std::endl;
        par.errorString = ss.str();
        continue;
      }

      // Assemble route.
      Route r;
      auto &path = r.path;
      auto &info = r.info;
      for (size_t i = 0; i < route_idx.size() - 1; ++i) {
        size_t nodeIndex0 = route_idx[i];
        size_t nodeIndex1 = route_idx[i + 1];
        const auto &n2t0 = nodeToTransectList[nodeIndex0];
        info.emplace_back(n2t0.transectsIndex, n2t0.reversed);
        // Copy transect to route.
        const auto &t = transects[n2t0.transectsIndex];
        if (n2t0.reversed) { // transect reversal needed?
          for (auto it = t.end() - 1; it > t.begin(); --it) {
            path.push_back(*it);
          }
        } else {
          for (auto it = t.begin(); it < t.end() - 1; ++it) {
            path.push_back(*it);
          }
        }
        // Connect transects.
        std::vector<size_t> idxList;
        if (!shortestPathFromGraph(connectionGraph,
                                   nodeList[nodeIndex0].fromIndex,
                                   nodeList[nodeIndex1].toIndex, idxList)) {
          std::stringstream ss;
          ss << par.errorString
             << "Error while assembling route (solution = " << counter
             << ", run = " << vehicle << ")." << std::endl;
          par.errorString = ss.str();
          continue;
        }
        if (i != route_idx.size() - 2) {
          idxList.pop_back();
        }
        for (auto idx : idxList) {
          auto p = int2Float(vertices[idx]);
          path.push_back(p);
        }
      }
      // Append last transect info.
      const auto &n2t0 = nodeToTransectList.back();
      info.emplace_back(n2t0.transectsIndex, n2t0.reversed);

      if (path.size() < 2 || info.size() < 2) {
        std::stringstream ss;
        ss << par.errorString << "Route empty (solution = " << counter
           << ", run = " << vehicle << ")." << std::endl;
        par.errorString = ss.str();