From 559237c820b112913071683043885297eafd6ec1 Mon Sep 17 00:00:00 2001
From: Willem Sonke <willem@wimiso.nl>
Date: Tue, 22 Aug 2023 11:52:26 +0200
Subject: [PATCH 01/27] Add optimization demo application

---
 demos/flow_map/CMakeLists.txt                 |   1 +
 .../flow_map/optimization_demo/CMakeLists.txt |  19 +++
 .../optimization_demo/optimization_demo.cpp   | 109 ++++++++++++++++++
 .../optimization_demo/optimization_demo.h     |  50 ++++++++
 .../spiral_tree_demo/spiral_tree_demo.cpp     |  13 ---
 5 files changed, 179 insertions(+), 13 deletions(-)
 create mode 100644 demos/flow_map/optimization_demo/CMakeLists.txt
 create mode 100644 demos/flow_map/optimization_demo/optimization_demo.cpp
 create mode 100644 demos/flow_map/optimization_demo/optimization_demo.h

diff --git a/demos/flow_map/CMakeLists.txt b/demos/flow_map/CMakeLists.txt
index c14886126..135f5f0a7 100644
--- a/demos/flow_map/CMakeLists.txt
+++ b/demos/flow_map/CMakeLists.txt
@@ -1 +1,2 @@
+add_subdirectory(optimization_demo)
 add_subdirectory(spiral_tree_demo)
diff --git a/demos/flow_map/optimization_demo/CMakeLists.txt b/demos/flow_map/optimization_demo/CMakeLists.txt
new file mode 100644
index 000000000..1c9b3ba80
--- /dev/null
+++ b/demos/flow_map/optimization_demo/CMakeLists.txt
@@ -0,0 +1,19 @@
+set(SOURCES
+    optimization_demo.cpp
+)
+
+add_executable(optimization_demo ${SOURCES})
+
+target_link_libraries(
+    optimization_demo
+    PRIVATE
+    ${COMMON_CLA_TARGET}
+    core
+    flow_map
+    renderer
+    CGAL::CGAL
+    Qt5::Widgets
+    glog::glog
+)
+
+install(TARGETS optimization_demo DESTINATION ${INSTALL_BINARY_DIR})
diff --git a/demos/flow_map/optimization_demo/optimization_demo.cpp b/demos/flow_map/optimization_demo/optimization_demo.cpp
new file mode 100644
index 000000000..5fec0089e
--- /dev/null
+++ b/demos/flow_map/optimization_demo/optimization_demo.cpp
@@ -0,0 +1,109 @@
+/*
+The Necklace Map console application implements the algorithmic
+geo-visualization method by the same name, developed by
+Bettina Speckmann and Kevin Verbeek at TU Eindhoven
+(DOI: 10.1109/TVCG.2010.180 & 10.1142/S021819591550003X).
+Copyright (C) 2021  Netherlands eScience Center and TU Eindhoven
+
+This program is free software: you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation, either version 3 of the License, or
+(at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with this program.  If not, see <https://www.gnu.org/licenses/>.
+*/
+
+#include "optimization_demo.h"
+
+#include <QApplication>
+#include <QCheckBox>
+
+#include "cartocrow/core/core.h"
+#include "cartocrow/flow_map/painting.h"
+#include "cartocrow/flow_map/parameters.h"
+#include "cartocrow/flow_map/place.h"
+#include "cartocrow/flow_map/reachable_region_algorithm.h"
+#include "cartocrow/flow_map/spiral_tree.h"
+#include "cartocrow/flow_map/spiral_tree_obstructed_algorithm.h"
+#include "cartocrow/flow_map/spiral_tree_unobstructed_algorithm.h"
+#include "cartocrow/renderer/geometry_painting.h"
+#include "cartocrow/renderer/geometry_widget.h"
+
+using namespace cartocrow;
+using namespace cartocrow::flow_map;
+using namespace cartocrow::renderer;
+
+OptimizationDemo::OptimizationDemo() {
+	setWindowTitle("CartoCrow – Optimization demo");
+
+	m_places.push_back(std::make_shared<Point<Inexact>>(11.2121212, 17.0707070));
+	m_places.push_back(std::make_shared<Point<Inexact>>(13.9393939, -14.1414141));
+	m_places.push_back(std::make_shared<Point<Inexact>>(-4.5454545, -18.9898989));
+	m_places.push_back(std::make_shared<Point<Inexact>>(16.6666666, 6.1616161));
+	m_places.push_back(std::make_shared<Point<Inexact>>(-9.8989898, 13.9393939));
+	m_places.push_back(std::make_shared<Point<Inexact>>(-16.1616161, -2.6262626));
+
+	m_renderer = new GeometryWidget();
+	m_renderer->setMaxZoom(10000);
+	m_renderer->setGridMode(GeometryWidget::GridMode::POLAR);
+	setCentralWidget(m_renderer);
+
+	QToolBar* toolBar = new QToolBar();
+	toolBar->addSeparator();
+	toolBar->addWidget(new QLabel("α = "));
+	m_alphaSlider = new QSlider(Qt::Horizontal);
+	m_alphaSlider->setMinimum(0);
+	m_alphaSlider->setMaximum(499);
+	m_alphaSlider->setValue(139);
+	toolBar->addWidget(m_alphaSlider);
+	addToolBar(toolBar);
+	m_alphaLabel = new QLabel("0.139π");
+	toolBar->addWidget(m_alphaLabel);
+	connect(m_alphaSlider, &QSlider::valueChanged, [&](int value) {
+		m_alpha = M_PI * value / 1000.0;
+		m_alphaLabel->setText(QString("%1π").arg(value / 1000.0, 0, 'f', 3));
+		recalculate();
+	});
+	for (auto place : m_places) {
+		m_renderer->registerEditable(place);
+	}
+	connect(m_renderer, &GeometryWidget::edited, [&]() {
+		recalculate();
+	});
+	recalculate();
+}
+
+void OptimizationDemo::recalculate() {
+	auto tree = std::make_shared<SpiralTree>(Point<Inexact>(0, 0), m_alpha);
+	for (const auto& place : m_places) {
+		tree->addPlace("", *place, 1);
+	}
+	tree->addShields();
+
+	m_renderer->clear();
+	ReachableRegionAlgorithm reachableRegionAlg(tree);
+	ReachableRegionAlgorithm::ReachableRegion reachableRegion = reachableRegionAlg.run();
+
+	SpiralTreeObstructedAlgorithm spiralTreeAlg(tree, reachableRegion);
+	spiralTreeAlg.run();
+
+	Painting::Options options;
+	Painting p(nullptr, tree, options);
+	auto painting = std::make_shared<Painting>(nullptr, tree, options);
+	m_renderer->addPainting(painting, "Spiral tree");
+
+	m_renderer->update();
+}
+
+int main(int argc, char* argv[]) {
+	QApplication app(argc, argv);
+	OptimizationDemo demo;
+	demo.show();
+	app.exec();
+}
diff --git a/demos/flow_map/optimization_demo/optimization_demo.h b/demos/flow_map/optimization_demo/optimization_demo.h
new file mode 100644
index 000000000..8b331906e
--- /dev/null
+++ b/demos/flow_map/optimization_demo/optimization_demo.h
@@ -0,0 +1,50 @@
+/*
+The Necklace Map console application implements the algorithmic
+geo-visualization method by the same name, developed by
+Bettina Speckmann and Kevin Verbeek at TU Eindhoven
+(DOI: 10.1109/TVCG.2010.180 & 10.1142/S021819591550003X).
+Copyright (C) 2021  Netherlands eScience Center and TU Eindhoven
+
+This program is free software: you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation, either version 3 of the License, or
+(at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with this program.  If not, see <https://www.gnu.org/licenses/>.
+*/
+
+#include <QLabel>
+#include <QMainWindow>
+#include <QSlider>
+
+#include <memory>
+#include <vector>
+
+#include <cartocrow/core/core.h>
+#include <cartocrow/renderer/geometry_widget.h>
+
+using namespace cartocrow;
+using namespace cartocrow::renderer;
+
+class OptimizationDemo : public QMainWindow {
+	Q_OBJECT;
+
+  public:
+	OptimizationDemo();
+
+  private:
+	void recalculate();
+	Number<Inexact> m_alpha = 25 * M_PI / 180;
+
+	std::vector<std::shared_ptr<Point<Inexact>>> m_places;
+
+	GeometryWidget* m_renderer;
+	QSlider* m_alphaSlider;
+	QLabel* m_alphaLabel;
+};
diff --git a/demos/flow_map/spiral_tree_demo/spiral_tree_demo.cpp b/demos/flow_map/spiral_tree_demo/spiral_tree_demo.cpp
index 8abe10ac0..ef8907535 100644
--- a/demos/flow_map/spiral_tree_demo/spiral_tree_demo.cpp
+++ b/demos/flow_map/spiral_tree_demo/spiral_tree_demo.cpp
@@ -201,19 +201,6 @@ void SpiralTreeDemo::recalculate() {
 	m_renderer->update();
 }
 
-/*Point<Inexact>* SpiralTreeDemo::findClosestPoint(Point<Inexact> p, Number<Inexact> radius) {
-	Point<Inexact>* closest = nullptr;
-	Number<Inexact> minSquaredDistance = radius * radius;
-	for (auto& vertex : m_obstacle) {
-		Number<Inexact> squaredDistance = (vertex - p).squared_length();
-		if (squaredDistance < minSquaredDistance) {
-			minSquaredDistance = squaredDistance;
-			closest = &vertex;
-		}
-	}
-	return closest;
-}*/
-
 int main(int argc, char* argv[]) {
 	QApplication app(argc, argv);
 	SpiralTreeDemo demo;

From 4d1712f255a3c617d062a48ca0b01851da3a0409 Mon Sep 17 00:00:00 2001
From: Willem Sonke <willem@wimiso.nl>
Date: Tue, 22 Aug 2023 18:18:43 +0200
Subject: [PATCH 02/27] Set up data structure for the smoothed tree

---
 cartocrow/flow_map/CMakeLists.txt             |  4 ++
 cartocrow/flow_map/smooth_tree.cpp            | 33 +++++++++++
 cartocrow/flow_map/smooth_tree.h              | 55 +++++++++++++++++++
 cartocrow/flow_map/smooth_tree_painting.cpp   | 39 +++++++++++++
 cartocrow/flow_map/smooth_tree_painting.h     | 38 +++++++++++++
 .../optimization_demo/optimization_demo.cpp   | 18 +++---
 .../spiral_tree_demo/spiral_tree_demo.cpp     |  1 -
 7 files changed, 177 insertions(+), 11 deletions(-)
 create mode 100644 cartocrow/flow_map/smooth_tree.cpp
 create mode 100644 cartocrow/flow_map/smooth_tree.h
 create mode 100644 cartocrow/flow_map/smooth_tree_painting.cpp
 create mode 100644 cartocrow/flow_map/smooth_tree_painting.h

diff --git a/cartocrow/flow_map/CMakeLists.txt b/cartocrow/flow_map/CMakeLists.txt
index 25466ee47..f9b5fd538 100644
--- a/cartocrow/flow_map/CMakeLists.txt
+++ b/cartocrow/flow_map/CMakeLists.txt
@@ -8,6 +8,8 @@ set(SOURCES
 	polar_point.cpp
 	polar_segment.cpp
 	reachable_region_algorithm.cpp
+	smooth_tree.cpp
+	smooth_tree_painting.cpp
 	spiral.cpp
 	spiral_segment.cpp
 	spiral_tree.cpp
@@ -28,6 +30,8 @@ set(HEADERS
 	polar_point.h
 	polar_segment.h
 	reachable_region_algorithm.h
+	smooth_tree.h
+	smooth_tree_painting.h
 	spiral.h
 	spiral_segment.h
 	spiral_tree.h
diff --git a/cartocrow/flow_map/smooth_tree.cpp b/cartocrow/flow_map/smooth_tree.cpp
new file mode 100644
index 000000000..7dfc81090
--- /dev/null
+++ b/cartocrow/flow_map/smooth_tree.cpp
@@ -0,0 +1,33 @@
+/*
+The Flow Map library implements the algorithmic geo-visualization
+method by the same name, developed by Kevin Verbeek, Kevin Buchin,
+and Bettina Speckmann at TU Eindhoven
+(DOI: 10.1007/s00453-013-9867-z & 10.1109/TVCG.2011.202).
+Copyright (C) 2021  Netherlands eScience Center and TU Eindhoven
+
+This program is free software: you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation, either version 3 of the License, or
+(at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with this program.  If not, see <https://www.gnu.org/licenses/>.
+*/
+
+#include "smooth_tree.h"
+#include "cartocrow/flow_map/smooth_tree_painting.h"
+
+namespace cartocrow::flow_map {
+
+SmoothTree::SmoothTree(const std::shared_ptr<SpiralTree> spiralTree) : m_tree(spiralTree) {}
+
+const std::vector<std::shared_ptr<Node>>& SmoothTree::nodes() const {
+	return m_nodes;
+}
+
+} // namespace cartocrow::flow_map
diff --git a/cartocrow/flow_map/smooth_tree.h b/cartocrow/flow_map/smooth_tree.h
new file mode 100644
index 000000000..d5d901e08
--- /dev/null
+++ b/cartocrow/flow_map/smooth_tree.h
@@ -0,0 +1,55 @@
+/*
+The Flow Map library implements the algorithmic geo-visualization
+method by the same name, developed by Kevin Verbeek, Kevin Buchin,
+and Bettina Speckmann at TU Eindhoven
+(DOI: 10.1007/s00453-013-9867-z & 10.1109/TVCG.2011.202).
+Copyright (C) 2021  Netherlands eScience Center and TU Eindhoven
+
+This program is free software: you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation, either version 3 of the License, or
+(at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with this program.  If not, see <https://www.gnu.org/licenses/>.
+*/
+
+#ifndef CARTOCROW_FLOW_MAP_SMOOTH_TREE_H
+#define CARTOCROW_FLOW_MAP_SMOOTH_TREE_H
+
+#include <deque>
+#include <map>
+#include <memory>
+#include <queue>
+#include <vector>
+
+#include "spiral_tree.h"
+
+namespace cartocrow::flow_map {
+
+/// A smoothed tree.
+class SmoothTree {
+
+  public:
+	/// Constructs a smooth tree.
+	SmoothTree(const std::shared_ptr<SpiralTree> spiralTree);
+
+	/// Returns a list of the nodes in this smooth tree.
+	const std::vector<std::shared_ptr<Node>>& nodes() const;
+
+  private:
+	/// The spiral tree underlying this smooth tree.
+	std::shared_ptr<SpiralTree> m_tree;
+
+	/// List of nodes in this tree.
+	std::vector<std::shared_ptr<Node>> m_nodes;
+};
+
+} // namespace cartocrow::flow_map
+
+#endif
diff --git a/cartocrow/flow_map/smooth_tree_painting.cpp b/cartocrow/flow_map/smooth_tree_painting.cpp
new file mode 100644
index 000000000..5d61cd303
--- /dev/null
+++ b/cartocrow/flow_map/smooth_tree_painting.cpp
@@ -0,0 +1,39 @@
+#include "smooth_tree_painting.h"
+
+#include "../core/core.h"
+#include "../renderer/geometry_renderer.h"
+#include <CGAL/Origin.h>
+
+namespace cartocrow::flow_map {
+
+SmoothTreePainting::SmoothTreePainting(std::shared_ptr<SmoothTree> tree, const Options options)
+    : m_tree(tree), m_options(options) {}
+
+void SmoothTreePainting::paint(renderer::GeometryRenderer& renderer) const {
+	paintFlow(renderer);
+	paintNodes(renderer);
+}
+
+void SmoothTreePainting::paintFlow(renderer::GeometryRenderer& renderer) const {
+	renderer.setMode(renderer::GeometryRenderer::stroke);
+	renderer.setStroke(Color{255, 84, 32}, 4);
+	for (const auto& node : m_tree->nodes()) {
+		if (node->m_parent == nullptr) {
+			continue;
+		}
+		renderer.draw(Segment<Inexact>(node->m_parent->m_position.toCartesian(),
+		                               node->m_position.toCartesian()));
+	}
+}
+
+void SmoothTreePainting::paintNodes(renderer::GeometryRenderer& renderer) const {
+	renderer.setMode(renderer::GeometryRenderer::vertices);
+	renderer.setStroke(Color{100, 100, 100}, 4);
+	for (const auto& node : m_tree->nodes()) {
+		if (!node->isSteiner()) {
+			renderer.draw(node->m_position.toCartesian());
+		}
+	}
+}
+
+} // namespace cartocrow::flow_map
diff --git a/cartocrow/flow_map/smooth_tree_painting.h b/cartocrow/flow_map/smooth_tree_painting.h
new file mode 100644
index 000000000..ecb121f91
--- /dev/null
+++ b/cartocrow/flow_map/smooth_tree_painting.h
@@ -0,0 +1,38 @@
+#ifndef CARTOCROW_FLOW_MAP_SMOOTH_TREE_PAINTING
+#define CARTOCROW_FLOW_MAP_SMOOTH_TREE_PAINTING
+
+#include "../renderer/geometry_painting.h"
+#include "../renderer/geometry_renderer.h"
+#include "spiral.h"
+
+#include "smooth_tree.h"
+
+namespace cartocrow::flow_map {
+
+/// The \ref renderer::GeometryPainting "GeometryPainting" for a \ref SmoothTree.
+class SmoothTreePainting : public renderer::GeometryPainting {
+
+  public:
+	/// Options that determine what to draw in the painting.
+	struct Options {
+		double spiralStep = 0.01;
+		double spiralMax = 6.0;
+	};
+
+	/// Creates a new painting with the given map and tree.
+	SmoothTreePainting(std::shared_ptr<SmoothTree> tree, const Options options);
+
+  protected:
+	void paint(renderer::GeometryRenderer& renderer) const override;
+
+  private:
+	void paintNodes(renderer::GeometryRenderer& renderer) const;
+	void paintFlow(renderer::GeometryRenderer& renderer) const;
+
+	std::shared_ptr<SmoothTree> m_tree;
+	Options m_options;
+};
+
+} // namespace cartocrow::flow_map
+
+#endif //CARTOCROW_FLOW_MAP_PAINTING
diff --git a/demos/flow_map/optimization_demo/optimization_demo.cpp b/demos/flow_map/optimization_demo/optimization_demo.cpp
index 5fec0089e..414020289 100644
--- a/demos/flow_map/optimization_demo/optimization_demo.cpp
+++ b/demos/flow_map/optimization_demo/optimization_demo.cpp
@@ -29,6 +29,8 @@ along with this program.  If not, see <https://www.gnu.org/licenses/>.
 #include "cartocrow/flow_map/parameters.h"
 #include "cartocrow/flow_map/place.h"
 #include "cartocrow/flow_map/reachable_region_algorithm.h"
+#include "cartocrow/flow_map/smooth_tree.h"
+#include "cartocrow/flow_map/smooth_tree_painting.h"
 #include "cartocrow/flow_map/spiral_tree.h"
 #include "cartocrow/flow_map/spiral_tree_obstructed_algorithm.h"
 #include "cartocrow/flow_map/spiral_tree_unobstructed_algorithm.h"
@@ -87,16 +89,12 @@ void OptimizationDemo::recalculate() {
 	tree->addShields();
 
 	m_renderer->clear();
-	ReachableRegionAlgorithm reachableRegionAlg(tree);
-	ReachableRegionAlgorithm::ReachableRegion reachableRegion = reachableRegionAlg.run();
-
-	SpiralTreeObstructedAlgorithm spiralTreeAlg(tree, reachableRegion);
-	spiralTreeAlg.run();
-
-	Painting::Options options;
-	Painting p(nullptr, tree, options);
-	auto painting = std::make_shared<Painting>(nullptr, tree, options);
-	m_renderer->addPainting(painting, "Spiral tree");
+	ReachableRegionAlgorithm::ReachableRegion reachableRegion = ReachableRegionAlgorithm(tree).run();
+	SpiralTreeObstructedAlgorithm(tree, reachableRegion).run();
+	auto smoothTree = std::make_shared<SmoothTree>(tree);
+	SmoothTreePainting::Options options;
+	auto painting = std::make_shared<SmoothTreePainting>(smoothTree, options);
+	m_renderer->addPainting(painting, "Smooth tree");
 
 	m_renderer->update();
 }
diff --git a/demos/flow_map/spiral_tree_demo/spiral_tree_demo.cpp b/demos/flow_map/spiral_tree_demo/spiral_tree_demo.cpp
index ef8907535..e9b349040 100644
--- a/demos/flow_map/spiral_tree_demo/spiral_tree_demo.cpp
+++ b/demos/flow_map/spiral_tree_demo/spiral_tree_demo.cpp
@@ -194,7 +194,6 @@ void SpiralTreeDemo::recalculate() {
 	t.output();
 
 	Painting::Options options;
-	Painting p(nullptr, tree, options);
 	auto painting = std::make_shared<Painting>(nullptr, tree, options);
 	m_renderer->addPainting(painting, "Spiral tree");
 

From 51165fd2239dd52fa185f877063ea29412c5c529 Mon Sep 17 00:00:00 2001
From: Willem Sonke <willem@wimiso.nl>
Date: Wed, 23 Aug 2023 17:03:53 +0200
Subject: [PATCH 03/27] Implement edge subdivision for the smooth tree

---
 cartocrow/flow_map/smooth_tree.cpp          | 37 ++++++++++++++++++++-
 cartocrow/flow_map/smooth_tree.h            |  2 ++
 cartocrow/flow_map/smooth_tree_painting.cpp |  4 +--
 cartocrow/flow_map/spiral_tree.cpp          |  4 +++
 cartocrow/flow_map/spiral_tree.h            |  2 ++
 5 files changed, 45 insertions(+), 4 deletions(-)

diff --git a/cartocrow/flow_map/smooth_tree.cpp b/cartocrow/flow_map/smooth_tree.cpp
index 7dfc81090..e33f6a1b0 100644
--- a/cartocrow/flow_map/smooth_tree.cpp
+++ b/cartocrow/flow_map/smooth_tree.cpp
@@ -24,7 +24,42 @@ along with this program.  If not, see <https://www.gnu.org/licenses/>.
 
 namespace cartocrow::flow_map {
 
-SmoothTree::SmoothTree(const std::shared_ptr<SpiralTree> spiralTree) : m_tree(spiralTree) {}
+SmoothTree::SmoothTree(const std::shared_ptr<SpiralTree> spiralTree) : m_tree(spiralTree) {
+	Number<Inexact> rMax = 0;
+	for (const auto& node : m_tree->nodes()) {
+		if (node->m_position.r() > rMax) {
+			rMax = node->m_position.r();
+		}
+	}
+	constructSmoothTree(spiralTree->root(), rMax / 10); // TODO make 10 configurable
+}
+
+std::shared_ptr<Node> SmoothTree::constructSmoothTree(const std::shared_ptr<Node>& node, Number<Inexact> maxRStep) {
+	auto smoothNode = std::make_shared<Node>(node->m_position);
+	m_nodes.push_back(smoothNode);
+	for (const auto& child : node->m_children) {
+		auto smoothChild = constructSmoothTree(child, maxRStep);
+		const Spiral spiral(node->m_position, child->m_position);
+		Number<Inexact> rMin = node->m_position.r();
+		Number<Inexact> rMax = child->m_position.r();
+		Number<Inexact> rRange = rMax - rMin;
+		int subdivisionCount = std::ceil(rRange / maxRStep);
+		auto previous = smoothNode;
+		for (double i = 1; i < subdivisionCount; ++i) {
+			const PolarPoint position =
+			    spiral.evaluate(spiral.parameterForR(rMin + i * (rMax - rMin) / subdivisionCount));
+			auto subdivision = std::make_shared<Node>(position);
+			m_nodes.push_back(subdivision);
+			subdivision->m_parent = previous;
+			previous->m_children.push_back(subdivision);
+			previous = subdivision;
+		}
+		smoothChild->m_parent = previous;
+		previous->m_children.push_back(smoothChild->m_parent);
+	}
+
+	return smoothNode;
+}
 
 const std::vector<std::shared_ptr<Node>>& SmoothTree::nodes() const {
 	return m_nodes;
diff --git a/cartocrow/flow_map/smooth_tree.h b/cartocrow/flow_map/smooth_tree.h
index d5d901e08..12a87803b 100644
--- a/cartocrow/flow_map/smooth_tree.h
+++ b/cartocrow/flow_map/smooth_tree.h
@@ -48,6 +48,8 @@ class SmoothTree {
 
 	/// List of nodes in this tree.
 	std::vector<std::shared_ptr<Node>> m_nodes;
+
+	std::shared_ptr<Node> constructSmoothTree(const std::shared_ptr<Node>& node, Number<Inexact> maxRStep);
 };
 
 } // namespace cartocrow::flow_map
diff --git a/cartocrow/flow_map/smooth_tree_painting.cpp b/cartocrow/flow_map/smooth_tree_painting.cpp
index 5d61cd303..21bef58b7 100644
--- a/cartocrow/flow_map/smooth_tree_painting.cpp
+++ b/cartocrow/flow_map/smooth_tree_painting.cpp
@@ -30,9 +30,7 @@ void SmoothTreePainting::paintNodes(renderer::GeometryRenderer& renderer) const
 	renderer.setMode(renderer::GeometryRenderer::vertices);
 	renderer.setStroke(Color{100, 100, 100}, 4);
 	for (const auto& node : m_tree->nodes()) {
-		if (!node->isSteiner()) {
-			renderer.draw(node->m_position.toCartesian());
-		}
+		renderer.draw(node->m_position.toCartesian());
 	}
 }
 
diff --git a/cartocrow/flow_map/spiral_tree.cpp b/cartocrow/flow_map/spiral_tree.cpp
index f81692084..756bb9a66 100644
--- a/cartocrow/flow_map/spiral_tree.cpp
+++ b/cartocrow/flow_map/spiral_tree.cpp
@@ -58,6 +58,10 @@ std::vector<SpiralTree::Obstacle>& SpiralTree::obstacles() {
 	return m_obstacles;
 }
 
+std::shared_ptr<Node> SpiralTree::root() {
+	return m_root;
+}
+
 void SpiralTree::addPlace(const std::string& name, const Point<Inexact>& position,
                           Number<Inexact> flow) {
 	auto newPlace = std::make_shared<Place>(name, position, flow);
diff --git a/cartocrow/flow_map/spiral_tree.h b/cartocrow/flow_map/spiral_tree.h
index 38c36228a..de2ae2269 100644
--- a/cartocrow/flow_map/spiral_tree.h
+++ b/cartocrow/flow_map/spiral_tree.h
@@ -96,6 +96,8 @@ class SpiralTree {
 	const std::vector<std::shared_ptr<Node>>& nodes() const;
 	/// Returns a list of the obstacles in this spiral tree.
 	std::vector<Obstacle>& obstacles();
+	/// Returns the root of this spiral tree.
+	std::shared_ptr<Node> root();
 
 	/// Adds a node to the spiral tree.
 	//void addNode(std::shared_ptr<Node> node);

From b19fc5d5f7a77d4f10a97816dbb0f583e21ab6a0 Mon Sep 17 00:00:00 2001
From: Willem Sonke <willem@wimiso.nl>
Date: Wed, 23 Aug 2023 17:25:18 +0200
Subject: [PATCH 04/27] Add a stub for the optimization routine

---
 cartocrow/flow_map/polar_point.cpp                   |  8 ++++++++
 cartocrow/flow_map/polar_point.h                     | 11 +++++------
 cartocrow/flow_map/smooth_tree.cpp                   |  6 ++++++
 cartocrow/flow_map/smooth_tree.h                     |  3 +++
 cartocrow/flow_map/smooth_tree_painting.cpp          |  4 +++-
 .../flow_map/optimization_demo/optimization_demo.cpp | 12 +++++++++---
 demos/flow_map/optimization_demo/optimization_demo.h |  8 ++++++--
 7 files changed, 40 insertions(+), 12 deletions(-)

diff --git a/cartocrow/flow_map/polar_point.cpp b/cartocrow/flow_map/polar_point.cpp
index bef54ad5f..b7eba3905 100644
--- a/cartocrow/flow_map/polar_point.cpp
+++ b/cartocrow/flow_map/polar_point.cpp
@@ -57,6 +57,14 @@ const Number<Inexact> PolarPoint::phi() const {
 	return m_phi;
 }
 
+void PolarPoint::setR(Number<Inexact> r) {
+	m_r = r;
+}
+
+void PolarPoint::setPhi(Number<Inexact> phi) {
+	m_phi = phi;
+}
+
 Point<Inexact> PolarPoint::toCartesian() const {
 	const Vector<Inexact> d = Vector<Inexact>(std::cos(phi()), std::sin(phi()));
 	return Point<Inexact>(CGAL::ORIGIN) + r() * d;
diff --git a/cartocrow/flow_map/polar_point.h b/cartocrow/flow_map/polar_point.h
index 9aa8437cf..505310d66 100644
--- a/cartocrow/flow_map/polar_point.h
+++ b/cartocrow/flow_map/polar_point.h
@@ -41,23 +41,17 @@ class PolarPoint {
   public:
 	/// Constructs a polar point at the origin.
 	PolarPoint();
-
 	/// Constructs a polar point at the origin.
 	[[deprecated]] PolarPoint(const CGAL::Origin& o);
-
 	/// Constructs a polar point with given \f$r\f$ and \f$\phi\f$.
 	PolarPoint(const Number<Inexact>& R, const Number<Inexact>& phi);
-
 	/// Copy constructor.
 	PolarPoint(const PolarPoint& p);
-
 	/// Constructs a polar point that corresponds to \f$p + t\f$, where \f$p\f$
 	/// is a polar point and \f$t\f$ is a vector in Cartesian coordinates.
 	PolarPoint(const PolarPoint& p, const Vector<Inexact>& t);
-
 	/// Constructs a polar point from a point in Cartesian coordinates.
 	explicit PolarPoint(const Point<Inexact>& p);
-
 	PolarPoint(const Point<Inexact>& p, const Vector<Inexact>& t);
 
 	/// Returns the distance \f$r\f$ from the origin.
@@ -67,6 +61,11 @@ class PolarPoint {
 	/// Returns the angle \f$\phi\f$ relative to the origin.
 	const Number<Inexact> phi() const;
 
+	/// Sets the distance \f$r\f$ from the origin.
+	void setR(Number<Inexact> r);
+	/// Sets the distance \f$\phi\f$ relative to the origin.
+	void setPhi(Number<Inexact> phi);
+
 	/// Returns the point in Cartesian coordinates corresponding to this polar
 	/// point.
 	Point<Inexact> toCartesian() const;
diff --git a/cartocrow/flow_map/smooth_tree.cpp b/cartocrow/flow_map/smooth_tree.cpp
index e33f6a1b0..58f2d8c8d 100644
--- a/cartocrow/flow_map/smooth_tree.cpp
+++ b/cartocrow/flow_map/smooth_tree.cpp
@@ -65,4 +65,10 @@ const std::vector<std::shared_ptr<Node>>& SmoothTree::nodes() const {
 	return m_nodes;
 }
 
+void SmoothTree::optimize() {
+	for (auto& node : m_nodes) {
+		node->m_position.setPhi(node->m_position.phi() + 0.1);
+	}
+}
+
 } // namespace cartocrow::flow_map
diff --git a/cartocrow/flow_map/smooth_tree.h b/cartocrow/flow_map/smooth_tree.h
index 12a87803b..cedab3cae 100644
--- a/cartocrow/flow_map/smooth_tree.h
+++ b/cartocrow/flow_map/smooth_tree.h
@@ -42,6 +42,9 @@ class SmoothTree {
 	/// Returns a list of the nodes in this smooth tree.
 	const std::vector<std::shared_ptr<Node>>& nodes() const;
 
+	/// Performs one optimization step.
+	void optimize();
+
   private:
 	/// The spiral tree underlying this smooth tree.
 	std::shared_ptr<SpiralTree> m_tree;
diff --git a/cartocrow/flow_map/smooth_tree_painting.cpp b/cartocrow/flow_map/smooth_tree_painting.cpp
index 21bef58b7..1a730164a 100644
--- a/cartocrow/flow_map/smooth_tree_painting.cpp
+++ b/cartocrow/flow_map/smooth_tree_painting.cpp
@@ -30,7 +30,9 @@ void SmoothTreePainting::paintNodes(renderer::GeometryRenderer& renderer) const
 	renderer.setMode(renderer::GeometryRenderer::vertices);
 	renderer.setStroke(Color{100, 100, 100}, 4);
 	for (const auto& node : m_tree->nodes()) {
-		renderer.draw(node->m_position.toCartesian());
+		if (node->getType() == Node::ConnectionType::kJoin) {
+			renderer.draw(node->m_position.toCartesian());
+		}
 	}
 }
 
diff --git a/demos/flow_map/optimization_demo/optimization_demo.cpp b/demos/flow_map/optimization_demo/optimization_demo.cpp
index 414020289..d9853283c 100644
--- a/demos/flow_map/optimization_demo/optimization_demo.cpp
+++ b/demos/flow_map/optimization_demo/optimization_demo.cpp
@@ -29,7 +29,6 @@ along with this program.  If not, see <https://www.gnu.org/licenses/>.
 #include "cartocrow/flow_map/parameters.h"
 #include "cartocrow/flow_map/place.h"
 #include "cartocrow/flow_map/reachable_region_algorithm.h"
-#include "cartocrow/flow_map/smooth_tree.h"
 #include "cartocrow/flow_map/smooth_tree_painting.h"
 #include "cartocrow/flow_map/spiral_tree.h"
 #include "cartocrow/flow_map/spiral_tree_obstructed_algorithm.h"
@@ -58,6 +57,13 @@ OptimizationDemo::OptimizationDemo() {
 
 	QToolBar* toolBar = new QToolBar();
 	toolBar->addSeparator();
+	m_optimizeButton = new QPushButton("Optimize");
+	toolBar->addWidget(m_optimizeButton);
+	connect(m_optimizeButton, &QPushButton::clicked, [&]() {
+		m_smoothTree->optimize();
+		m_renderer->update();
+	});
+	toolBar->addSeparator();
 	toolBar->addWidget(new QLabel("α = "));
 	m_alphaSlider = new QSlider(Qt::Horizontal);
 	m_alphaSlider->setMinimum(0);
@@ -91,9 +97,9 @@ void OptimizationDemo::recalculate() {
 	m_renderer->clear();
 	ReachableRegionAlgorithm::ReachableRegion reachableRegion = ReachableRegionAlgorithm(tree).run();
 	SpiralTreeObstructedAlgorithm(tree, reachableRegion).run();
-	auto smoothTree = std::make_shared<SmoothTree>(tree);
+	m_smoothTree = std::make_shared<SmoothTree>(tree);
 	SmoothTreePainting::Options options;
-	auto painting = std::make_shared<SmoothTreePainting>(smoothTree, options);
+	auto painting = std::make_shared<SmoothTreePainting>(m_smoothTree, options);
 	m_renderer->addPainting(painting, "Smooth tree");
 
 	m_renderer->update();
diff --git a/demos/flow_map/optimization_demo/optimization_demo.h b/demos/flow_map/optimization_demo/optimization_demo.h
index 8b331906e..9a60f4433 100644
--- a/demos/flow_map/optimization_demo/optimization_demo.h
+++ b/demos/flow_map/optimization_demo/optimization_demo.h
@@ -21,13 +21,15 @@ along with this program.  If not, see <https://www.gnu.org/licenses/>.
 
 #include <QLabel>
 #include <QMainWindow>
+#include <QPushButton>
 #include <QSlider>
 
 #include <memory>
 #include <vector>
 
-#include <cartocrow/core/core.h>
-#include <cartocrow/renderer/geometry_widget.h>
+#include "cartocrow/core/core.h"
+#include "cartocrow/flow_map/smooth_tree.h"
+#include "cartocrow/renderer/geometry_widget.h"
 
 using namespace cartocrow;
 using namespace cartocrow::renderer;
@@ -43,8 +45,10 @@ class OptimizationDemo : public QMainWindow {
 	Number<Inexact> m_alpha = 25 * M_PI / 180;
 
 	std::vector<std::shared_ptr<Point<Inexact>>> m_places;
+	std::shared_ptr<cartocrow::flow_map::SmoothTree> m_smoothTree;
 
 	GeometryWidget* m_renderer;
 	QSlider* m_alphaSlider;
 	QLabel* m_alphaLabel;
+	QPushButton* m_optimizeButton;
 };

From f6e0409277e58a725534aee93189874aaaa173f8 Mon Sep 17 00:00:00 2001
From: Willem Sonke <willem@wimiso.nl>
Date: Thu, 24 Aug 2023 12:06:44 +0200
Subject: [PATCH 05/27] Run optimization in a QTimer in the demo

---
 cartocrow/flow_map/smooth_tree.cpp                     | 2 +-
 demos/flow_map/optimization_demo/optimization_demo.cpp | 9 +++++++++
 demos/flow_map/optimization_demo/optimization_demo.h   | 2 ++
 3 files changed, 12 insertions(+), 1 deletion(-)

diff --git a/cartocrow/flow_map/smooth_tree.cpp b/cartocrow/flow_map/smooth_tree.cpp
index 58f2d8c8d..b64a0c4ae 100644
--- a/cartocrow/flow_map/smooth_tree.cpp
+++ b/cartocrow/flow_map/smooth_tree.cpp
@@ -67,7 +67,7 @@ const std::vector<std::shared_ptr<Node>>& SmoothTree::nodes() const {
 
 void SmoothTree::optimize() {
 	for (auto& node : m_nodes) {
-		node->m_position.setPhi(node->m_position.phi() + 0.1);
+		node->m_position.setPhi(node->m_position.phi() + 0.001);
 	}
 }
 
diff --git a/demos/flow_map/optimization_demo/optimization_demo.cpp b/demos/flow_map/optimization_demo/optimization_demo.cpp
index d9853283c..8e6873530 100644
--- a/demos/flow_map/optimization_demo/optimization_demo.cpp
+++ b/demos/flow_map/optimization_demo/optimization_demo.cpp
@@ -58,8 +58,17 @@ OptimizationDemo::OptimizationDemo() {
 	QToolBar* toolBar = new QToolBar();
 	toolBar->addSeparator();
 	m_optimizeButton = new QPushButton("Optimize");
+	m_optimizeButton->setCheckable(true);
 	toolBar->addWidget(m_optimizeButton);
+	m_optimizeTimer = new QTimer();
 	connect(m_optimizeButton, &QPushButton::clicked, [&]() {
+		if (m_optimizeButton->isChecked()) {
+			m_optimizeTimer->start(0);
+		} else {
+			m_optimizeTimer->stop();
+		}
+	});
+	connect(m_optimizeTimer, &QTimer::timeout, [&]() {
 		m_smoothTree->optimize();
 		m_renderer->update();
 	});
diff --git a/demos/flow_map/optimization_demo/optimization_demo.h b/demos/flow_map/optimization_demo/optimization_demo.h
index 9a60f4433..448e9fbb6 100644
--- a/demos/flow_map/optimization_demo/optimization_demo.h
+++ b/demos/flow_map/optimization_demo/optimization_demo.h
@@ -23,6 +23,7 @@ along with this program.  If not, see <https://www.gnu.org/licenses/>.
 #include <QMainWindow>
 #include <QPushButton>
 #include <QSlider>
+#include <QTimer>
 
 #include <memory>
 #include <vector>
@@ -51,4 +52,5 @@ class OptimizationDemo : public QMainWindow {
 	QSlider* m_alphaSlider;
 	QLabel* m_alphaLabel;
 	QPushButton* m_optimizeButton;
+	QTimer* m_optimizeTimer;
 };

From eca3d916f8e2cfd017e3bf9d6a4ded92b1876ab8 Mon Sep 17 00:00:00 2001
From: Willem Sonke <willem@wimiso.nl>
Date: Thu, 24 Aug 2023 17:47:49 +0200
Subject: [PATCH 06/27] Implement smoothing cost

---
 cartocrow/flow_map/smooth_tree.cpp            | 42 +++++++++++++++-
 cartocrow/flow_map/smooth_tree.h              |  5 +-
 .../optimization_demo/optimization_demo.cpp   | 48 +++++++++++++++----
 .../optimization_demo/optimization_demo.h     |  1 +
 4 files changed, 83 insertions(+), 13 deletions(-)

diff --git a/cartocrow/flow_map/smooth_tree.cpp b/cartocrow/flow_map/smooth_tree.cpp
index b64a0c4ae..abae32695 100644
--- a/cartocrow/flow_map/smooth_tree.cpp
+++ b/cartocrow/flow_map/smooth_tree.cpp
@@ -21,6 +21,7 @@ along with this program.  If not, see <https://www.gnu.org/licenses/>.
 
 #include "smooth_tree.h"
 #include "cartocrow/flow_map/smooth_tree_painting.h"
+#include <limits>
 
 namespace cartocrow::flow_map {
 
@@ -65,9 +66,46 @@ const std::vector<std::shared_ptr<Node>>& SmoothTree::nodes() const {
 	return m_nodes;
 }
 
+Number<Inexact> SmoothTree::computeSmoothFunction(const std::shared_ptr<Node>& node) {
+	Point<Inexact> n = node->m_position.toCartesian();
+	Point<Inexact> p = node->m_parent->m_position.toCartesian();
+	Point<Inexact> c = node->m_children[0]->m_position.toCartesian();
+	return std::pow(
+	    std::atan2(c.y() - n.y(), c.x() - n.x()) - std::atan2(n.y() - p.y(), n.x() - p.x()), 2);
+}
+
+Number<Inexact> SmoothTree::computeSmoothForce(const std::shared_ptr<Node>& node) {
+	Point<Inexact> n = node->m_position.toCartesian();
+	Point<Inexact> p = node->m_parent->m_position.toCartesian();
+	Point<Inexact> c = node->m_children[0]->m_position.toCartesian();
+	auto datan = [&](Point<Inexact> p1, Point<Inexact> p2) {
+		Number<Inexact> r = std::hypot(p1.x(), p1.y());
+		Number<Inexact> phi = std::atan2(p1.y(), p1.x());
+		return r * (r - p2.x() * std::cos(phi) - p2.y() * std::sin(phi)) /
+		       (std::pow(r, 2) - 2 * r * p2.x() * std::cos(phi) - 2 * r * p2.y() * std::sin(phi) +
+		        std::pow(p2.x(), 2) + std::pow(p2.y(), 2));
+	};
+	return -2 *
+	       (std::atan2(c.y() - n.y(), c.x() - n.x()) - std::atan2(n.y() - p.y(), n.x() - p.x())) *
+	       (datan(n, c) - datan(p, n));
+}
+
 void SmoothTree::optimize() {
-	for (auto& node : m_nodes) {
-		node->m_position.setPhi(node->m_position.phi() + 0.001);
+	std::vector<Number<Inexact>> forces(m_nodes.size(), 0);
+	for (int i = 0; i < m_nodes.size(); i++) {
+		const auto& node = m_nodes[i];
+		if (node->getType() == Node::ConnectionType::kSubdivision) {
+			{
+				std::cout << i << ":   " << computeSmoothFunction(node) << " -- "
+				          << computeSmoothForce(node) << std::endl;
+			}
+			forces[i] += computeSmoothForce(node);
+		}
+	}
+	for (int i = 0; i < m_nodes.size(); i++) {
+		if (forces[i] != std::numeric_limits<Number<Inexact>>::infinity()) {
+			m_nodes[i]->m_position.setPhi(m_nodes[i]->m_position.phi() + 0.001 * forces[i]);
+		}
 	}
 }
 
diff --git a/cartocrow/flow_map/smooth_tree.h b/cartocrow/flow_map/smooth_tree.h
index cedab3cae..6b116b2e3 100644
--- a/cartocrow/flow_map/smooth_tree.h
+++ b/cartocrow/flow_map/smooth_tree.h
@@ -45,7 +45,7 @@ class SmoothTree {
 	/// Performs one optimization step.
 	void optimize();
 
-  private:
+  //private:  // TODO temporary
 	/// The spiral tree underlying this smooth tree.
 	std::shared_ptr<SpiralTree> m_tree;
 
@@ -53,6 +53,9 @@ class SmoothTree {
 	std::vector<std::shared_ptr<Node>> m_nodes;
 
 	std::shared_ptr<Node> constructSmoothTree(const std::shared_ptr<Node>& node, Number<Inexact> maxRStep);
+
+	Number<Inexact> computeSmoothFunction(const std::shared_ptr<Node>& node);
+	Number<Inexact> computeSmoothForce(const std::shared_ptr<Node>& node);
 };
 
 } // namespace cartocrow::flow_map
diff --git a/demos/flow_map/optimization_demo/optimization_demo.cpp b/demos/flow_map/optimization_demo/optimization_demo.cpp
index 8e6873530..1e4de1f1d 100644
--- a/demos/flow_map/optimization_demo/optimization_demo.cpp
+++ b/demos/flow_map/optimization_demo/optimization_demo.cpp
@@ -43,12 +43,14 @@ using namespace cartocrow::renderer;
 OptimizationDemo::OptimizationDemo() {
 	setWindowTitle("CartoCrow – Optimization demo");
 
-	m_places.push_back(std::make_shared<Point<Inexact>>(11.2121212, 17.0707070));
-	m_places.push_back(std::make_shared<Point<Inexact>>(13.9393939, -14.1414141));
-	m_places.push_back(std::make_shared<Point<Inexact>>(-4.5454545, -18.9898989));
-	m_places.push_back(std::make_shared<Point<Inexact>>(16.6666666, 6.1616161));
-	m_places.push_back(std::make_shared<Point<Inexact>>(-9.8989898, 13.9393939));
-	m_places.push_back(std::make_shared<Point<Inexact>>(-16.1616161, -2.6262626));
+	//m_places.push_back(std::make_shared<Point<Inexact>>(11.2121212, 17.0707070));
+	//m_places.push_back(std::make_shared<Point<Inexact>>(13.9393939, -14.1414141));
+	//m_places.push_back(std::make_shared<Point<Inexact>>(-4.5454545, -18.9898989));
+	//m_places.push_back(std::make_shared<Point<Inexact>>(16.6666666, 6.1616161));
+	//m_places.push_back(std::make_shared<Point<Inexact>>(-9.8989898, 13.9393939));
+	//m_places.push_back(std::make_shared<Point<Inexact>>(-16.1616161, -2.6262626));
+	/*m_places.push_back(std::make_shared<Point<Inexact>>(2, 1));
+	m_places.push_back(std::make_shared<Point<Inexact>>(2, -1));*/
 
 	m_renderer = new GeometryWidget();
 	m_renderer->setMaxZoom(10000);
@@ -57,7 +59,7 @@ OptimizationDemo::OptimizationDemo() {
 
 	QToolBar* toolBar = new QToolBar();
 	toolBar->addSeparator();
-	m_optimizeButton = new QPushButton("Optimize");
+	m_optimizeButton = new QPushButton("Run optimization");
 	m_optimizeButton->setCheckable(true);
 	toolBar->addWidget(m_optimizeButton);
 	m_optimizeTimer = new QTimer();
@@ -72,6 +74,12 @@ OptimizationDemo::OptimizationDemo() {
 		m_smoothTree->optimize();
 		m_renderer->update();
 	});
+	m_optimizeOneStepButton = new QPushButton("Optimize one step");
+	toolBar->addWidget(m_optimizeOneStepButton);
+	connect(m_optimizeOneStepButton, &QPushButton::clicked, [&]() {
+		m_smoothTree->optimize();
+		m_renderer->update();
+	});
 	toolBar->addSeparator();
 	toolBar->addWidget(new QLabel("α = "));
 	m_alphaSlider = new QSlider(Qt::Horizontal);
@@ -97,16 +105,36 @@ OptimizationDemo::OptimizationDemo() {
 }
 
 void OptimizationDemo::recalculate() {
+	m_renderer->clear();
 	auto tree = std::make_shared<SpiralTree>(Point<Inexact>(0, 0), m_alpha);
-	for (const auto& place : m_places) {
+	/*for (const auto& place : m_places) {
 		tree->addPlace("", *place, 1);
 	}
 	tree->addShields();
 
-	m_renderer->clear();
 	ReachableRegionAlgorithm::ReachableRegion reachableRegion = ReachableRegionAlgorithm(tree).run();
-	SpiralTreeObstructedAlgorithm(tree, reachableRegion).run();
+	SpiralTreeObstructedAlgorithm(tree, reachableRegion).run();*/
 	m_smoothTree = std::make_shared<SmoothTree>(tree);
+	auto n1 = std::make_shared<Node>(PolarPoint(Point<Inexact>(1, 0.2)));
+	n1->m_parent = m_smoothTree->m_nodes[0];
+	m_smoothTree->m_nodes[0]->m_children.push_back(n1);
+	auto n2 = std::make_shared<Node>(PolarPoint(Point<Inexact>(2, 0)));
+	n2->m_parent = n1;
+	n1->m_children.push_back(n2);
+	auto n3 = std::make_shared<Node>(PolarPoint(Point<Inexact>(3, 0.2)));
+	n3->m_parent = n2;
+	n2->m_children.push_back(n3);
+	auto n4 = std::make_shared<Node>(PolarPoint(Point<Inexact>(4, 0)));
+	n4->m_parent = n3;
+	n3->m_children.push_back(n4);
+	auto n5 = std::make_shared<Node>(PolarPoint(Point<Inexact>(5, 0.2)));
+	n5->m_parent = n4;
+	n4->m_children.push_back(n5);
+	m_smoothTree->m_nodes.push_back(n1);
+	m_smoothTree->m_nodes.push_back(n2);
+	m_smoothTree->m_nodes.push_back(n3);
+	m_smoothTree->m_nodes.push_back(n4);
+	m_smoothTree->m_nodes.push_back(n5);
 	SmoothTreePainting::Options options;
 	auto painting = std::make_shared<SmoothTreePainting>(m_smoothTree, options);
 	m_renderer->addPainting(painting, "Smooth tree");
diff --git a/demos/flow_map/optimization_demo/optimization_demo.h b/demos/flow_map/optimization_demo/optimization_demo.h
index 448e9fbb6..e02b36f3e 100644
--- a/demos/flow_map/optimization_demo/optimization_demo.h
+++ b/demos/flow_map/optimization_demo/optimization_demo.h
@@ -52,5 +52,6 @@ class OptimizationDemo : public QMainWindow {
 	QSlider* m_alphaSlider;
 	QLabel* m_alphaLabel;
 	QPushButton* m_optimizeButton;
+	QPushButton* m_optimizeOneStepButton;
 	QTimer* m_optimizeTimer;
 };

From 3fca5e2e0d18deccb9326fb73f3527fb817685b5 Mon Sep 17 00:00:00 2001
From: Willem Sonke <willem@wimiso.nl>
Date: Thu, 24 Aug 2023 18:13:41 +0200
Subject: [PATCH 07/27] Fix bug in the smooth tree generation

The last subdivision node in an edge would have its child node set
incorrectly.
---
 cartocrow/flow_map/smooth_tree.cpp | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/cartocrow/flow_map/smooth_tree.cpp b/cartocrow/flow_map/smooth_tree.cpp
index abae32695..c6f932c97 100644
--- a/cartocrow/flow_map/smooth_tree.cpp
+++ b/cartocrow/flow_map/smooth_tree.cpp
@@ -56,7 +56,7 @@ std::shared_ptr<Node> SmoothTree::constructSmoothTree(const std::shared_ptr<Node
 			previous = subdivision;
 		}
 		smoothChild->m_parent = previous;
-		previous->m_children.push_back(smoothChild->m_parent);
+		previous->m_children.push_back(smoothChild);
 	}
 
 	return smoothNode;

From 60abbe0537e4756223868ac0637cadf79c4606bf Mon Sep 17 00:00:00 2001
From: Willem Sonke <willem@wimiso.nl>
Date: Thu, 24 Aug 2023 18:16:52 +0200
Subject: [PATCH 08/27] Various fixes to the smoothing cost

---
 cartocrow/flow_map/smooth_tree.cpp            |  6 +---
 cartocrow/flow_map/smooth_tree_painting.cpp   | 16 +++++++--
 .../optimization_demo/optimization_demo.cpp   | 36 +++++--------------
 3 files changed, 22 insertions(+), 36 deletions(-)

diff --git a/cartocrow/flow_map/smooth_tree.cpp b/cartocrow/flow_map/smooth_tree.cpp
index c6f932c97..58a335ca3 100644
--- a/cartocrow/flow_map/smooth_tree.cpp
+++ b/cartocrow/flow_map/smooth_tree.cpp
@@ -78,7 +78,7 @@ Number<Inexact> SmoothTree::computeSmoothForce(const std::shared_ptr<Node>& node
 	Point<Inexact> n = node->m_position.toCartesian();
 	Point<Inexact> p = node->m_parent->m_position.toCartesian();
 	Point<Inexact> c = node->m_children[0]->m_position.toCartesian();
-	auto datan = [&](Point<Inexact> p1, Point<Inexact> p2) {
+	auto datan = [](Point<Inexact> p1, Point<Inexact> p2) {
 		Number<Inexact> r = std::hypot(p1.x(), p1.y());
 		Number<Inexact> phi = std::atan2(p1.y(), p1.x());
 		return r * (r - p2.x() * std::cos(phi) - p2.y() * std::sin(phi)) /
@@ -95,10 +95,6 @@ void SmoothTree::optimize() {
 	for (int i = 0; i < m_nodes.size(); i++) {
 		const auto& node = m_nodes[i];
 		if (node->getType() == Node::ConnectionType::kSubdivision) {
-			{
-				std::cout << i << ":   " << computeSmoothFunction(node) << " -- "
-				          << computeSmoothForce(node) << std::endl;
-			}
 			forces[i] += computeSmoothForce(node);
 		}
 	}
diff --git a/cartocrow/flow_map/smooth_tree_painting.cpp b/cartocrow/flow_map/smooth_tree_painting.cpp
index 1a730164a..da5cc80d0 100644
--- a/cartocrow/flow_map/smooth_tree_painting.cpp
+++ b/cartocrow/flow_map/smooth_tree_painting.cpp
@@ -27,12 +27,22 @@ void SmoothTreePainting::paintFlow(renderer::GeometryRenderer& renderer) const {
 }
 
 void SmoothTreePainting::paintNodes(renderer::GeometryRenderer& renderer) const {
-	renderer.setMode(renderer::GeometryRenderer::vertices);
+	renderer.setMode(renderer::GeometryRenderer::vertices | renderer::GeometryRenderer::stroke);
 	renderer.setStroke(Color{100, 100, 100}, 4);
-	for (const auto& node : m_tree->nodes()) {
+	for (int i = 0; i < m_tree->nodes().size(); i++) {
+		const auto& node = m_tree->nodes()[i];
 		if (node->getType() == Node::ConnectionType::kJoin) {
-			renderer.draw(node->m_position.toCartesian());
+			renderer.setStroke(Color{100, 100, 100}, 4);
+		} else if (node->getType() == Node::ConnectionType::kSubdivision) {
+			renderer.setStroke(Color{255, 84, 32}, 4);
+		} else if (node->getType() == Node::ConnectionType::kLeaf) {
+			renderer.setStroke(Color{84, 160, 32}, 4);
+		} else if (node->getType() == Node::ConnectionType::kRoot) {
+			renderer.setStroke(Color{0, 0, 0}, 4);
 		}
+		renderer.draw(node->m_position.toCartesian());
+		renderer.setStroke(Color{50, 50, 50}, 4);
+		renderer.drawText(node->m_position.toCartesian(), std::to_string(i));
 	}
 }
 
diff --git a/demos/flow_map/optimization_demo/optimization_demo.cpp b/demos/flow_map/optimization_demo/optimization_demo.cpp
index 1e4de1f1d..2e5e158ff 100644
--- a/demos/flow_map/optimization_demo/optimization_demo.cpp
+++ b/demos/flow_map/optimization_demo/optimization_demo.cpp
@@ -43,12 +43,12 @@ using namespace cartocrow::renderer;
 OptimizationDemo::OptimizationDemo() {
 	setWindowTitle("CartoCrow – Optimization demo");
 
-	//m_places.push_back(std::make_shared<Point<Inexact>>(11.2121212, 17.0707070));
-	//m_places.push_back(std::make_shared<Point<Inexact>>(13.9393939, -14.1414141));
-	//m_places.push_back(std::make_shared<Point<Inexact>>(-4.5454545, -18.9898989));
-	//m_places.push_back(std::make_shared<Point<Inexact>>(16.6666666, 6.1616161));
-	//m_places.push_back(std::make_shared<Point<Inexact>>(-9.8989898, 13.9393939));
-	//m_places.push_back(std::make_shared<Point<Inexact>>(-16.1616161, -2.6262626));
+	m_places.push_back(std::make_shared<Point<Inexact>>(11.2121212, 17.0707070));
+	m_places.push_back(std::make_shared<Point<Inexact>>(13.9393939, -14.1414141));
+	m_places.push_back(std::make_shared<Point<Inexact>>(-4.5454545, -18.9898989));
+	m_places.push_back(std::make_shared<Point<Inexact>>(16.6666666, 6.1616161));
+	m_places.push_back(std::make_shared<Point<Inexact>>(-9.8989898, 13.9393939));
+	m_places.push_back(std::make_shared<Point<Inexact>>(-16.1616161, -2.6262626));
 	/*m_places.push_back(std::make_shared<Point<Inexact>>(2, 1));
 	m_places.push_back(std::make_shared<Point<Inexact>>(2, -1));*/
 
@@ -107,34 +107,14 @@ OptimizationDemo::OptimizationDemo() {
 void OptimizationDemo::recalculate() {
 	m_renderer->clear();
 	auto tree = std::make_shared<SpiralTree>(Point<Inexact>(0, 0), m_alpha);
-	/*for (const auto& place : m_places) {
+	for (const auto& place : m_places) {
 		tree->addPlace("", *place, 1);
 	}
 	tree->addShields();
 
 	ReachableRegionAlgorithm::ReachableRegion reachableRegion = ReachableRegionAlgorithm(tree).run();
-	SpiralTreeObstructedAlgorithm(tree, reachableRegion).run();*/
+	SpiralTreeObstructedAlgorithm(tree, reachableRegion).run();
 	m_smoothTree = std::make_shared<SmoothTree>(tree);
-	auto n1 = std::make_shared<Node>(PolarPoint(Point<Inexact>(1, 0.2)));
-	n1->m_parent = m_smoothTree->m_nodes[0];
-	m_smoothTree->m_nodes[0]->m_children.push_back(n1);
-	auto n2 = std::make_shared<Node>(PolarPoint(Point<Inexact>(2, 0)));
-	n2->m_parent = n1;
-	n1->m_children.push_back(n2);
-	auto n3 = std::make_shared<Node>(PolarPoint(Point<Inexact>(3, 0.2)));
-	n3->m_parent = n2;
-	n2->m_children.push_back(n3);
-	auto n4 = std::make_shared<Node>(PolarPoint(Point<Inexact>(4, 0)));
-	n4->m_parent = n3;
-	n3->m_children.push_back(n4);
-	auto n5 = std::make_shared<Node>(PolarPoint(Point<Inexact>(5, 0.2)));
-	n5->m_parent = n4;
-	n4->m_children.push_back(n5);
-	m_smoothTree->m_nodes.push_back(n1);
-	m_smoothTree->m_nodes.push_back(n2);
-	m_smoothTree->m_nodes.push_back(n3);
-	m_smoothTree->m_nodes.push_back(n4);
-	m_smoothTree->m_nodes.push_back(n5);
 	SmoothTreePainting::Options options;
 	auto painting = std::make_shared<SmoothTreePainting>(m_smoothTree, options);
 	m_renderer->addPainting(painting, "Smooth tree");

From 0b818aedbb9a52b91f35bc34ee0391b053dbcc3d Mon Sep 17 00:00:00 2001
From: Willem Sonke <willem@wimiso.nl>
Date: Thu, 24 Aug 2023 19:03:06 +0200
Subject: [PATCH 09/27] Implement smoothing costs properly

I misunderstood the paper first and implemented it based on the
straight-edge angles, but rather it's supposed to work based on the
spiral angle parameters.

This works quite well but some unstabilities occur.
---
 cartocrow/flow_map/smooth_tree.cpp            | 28 +++++++------------
 cartocrow/flow_map/spiral.cpp                 | 27 ++++++++++++++----
 cartocrow/flow_map/spiral.h                   |  8 ++++++
 .../optimization_demo/optimization_demo.cpp   |  5 ++++
 4 files changed, 44 insertions(+), 24 deletions(-)

diff --git a/cartocrow/flow_map/smooth_tree.cpp b/cartocrow/flow_map/smooth_tree.cpp
index 58a335ca3..6ed7149a5 100644
--- a/cartocrow/flow_map/smooth_tree.cpp
+++ b/cartocrow/flow_map/smooth_tree.cpp
@@ -67,27 +67,19 @@ const std::vector<std::shared_ptr<Node>>& SmoothTree::nodes() const {
 }
 
 Number<Inexact> SmoothTree::computeSmoothFunction(const std::shared_ptr<Node>& node) {
-	Point<Inexact> n = node->m_position.toCartesian();
-	Point<Inexact> p = node->m_parent->m_position.toCartesian();
-	Point<Inexact> c = node->m_children[0]->m_position.toCartesian();
-	return std::pow(
-	    std::atan2(c.y() - n.y(), c.x() - n.x()) - std::atan2(n.y() - p.y(), n.x() - p.x()), 2);
+	PolarPoint n = node->m_position;
+	PolarPoint p = node->m_parent->m_position;
+	PolarPoint c = node->m_children[0]->m_position;
+	return std::pow(Spiral::alphaBetweenPoints(p, n) - Spiral::alphaBetweenPoints(n, c), 2);
 }
 
 Number<Inexact> SmoothTree::computeSmoothForce(const std::shared_ptr<Node>& node) {
-	Point<Inexact> n = node->m_position.toCartesian();
-	Point<Inexact> p = node->m_parent->m_position.toCartesian();
-	Point<Inexact> c = node->m_children[0]->m_position.toCartesian();
-	auto datan = [](Point<Inexact> p1, Point<Inexact> p2) {
-		Number<Inexact> r = std::hypot(p1.x(), p1.y());
-		Number<Inexact> phi = std::atan2(p1.y(), p1.x());
-		return r * (r - p2.x() * std::cos(phi) - p2.y() * std::sin(phi)) /
-		       (std::pow(r, 2) - 2 * r * p2.x() * std::cos(phi) - 2 * r * p2.y() * std::sin(phi) +
-		        std::pow(p2.x(), 2) + std::pow(p2.y(), 2));
-	};
-	return -2 *
-	       (std::atan2(c.y() - n.y(), c.x() - n.x()) - std::atan2(n.y() - p.y(), n.x() - p.x())) *
-	       (datan(n, c) - datan(p, n));
+	PolarPoint n = node->m_position;
+	PolarPoint p = node->m_parent->m_position;
+	PolarPoint c = node->m_children[0]->m_position;
+	// chain rule: the derivative of (β_1 - β_2)² is 2(β_1 - β_2) * [β_1 - β_2]'
+	return 2 * (Spiral::alphaBetweenPoints(p, n) - Spiral::alphaBetweenPoints(n, c)) *
+	       (Spiral::alphaBetweenPointsDerivative(p, n) - Spiral::alphaBetweenPointsDerivative(n, c));
 }
 
 void SmoothTree::optimize() {
diff --git a/cartocrow/flow_map/spiral.cpp b/cartocrow/flow_map/spiral.cpp
index 6b6d107f0..e1263fc21 100644
--- a/cartocrow/flow_map/spiral.cpp
+++ b/cartocrow/flow_map/spiral.cpp
@@ -34,17 +34,18 @@ Spiral::Spiral(const PolarPoint& anchor, const Number<Inexact>& angle)
 	}
 }
 
-Spiral::Spiral(const PolarPoint& p1, const PolarPoint& p2) {
+Spiral::Spiral(const PolarPoint& p1, const PolarPoint& p2)
+    : m_anchor(p1.r() < p2.r() ? p2 : p1), m_angle(alphaBetweenPoints(p1, p2)) {}
+
+Number<Inexact> Spiral::alphaBetweenPoints(const PolarPoint& p1, const PolarPoint& p2) {
 	const PolarPoint& source = p1.r() < p2.r() ? p2 : p1;
 	const PolarPoint& target = p1.r() < p2.r() ? p1 : p2;
 
 	if (source.r() == target.r()) {
 		throw std::runtime_error(
-		    "Tried to construct a spiral containing two points equidistant to the root");
+		    "Cannot compute α for a spiral connecting two points equidistant to the root");
 	}
 
-	m_anchor = source;
-
 	// Computing a spiral through two points (p, q), with unknown angle (alpha):
 	//
 	// p = (R_p, phi_p)  ->  R(0) = R_p; phi(0) = phi_p
@@ -59,13 +60,27 @@ Spiral::Spiral(const PolarPoint& p1, const PolarPoint& p2) {
 	//    =>  alpha = tan^-1((phi_q - phi_p) / -ln(R_q / R_p))
 
 	if (target.r() == 0) {
-		m_angle = 0;
+		return 0;
 	} else {
 		Number<Inexact> diff_phi = wrapAngle(target.phi() - source.phi(), -M_PI);
-		m_angle = std::atan(diff_phi / -std::log(target.r() / source.r()));
+		return std::atan(diff_phi / -std::log(target.r() / source.r()));
 	}
 }
 
+Number<Inexact> Spiral::alphaBetweenPointsDerivative(const PolarPoint& p1, const PolarPoint& p2) {
+	const PolarPoint& source = p1.r() < p2.r() ? p2 : p1;
+	const PolarPoint& target = p1.r() < p2.r() ? p1 : p2;
+
+	if (source.r() == target.r()) {
+		throw std::runtime_error(
+		    "Cannot compute α for a spiral connecting two points equidistant to the root");
+	}
+
+	Number<Inexact> phiDiff = wrapAngle(target.phi() - source.phi(), -M_PI);
+	Number<Inexact> rDiffLog = std::log(target.r() / source.r());
+	return rDiffLog / (phiDiff * phiDiff + rDiffLog * rDiffLog);
+}
+
 const PolarPoint& Spiral::anchor() const {
 	return m_anchor;
 }
diff --git a/cartocrow/flow_map/spiral.h b/cartocrow/flow_map/spiral.h
index af8b5f497..7173d6557 100644
--- a/cartocrow/flow_map/spiral.h
+++ b/cartocrow/flow_map/spiral.h
@@ -68,6 +68,14 @@ class Spiral {
 	/// segment.
 	Spiral(const PolarPoint& p1, const PolarPoint& p2);
 
+	/// Computes the \f$\alpha\f$ of the shortest logarithmic spiral connecting
+	/// the two given points.
+	static Number<Inexact> alphaBetweenPoints(const PolarPoint& p1, const PolarPoint& p2);
+	/// Computes \f$\frac{d\alpha}{d\phi}\f$, where \f$\alpha\f$ is the angle of
+	/// the shortest logarithmic spiral connecting the two given points (see
+	/// \ref alphaBetweenPoints) and \f$\phi\f$ is the polar angle of `p1`.
+	static Number<Inexact> alphaBetweenPointsDerivative(const PolarPoint& p1, const PolarPoint& p2);
+
 	/// Returns the anchor of this spiral.
 	const PolarPoint& anchor() const;
 	/// Returns the angle of this spiral.
diff --git a/demos/flow_map/optimization_demo/optimization_demo.cpp b/demos/flow_map/optimization_demo/optimization_demo.cpp
index 2e5e158ff..a625a07e9 100644
--- a/demos/flow_map/optimization_demo/optimization_demo.cpp
+++ b/demos/flow_map/optimization_demo/optimization_demo.cpp
@@ -23,6 +23,7 @@ along with this program.  If not, see <https://www.gnu.org/licenses/>.
 
 #include <QApplication>
 #include <QCheckBox>
+#include <cstdlib>
 
 #include "cartocrow/core/core.h"
 #include "cartocrow/flow_map/painting.h"
@@ -115,6 +116,10 @@ void OptimizationDemo::recalculate() {
 	ReachableRegionAlgorithm::ReachableRegion reachableRegion = ReachableRegionAlgorithm(tree).run();
 	SpiralTreeObstructedAlgorithm(tree, reachableRegion).run();
 	m_smoothTree = std::make_shared<SmoothTree>(tree);
+	for (auto& node : m_smoothTree->m_nodes) {
+		node->m_position.setPhi(wrapAngle(node->m_position.phi() +
+		                        static_cast<float>(std::rand()) / RAND_MAX * 0.05, -M_PI));
+	}
 	SmoothTreePainting::Options options;
 	auto painting = std::make_shared<SmoothTreePainting>(m_smoothTree, options);
 	m_renderer->addPainting(painting, "Smooth tree");

From 656c635d05f6e54366551d8a1d278ce59878f001 Mon Sep 17 00:00:00 2001
From: Willem Sonke <willem@wimiso.nl>
Date: Fri, 25 Aug 2023 12:17:03 +0200
Subject: [PATCH 10/27] Fix smoothing cost sign issue causing instability

---
 cartocrow/flow_map/smooth_tree.cpp                     | 8 ++++----
 cartocrow/flow_map/spiral.cpp                          | 7 ++++++-
 demos/flow_map/optimization_demo/optimization_demo.cpp | 7 +++++--
 3 files changed, 15 insertions(+), 7 deletions(-)

diff --git a/cartocrow/flow_map/smooth_tree.cpp b/cartocrow/flow_map/smooth_tree.cpp
index 6ed7149a5..bef1643d0 100644
--- a/cartocrow/flow_map/smooth_tree.cpp
+++ b/cartocrow/flow_map/smooth_tree.cpp
@@ -70,7 +70,7 @@ Number<Inexact> SmoothTree::computeSmoothFunction(const std::shared_ptr<Node>& n
 	PolarPoint n = node->m_position;
 	PolarPoint p = node->m_parent->m_position;
 	PolarPoint c = node->m_children[0]->m_position;
-	return std::pow(Spiral::alphaBetweenPoints(p, n) - Spiral::alphaBetweenPoints(n, c), 2);
+	return std::pow(Spiral::alphaBetweenPoints(n, p) - Spiral::alphaBetweenPoints(n, c), 2);
 }
 
 Number<Inexact> SmoothTree::computeSmoothForce(const std::shared_ptr<Node>& node) {
@@ -78,8 +78,8 @@ Number<Inexact> SmoothTree::computeSmoothForce(const std::shared_ptr<Node>& node
 	PolarPoint p = node->m_parent->m_position;
 	PolarPoint c = node->m_children[0]->m_position;
 	// chain rule: the derivative of (β_1 - β_2)² is 2(β_1 - β_2) * [β_1 - β_2]'
-	return 2 * (Spiral::alphaBetweenPoints(p, n) - Spiral::alphaBetweenPoints(n, c)) *
-	       (Spiral::alphaBetweenPointsDerivative(p, n) - Spiral::alphaBetweenPointsDerivative(n, c));
+	return 2 * (Spiral::alphaBetweenPoints(n, p) - Spiral::alphaBetweenPoints(n, c)) *
+	       (Spiral::alphaBetweenPointsDerivative(n, p) - Spiral::alphaBetweenPointsDerivative(n, c));
 }
 
 void SmoothTree::optimize() {
@@ -92,7 +92,7 @@ void SmoothTree::optimize() {
 	}
 	for (int i = 0; i < m_nodes.size(); i++) {
 		if (forces[i] != std::numeric_limits<Number<Inexact>>::infinity()) {
-			m_nodes[i]->m_position.setPhi(m_nodes[i]->m_position.phi() + 0.001 * forces[i]);
+			m_nodes[i]->m_position.setPhi(m_nodes[i]->m_position.phi() - 0.001 * forces[i]);
 		}
 	}
 }
diff --git a/cartocrow/flow_map/spiral.cpp b/cartocrow/flow_map/spiral.cpp
index e1263fc21..b17d4b0ee 100644
--- a/cartocrow/flow_map/spiral.cpp
+++ b/cartocrow/flow_map/spiral.cpp
@@ -78,7 +78,12 @@ Number<Inexact> Spiral::alphaBetweenPointsDerivative(const PolarPoint& p1, const
 
 	Number<Inexact> phiDiff = wrapAngle(target.phi() - source.phi(), -M_PI);
 	Number<Inexact> rDiffLog = std::log(target.r() / source.r());
-	return rDiffLog / (phiDiff * phiDiff + rDiffLog * rDiffLog);
+	Number<Inexact> derivative = -rDiffLog / (phiDiff * phiDiff + rDiffLog * rDiffLog);
+
+	// swap the sign if p1 was further from the origin than p2
+	int sign = p1.r() < p2.r() ? 1 : -1;
+
+	return sign * derivative;
 }
 
 const PolarPoint& Spiral::anchor() const {
diff --git a/demos/flow_map/optimization_demo/optimization_demo.cpp b/demos/flow_map/optimization_demo/optimization_demo.cpp
index a625a07e9..5b8df9249 100644
--- a/demos/flow_map/optimization_demo/optimization_demo.cpp
+++ b/demos/flow_map/optimization_demo/optimization_demo.cpp
@@ -117,8 +117,11 @@ void OptimizationDemo::recalculate() {
 	SpiralTreeObstructedAlgorithm(tree, reachableRegion).run();
 	m_smoothTree = std::make_shared<SmoothTree>(tree);
 	for (auto& node : m_smoothTree->m_nodes) {
-		node->m_position.setPhi(wrapAngle(node->m_position.phi() +
-		                        static_cast<float>(std::rand()) / RAND_MAX * 0.05, -M_PI));
+		if (node->getType() == Node::ConnectionType::kSubdivision) {
+			node->m_position.setPhi(wrapAngle(
+			    node->m_position.phi() + static_cast<float>(std::rand()) / RAND_MAX * 0.5 - 0.25,
+			    -M_PI));
+		}
 	}
 	SmoothTreePainting::Options options;
 	auto painting = std::make_shared<SmoothTreePainting>(m_smoothTree, options);

From f0cc83a1cc87823fcfcf6f89ef003167b5a0f8c7 Mon Sep 17 00:00:00 2001
From: Willem Sonke <willem@wimiso.nl>
Date: Fri, 25 Aug 2023 15:03:08 +0200
Subject: [PATCH 11/27] Use random node locations in the optimization demo

---
 .../optimization_demo/optimization_demo.cpp    | 18 ++++++++++++------
 1 file changed, 12 insertions(+), 6 deletions(-)

diff --git a/demos/flow_map/optimization_demo/optimization_demo.cpp b/demos/flow_map/optimization_demo/optimization_demo.cpp
index 5b8df9249..fd5ab2492 100644
--- a/demos/flow_map/optimization_demo/optimization_demo.cpp
+++ b/demos/flow_map/optimization_demo/optimization_demo.cpp
@@ -44,15 +44,21 @@ using namespace cartocrow::renderer;
 OptimizationDemo::OptimizationDemo() {
 	setWindowTitle("CartoCrow – Optimization demo");
 
-	m_places.push_back(std::make_shared<Point<Inexact>>(11.2121212, 17.0707070));
+	/*m_places.push_back(std::make_shared<Point<Inexact>>(11.2121212, 17.0707070));
 	m_places.push_back(std::make_shared<Point<Inexact>>(13.9393939, -14.1414141));
 	m_places.push_back(std::make_shared<Point<Inexact>>(-4.5454545, -18.9898989));
 	m_places.push_back(std::make_shared<Point<Inexact>>(16.6666666, 6.1616161));
 	m_places.push_back(std::make_shared<Point<Inexact>>(-9.8989898, 13.9393939));
-	m_places.push_back(std::make_shared<Point<Inexact>>(-16.1616161, -2.6262626));
+	m_places.push_back(std::make_shared<Point<Inexact>>(-16.1616161, -2.6262626));*/
 	/*m_places.push_back(std::make_shared<Point<Inexact>>(2, 1));
 	m_places.push_back(std::make_shared<Point<Inexact>>(2, -1));*/
 
+	for (int i = 0; i < 40; i++) {
+		m_places.push_back(std::make_shared<Point<Inexact>>(
+		    static_cast<float>(std::rand()) / RAND_MAX * 200 - 100,
+		    static_cast<float>(std::rand()) / RAND_MAX * 200 - 100));
+	}
+
 	m_renderer = new GeometryWidget();
 	m_renderer->setMaxZoom(10000);
 	m_renderer->setGridMode(GeometryWidget::GridMode::POLAR);
@@ -111,18 +117,18 @@ void OptimizationDemo::recalculate() {
 	for (const auto& place : m_places) {
 		tree->addPlace("", *place, 1);
 	}
-	tree->addShields();
+	//tree->addShields();
 
 	ReachableRegionAlgorithm::ReachableRegion reachableRegion = ReachableRegionAlgorithm(tree).run();
 	SpiralTreeObstructedAlgorithm(tree, reachableRegion).run();
 	m_smoothTree = std::make_shared<SmoothTree>(tree);
-	for (auto& node : m_smoothTree->m_nodes) {
+	/*for (auto& node : m_smoothTree->m_nodes) {
 		if (node->getType() == Node::ConnectionType::kSubdivision) {
 			node->m_position.setPhi(wrapAngle(
-			    node->m_position.phi() + static_cast<float>(std::rand()) / RAND_MAX * 0.5 - 0.25,
+			    node->m_position.phi() + static_cast<float>(std::rand()) / RAND_MAX * 0.2 - 0.1,
 			    -M_PI));
 		}
-	}
+	}*/
 	SmoothTreePainting::Options options;
 	auto painting = std::make_shared<SmoothTreePainting>(m_smoothTree, options);
 	m_renderer->addPainting(painting, "Smooth tree");

From 5e2a11d5ae628ad83287f93604cf0bc116e0cc60 Mon Sep 17 00:00:00 2001
From: Willem Sonke <willem@wimiso.nl>
Date: Thu, 31 Aug 2023 13:48:50 +0200
Subject: [PATCH 12/27] Store each node's flow in the smooth tree

---
 cartocrow/flow_map/node.h                   |  3 ++
 cartocrow/flow_map/smooth_tree.cpp          | 45 ++++++++++++---------
 cartocrow/flow_map/smooth_tree_painting.cpp | 13 ++----
 cartocrow/renderer/geometry_widget.cpp      |  6 +--
 4 files changed, 36 insertions(+), 31 deletions(-)

diff --git a/cartocrow/flow_map/node.h b/cartocrow/flow_map/node.h
index a5c407f4e..51f762bf4 100644
--- a/cartocrow/flow_map/node.h
+++ b/cartocrow/flow_map/node.h
@@ -57,6 +57,9 @@ struct Node {
 	/// leaf nodes, a node with the leaf type may have children if it is located
 	/// inside the spiral region of another node.
 	std::vector<std::shared_ptr<Node>> m_children;
+
+	/// The amount of flow through this node.
+	Number<Inexact> m_flow;
 };
 
 } // namespace cartocrow::flow_map
diff --git a/cartocrow/flow_map/smooth_tree.cpp b/cartocrow/flow_map/smooth_tree.cpp
index bef1643d0..20250a75d 100644
--- a/cartocrow/flow_map/smooth_tree.cpp
+++ b/cartocrow/flow_map/smooth_tree.cpp
@@ -38,25 +38,32 @@ SmoothTree::SmoothTree(const std::shared_ptr<SpiralTree> spiralTree) : m_tree(sp
 std::shared_ptr<Node> SmoothTree::constructSmoothTree(const std::shared_ptr<Node>& node, Number<Inexact> maxRStep) {
 	auto smoothNode = std::make_shared<Node>(node->m_position);
 	m_nodes.push_back(smoothNode);
-	for (const auto& child : node->m_children) {
-		auto smoothChild = constructSmoothTree(child, maxRStep);
-		const Spiral spiral(node->m_position, child->m_position);
-		Number<Inexact> rMin = node->m_position.r();
-		Number<Inexact> rMax = child->m_position.r();
-		Number<Inexact> rRange = rMax - rMin;
-		int subdivisionCount = std::ceil(rRange / maxRStep);
-		auto previous = smoothNode;
-		for (double i = 1; i < subdivisionCount; ++i) {
-			const PolarPoint position =
-			    spiral.evaluate(spiral.parameterForR(rMin + i * (rMax - rMin) / subdivisionCount));
-			auto subdivision = std::make_shared<Node>(position);
-			m_nodes.push_back(subdivision);
-			subdivision->m_parent = previous;
-			previous->m_children.push_back(subdivision);
-			previous = subdivision;
+	if (node->getType() == Node::ConnectionType::kLeaf) {
+		smoothNode->m_flow = node->m_place->m_flow;
+	} else {
+		smoothNode->m_flow = 0;
+		for (const auto& child : node->m_children) {
+			auto smoothChild = constructSmoothTree(child, maxRStep);
+			smoothNode->m_flow += smoothChild->m_flow;
+			const Spiral spiral(node->m_position, child->m_position);
+			Number<Inexact> rMin = node->m_position.r();
+			Number<Inexact> rMax = child->m_position.r();
+			Number<Inexact> rRange = rMax - rMin;
+			int subdivisionCount = std::ceil(rRange / maxRStep);
+			auto previous = smoothNode;
+			for (double i = 1; i < subdivisionCount; ++i) {
+				const PolarPoint position =
+					spiral.evaluate(spiral.parameterForR(rMin + i * (rMax - rMin) / subdivisionCount));
+				auto subdivision = std::make_shared<Node>(position);
+				subdivision->m_flow = smoothChild->m_flow;
+				m_nodes.push_back(subdivision);
+				subdivision->m_parent = previous;
+				previous->m_children.push_back(subdivision);
+				previous = subdivision;
+			}
+			smoothChild->m_parent = previous;
+			previous->m_children.push_back(smoothChild);
 		}
-		smoothChild->m_parent = previous;
-		previous->m_children.push_back(smoothChild);
 	}
 
 	return smoothNode;
@@ -87,7 +94,7 @@ void SmoothTree::optimize() {
 	for (int i = 0; i < m_nodes.size(); i++) {
 		const auto& node = m_nodes[i];
 		if (node->getType() == Node::ConnectionType::kSubdivision) {
-			forces[i] += computeSmoothForce(node);
+			forces[i] += 0.4 * computeSmoothForce(node); // TODO c_S = 0.4
 		}
 	}
 	for (int i = 0; i < m_nodes.size(); i++) {
diff --git a/cartocrow/flow_map/smooth_tree_painting.cpp b/cartocrow/flow_map/smooth_tree_painting.cpp
index da5cc80d0..fd6db0928 100644
--- a/cartocrow/flow_map/smooth_tree_painting.cpp
+++ b/cartocrow/flow_map/smooth_tree_painting.cpp
@@ -16,11 +16,11 @@ void SmoothTreePainting::paint(renderer::GeometryRenderer& renderer) const {
 
 void SmoothTreePainting::paintFlow(renderer::GeometryRenderer& renderer) const {
 	renderer.setMode(renderer::GeometryRenderer::stroke);
-	renderer.setStroke(Color{255, 84, 32}, 4);
 	for (const auto& node : m_tree->nodes()) {
 		if (node->m_parent == nullptr) {
 			continue;
 		}
+		renderer.setStroke(Color{255, 84, 32}, node->m_flow);
 		renderer.draw(Segment<Inexact>(node->m_parent->m_position.toCartesian(),
 		                               node->m_position.toCartesian()));
 	}
@@ -31,18 +31,13 @@ void SmoothTreePainting::paintNodes(renderer::GeometryRenderer& renderer) const
 	renderer.setStroke(Color{100, 100, 100}, 4);
 	for (int i = 0; i < m_tree->nodes().size(); i++) {
 		const auto& node = m_tree->nodes()[i];
-		if (node->getType() == Node::ConnectionType::kJoin) {
-			renderer.setStroke(Color{100, 100, 100}, 4);
-		} else if (node->getType() == Node::ConnectionType::kSubdivision) {
-			renderer.setStroke(Color{255, 84, 32}, 4);
-		} else if (node->getType() == Node::ConnectionType::kLeaf) {
+		if (node->getType() == Node::ConnectionType::kLeaf) {
 			renderer.setStroke(Color{84, 160, 32}, 4);
+			renderer.draw(node->m_position.toCartesian());
 		} else if (node->getType() == Node::ConnectionType::kRoot) {
 			renderer.setStroke(Color{0, 0, 0}, 4);
+			renderer.draw(node->m_position.toCartesian());
 		}
-		renderer.draw(node->m_position.toCartesian());
-		renderer.setStroke(Color{50, 50, 50}, 4);
-		renderer.drawText(node->m_position.toCartesian(), std::to_string(i));
 	}
 }
 
diff --git a/cartocrow/renderer/geometry_widget.cpp b/cartocrow/renderer/geometry_widget.cpp
index 58d956813..d482f255f 100644
--- a/cartocrow/renderer/geometry_widget.cpp
+++ b/cartocrow/renderer/geometry_widget.cpp
@@ -31,6 +31,7 @@ along with this program.  If not, see <https://www.gnu.org/licenses/>.
 
 #include <cmath>
 #include <limits>
+#include <qnamespace.h>
 
 namespace cartocrow::renderer {
 
@@ -475,8 +476,7 @@ void GeometryWidget::draw(const Point<Inexact>& p) {
 }
 
 void GeometryWidget::draw(const Segment<Inexact>& s) {
-	m_painter->setPen(QPen(m_style.m_strokeColor, m_style.m_strokeWidth));
-	m_painter->setBrush(Qt::NoBrush);
+	setupPainter();
 	QPointF p1 = convertPoint(s.start());
 	QPointF p2 = convertPoint(s.end());
 	m_painter->drawLine(p1, p2);
@@ -546,7 +546,7 @@ void GeometryWidget::setupPainter() {
 		m_painter->setBrush(Qt::NoBrush);
 	}
 	if (m_style.m_mode & GeometryRenderer::stroke) {
-		m_painter->setPen(QPen(m_style.m_strokeColor, m_style.m_strokeWidth));
+		m_painter->setPen(QPen(m_style.m_strokeColor, m_style.m_strokeWidth, Qt::SolidLine, Qt::FlatCap, Qt::RoundJoin));
 	} else {
 		m_painter->setPen(Qt::NoPen);
 	}

From 009424f9adeeb66e5582dca0f0380583304f3990 Mon Sep 17 00:00:00 2001
From: Willem Sonke <willem@wimiso.nl>
Date: Fri, 1 Sep 2023 12:00:51 +0200
Subject: [PATCH 13/27] Draw an outline for the flow map

---
 cartocrow/flow_map/smooth_tree.cpp          |  6 +++---
 cartocrow/flow_map/smooth_tree.h            |  4 ++--
 cartocrow/flow_map/smooth_tree_painting.cpp | 12 ++++++++++++
 3 files changed, 17 insertions(+), 5 deletions(-)

diff --git a/cartocrow/flow_map/smooth_tree.cpp b/cartocrow/flow_map/smooth_tree.cpp
index 20250a75d..1fd90a00a 100644
--- a/cartocrow/flow_map/smooth_tree.cpp
+++ b/cartocrow/flow_map/smooth_tree.cpp
@@ -73,14 +73,14 @@ const std::vector<std::shared_ptr<Node>>& SmoothTree::nodes() const {
 	return m_nodes;
 }
 
-Number<Inexact> SmoothTree::computeSmoothFunction(const std::shared_ptr<Node>& node) {
+Number<Inexact> SmoothTree::computeSmoothingFunction(const std::shared_ptr<Node>& node) {
 	PolarPoint n = node->m_position;
 	PolarPoint p = node->m_parent->m_position;
 	PolarPoint c = node->m_children[0]->m_position;
 	return std::pow(Spiral::alphaBetweenPoints(n, p) - Spiral::alphaBetweenPoints(n, c), 2);
 }
 
-Number<Inexact> SmoothTree::computeSmoothForce(const std::shared_ptr<Node>& node) {
+Number<Inexact> SmoothTree::computeSmoothingForce(const std::shared_ptr<Node>& node) {
 	PolarPoint n = node->m_position;
 	PolarPoint p = node->m_parent->m_position;
 	PolarPoint c = node->m_children[0]->m_position;
@@ -94,7 +94,7 @@ void SmoothTree::optimize() {
 	for (int i = 0; i < m_nodes.size(); i++) {
 		const auto& node = m_nodes[i];
 		if (node->getType() == Node::ConnectionType::kSubdivision) {
-			forces[i] += 0.4 * computeSmoothForce(node); // TODO c_S = 0.4
+			forces[i] += 0.4 * computeSmoothingForce(node); // TODO c_S = 0.4
 		}
 	}
 	for (int i = 0; i < m_nodes.size(); i++) {
diff --git a/cartocrow/flow_map/smooth_tree.h b/cartocrow/flow_map/smooth_tree.h
index 6b116b2e3..8d27b845e 100644
--- a/cartocrow/flow_map/smooth_tree.h
+++ b/cartocrow/flow_map/smooth_tree.h
@@ -54,8 +54,8 @@ class SmoothTree {
 
 	std::shared_ptr<Node> constructSmoothTree(const std::shared_ptr<Node>& node, Number<Inexact> maxRStep);
 
-	Number<Inexact> computeSmoothFunction(const std::shared_ptr<Node>& node);
-	Number<Inexact> computeSmoothForce(const std::shared_ptr<Node>& node);
+	Number<Inexact> computeSmoothingFunction(const std::shared_ptr<Node>& node);
+	Number<Inexact> computeSmoothingForce(const std::shared_ptr<Node>& node);
 };
 
 } // namespace cartocrow::flow_map
diff --git a/cartocrow/flow_map/smooth_tree_painting.cpp b/cartocrow/flow_map/smooth_tree_painting.cpp
index fd6db0928..e1fd06f4e 100644
--- a/cartocrow/flow_map/smooth_tree_painting.cpp
+++ b/cartocrow/flow_map/smooth_tree_painting.cpp
@@ -16,6 +16,18 @@ void SmoothTreePainting::paint(renderer::GeometryRenderer& renderer) const {
 
 void SmoothTreePainting::paintFlow(renderer::GeometryRenderer& renderer) const {
 	renderer.setMode(renderer::GeometryRenderer::stroke);
+
+	// boundary
+	for (const auto& node : m_tree->nodes()) {
+		if (node->m_parent == nullptr) {
+			continue;
+		}
+		renderer.setStroke(Color{0, 0, 0}, node->m_flow + 2.0f);
+		renderer.draw(Segment<Inexact>(node->m_parent->m_position.toCartesian(),
+		                               node->m_position.toCartesian()));
+	}
+
+	// interior
 	for (const auto& node : m_tree->nodes()) {
 		if (node->m_parent == nullptr) {
 			continue;

From ad68d3428a5484ab91e21cce8e6d71de0ef59164 Mon Sep 17 00:00:00 2001
From: Willem Sonke <willem@wimiso.nl>
Date: Tue, 5 Sep 2023 19:12:22 +0200
Subject: [PATCH 14/27] Make setup for gradient computations more general

I redid the derivative computations and changed the derivative functions
in Spiral correspondingly.

Also fix the smoothing cost gradient to also apply to the parent and
child, which was missing before.
---
 cartocrow/flow_map/node.h          |  3 ++
 cartocrow/flow_map/smooth_tree.cpp | 43 +++++++++++++++++------
 cartocrow/flow_map/smooth_tree.h   |  8 ++++-
 cartocrow/flow_map/spiral.cpp      | 34 ++++++++++++------
 cartocrow/flow_map/spiral.h        | 55 ++++++++++++++++++++++++++----
 5 files changed, 114 insertions(+), 29 deletions(-)

diff --git a/cartocrow/flow_map/node.h b/cartocrow/flow_map/node.h
index 51f762bf4..9381d6a2b 100644
--- a/cartocrow/flow_map/node.h
+++ b/cartocrow/flow_map/node.h
@@ -60,6 +60,9 @@ struct Node {
 
 	/// The amount of flow through this node.
 	Number<Inexact> m_flow;
+
+	/// The index of this node in its tree.
+	int m_id = -1;
 };
 
 } // namespace cartocrow::flow_map
diff --git a/cartocrow/flow_map/smooth_tree.cpp b/cartocrow/flow_map/smooth_tree.cpp
index 1fd90a00a..285632a5d 100644
--- a/cartocrow/flow_map/smooth_tree.cpp
+++ b/cartocrow/flow_map/smooth_tree.cpp
@@ -77,29 +77,50 @@ Number<Inexact> SmoothTree::computeSmoothingFunction(const std::shared_ptr<Node>
 	PolarPoint n = node->m_position;
 	PolarPoint p = node->m_parent->m_position;
 	PolarPoint c = node->m_children[0]->m_position;
-	return std::pow(Spiral::alphaBetweenPoints(n, p) - Spiral::alphaBetweenPoints(n, c), 2);
+	return std::pow(Spiral::alpha(p, n) - Spiral::alpha(n, c), 2);
 }
 
-Number<Inexact> SmoothTree::computeSmoothingForce(const std::shared_ptr<Node>& node) {
-	PolarPoint n = node->m_position;
-	PolarPoint p = node->m_parent->m_position;
-	PolarPoint c = node->m_children[0]->m_position;
+void SmoothTree::applySmoothingForce(int i, int iParent, int iChild) {
+	PolarPoint n = m_nodes[i]->m_position;
+	PolarPoint p = m_nodes[iParent]->m_position;
+	PolarPoint c = m_nodes[iChild]->m_position;
+
 	// chain rule: the derivative of (β_1 - β_2)² is 2(β_1 - β_2) * [β_1 - β_2]'
-	return 2 * (Spiral::alphaBetweenPoints(n, p) - Spiral::alphaBetweenPoints(n, c)) *
-	       (Spiral::alphaBetweenPointsDerivative(n, p) - Spiral::alphaBetweenPointsDerivative(n, c));
+	m_forces[i].r += 2 * (Spiral::alpha(p, n) - Spiral::alpha(n, c)) *
+	                 (Spiral::dAlphaDR2(p, n) - Spiral::dAlphaDR1(n, c));
+	m_forces[i].phi += 2 * (Spiral::alpha(p, n) - Spiral::alpha(n, c)) *
+	                   (Spiral::dAlphaDPhi2(p, n) - Spiral::dAlphaDPhi1(n, c));
+
+	m_forces[iParent].r += 2 * (Spiral::alpha(p, n) - Spiral::alpha(n, c)) *
+	                 (Spiral::dAlphaDR1(p, n));
+	m_forces[iParent].phi += 2 * (Spiral::alpha(p, n) - Spiral::alpha(n, c)) *
+	                   (Spiral::dAlphaDPhi1(p, n));
+
+	m_forces[iChild].r += 2 * (Spiral::alpha(p, n) - Spiral::alpha(n, c)) *
+	                 (Spiral::dAlphaDR2(n, c));
+	m_forces[iChild].phi += 2 * (Spiral::alpha(p, n) - Spiral::alpha(n, c)) *
+	                   (Spiral::dAlphaDPhi2(n, c));
 }
 
 void SmoothTree::optimize() {
-	std::vector<Number<Inexact>> forces(m_nodes.size(), 0);
+	m_forces = std::vector<PolarForce>(m_nodes.size(), PolarForce{});
+	for (int i = 0; i < m_nodes.size(); i++) {
+		m_nodes[i]->m_id = i;
+	}
 	for (int i = 0; i < m_nodes.size(); i++) {
 		const auto& node = m_nodes[i];
 		if (node->getType() == Node::ConnectionType::kSubdivision) {
-			forces[i] += 0.4 * computeSmoothingForce(node); // TODO c_S = 0.4
+			applySmoothingForce(i, m_nodes[i]->m_parent->m_id, m_nodes[i]->m_children[0]->m_id);
 		}
 	}
 	for (int i = 0; i < m_nodes.size(); i++) {
-		if (forces[i] != std::numeric_limits<Number<Inexact>>::infinity()) {
-			m_nodes[i]->m_position.setPhi(m_nodes[i]->m_position.phi() - 0.001 * forces[i]);
+		Number<Inexact> epsilon = 0.001;
+		if (m_nodes[i]->getType() == Node::ConnectionType::kJoin) {
+			m_nodes[i]->m_position.setR(m_nodes[i]->m_position.r() + epsilon * m_forces[i].r);
+		}
+		if (m_nodes[i]->getType() == Node::ConnectionType::kJoin ||
+		    m_nodes[i]->getType() == Node::ConnectionType::kSubdivision) {
+			m_nodes[i]->m_position.setPhi(m_nodes[i]->m_position.phi() + epsilon * m_forces[i].phi);
 		}
 	}
 }
diff --git a/cartocrow/flow_map/smooth_tree.h b/cartocrow/flow_map/smooth_tree.h
index 8d27b845e..8ca316372 100644
--- a/cartocrow/flow_map/smooth_tree.h
+++ b/cartocrow/flow_map/smooth_tree.h
@@ -54,8 +54,14 @@ class SmoothTree {
 
 	std::shared_ptr<Node> constructSmoothTree(const std::shared_ptr<Node>& node, Number<Inexact> maxRStep);
 
+	struct PolarForce {
+		Number<Inexact> r = 0;
+		Number<Inexact> phi = 0;
+	};
+	std::vector<PolarForce> m_forces;
+
 	Number<Inexact> computeSmoothingFunction(const std::shared_ptr<Node>& node);
-	Number<Inexact> computeSmoothingForce(const std::shared_ptr<Node>& node);
+	void applySmoothingForce(int i, int iParent, int iChild);
 };
 
 } // namespace cartocrow::flow_map
diff --git a/cartocrow/flow_map/spiral.cpp b/cartocrow/flow_map/spiral.cpp
index b17d4b0ee..a9c4a6ef9 100644
--- a/cartocrow/flow_map/spiral.cpp
+++ b/cartocrow/flow_map/spiral.cpp
@@ -35,9 +35,9 @@ Spiral::Spiral(const PolarPoint& anchor, const Number<Inexact>& angle)
 }
 
 Spiral::Spiral(const PolarPoint& p1, const PolarPoint& p2)
-    : m_anchor(p1.r() < p2.r() ? p2 : p1), m_angle(alphaBetweenPoints(p1, p2)) {}
+    : m_anchor(p1.r() < p2.r() ? p2 : p1), m_angle(alpha(p1, p2)) {}
 
-Number<Inexact> Spiral::alphaBetweenPoints(const PolarPoint& p1, const PolarPoint& p2) {
+Number<Inexact> Spiral::alpha(const PolarPoint& p1, const PolarPoint& p2) {
 	const PolarPoint& source = p1.r() < p2.r() ? p2 : p1;
 	const PolarPoint& target = p1.r() < p2.r() ? p1 : p2;
 
@@ -67,23 +67,35 @@ Number<Inexact> Spiral::alphaBetweenPoints(const PolarPoint& p1, const PolarPoin
 	}
 }
 
-Number<Inexact> Spiral::alphaBetweenPointsDerivative(const PolarPoint& p1, const PolarPoint& p2) {
+Number<Inexact> Spiral::dAlphaDPhi1(const PolarPoint& p1, const PolarPoint& p2) {
 	const PolarPoint& source = p1.r() < p2.r() ? p2 : p1;
 	const PolarPoint& target = p1.r() < p2.r() ? p1 : p2;
 
-	if (source.r() == target.r()) {
-		throw std::runtime_error(
-		    "Cannot compute α for a spiral connecting two points equidistant to the root");
-	}
+	Number<Inexact> phiDiff = wrapAngle(target.phi() - source.phi(), -M_PI);
+	Number<Inexact> rDiffLog = std::log(target.r() / source.r());
+	return rDiffLog / (rDiffLog * rDiffLog + phiDiff * phiDiff);
+}
+
+Number<Inexact> Spiral::dAlphaDPhi2(const PolarPoint& p1, const PolarPoint& p2) {
+	return -dAlphaDPhi1(p1, p2);
+}
+
+Number<Inexact> Spiral::dAlphaDR1(const PolarPoint& p1, const PolarPoint& p2) {
+	const PolarPoint& source = p1.r() < p2.r() ? p2 : p1;
+	const PolarPoint& target = p1.r() < p2.r() ? p1 : p2;
 
 	Number<Inexact> phiDiff = wrapAngle(target.phi() - source.phi(), -M_PI);
 	Number<Inexact> rDiffLog = std::log(target.r() / source.r());
-	Number<Inexact> derivative = -rDiffLog / (phiDiff * phiDiff + rDiffLog * rDiffLog);
+	return (-phiDiff / p1.r()) / (rDiffLog * rDiffLog + phiDiff * phiDiff);
+}
 
-	// swap the sign if p1 was further from the origin than p2
-	int sign = p1.r() < p2.r() ? 1 : -1;
+Number<Inexact> Spiral::dAlphaDR2(const PolarPoint& p1, const PolarPoint& p2) {
+	const PolarPoint& source = p1.r() < p2.r() ? p2 : p1;
+	const PolarPoint& target = p1.r() < p2.r() ? p1 : p2;
 
-	return sign * derivative;
+	Number<Inexact> phiDiff = wrapAngle(target.phi() - source.phi(), -M_PI);
+	Number<Inexact> rDiffLog = std::log(target.r() / source.r());
+	return (phiDiff / p2.r()) / (rDiffLog * rDiffLog + phiDiff * phiDiff);
 }
 
 const PolarPoint& Spiral::anchor() const {
diff --git a/cartocrow/flow_map/spiral.h b/cartocrow/flow_map/spiral.h
index 7173d6557..1d7a91e8c 100644
--- a/cartocrow/flow_map/spiral.h
+++ b/cartocrow/flow_map/spiral.h
@@ -69,12 +69,55 @@ class Spiral {
 	Spiral(const PolarPoint& p1, const PolarPoint& p2);
 
 	/// Computes the \f$\alpha\f$ of the shortest logarithmic spiral connecting
-	/// the two given points.
-	static Number<Inexact> alphaBetweenPoints(const PolarPoint& p1, const PolarPoint& p2);
-	/// Computes \f$\frac{d\alpha}{d\phi}\f$, where \f$\alpha\f$ is the angle of
-	/// the shortest logarithmic spiral connecting the two given points (see
-	/// \ref alphaBetweenPoints) and \f$\phi\f$ is the polar angle of `p1`.
-	static Number<Inexact> alphaBetweenPointsDerivative(const PolarPoint& p1, const PolarPoint& p2);
+	/// the two given points \f$(r_1, \phi_1)\f$ and \f$(r_2, \phi_2)\f$.
+	///
+	/// This value can be computed as follows:
+	/// \f[
+	///     \alpha (r_1, \phi_1, r_2, \phi_2) = \arctan \left(
+	///         \frac{\phi_2 - \phi_1}{-\ln(r_2 / r_1)}
+	///     \right) \text{.}
+	/// \f]
+	static Number<Inexact> alpha(const PolarPoint& p1, const PolarPoint& p2);
+	/// Computes \f$\frac{\partial\alpha}{\partial r_1}\f$ for the \f$\alpha\f$
+	/// function (see \ref alpha).
+	///
+	/// This value can be computed as follows:
+	/// \f[
+	///     \frac{\partial\alpha}{\partial r_1}(r_1, \phi_1, r_2, \phi_2) =
+	///     \frac{-(\phi_2 - \phi_1) \ / \ r_1}{\ln^2(r_2 / r_1) + (\phi_2 - \phi_1)^2}
+	///     \text{.}
+	/// \f]
+	static Number<Inexact> dAlphaDR1(const PolarPoint& p1, const PolarPoint& p2);
+	/// Computes \f$\frac{\partial\alpha}{\partial r_2}\f$ for the \f$\alpha\f$
+	/// function (see \ref alpha).
+	///
+	/// This value can be computed as follows:
+	/// \f[
+	///     \frac{\partial\alpha}{\partial r_2}(r_1, \phi_1, r_2, \phi_2) =
+	///     \frac{(\phi_2 - \phi_1) \ / \ r_2}{\ln^2(r_2 / r_1) + (\phi_2 - \phi_1)^2}
+	///     \text{.}
+	/// \f]
+	static Number<Inexact> dAlphaDR2(const PolarPoint& p1, const PolarPoint& p2);
+	/// Computes \f$\frac{\partial\alpha}{\partial\phi_1}\f$ for the
+	/// \f$\alpha\f$ function (see \ref alpha).
+	///
+	/// This value can be computed as follows:
+	/// \f[
+	///     \frac{\partial\alpha}{\partial\phi_1}(r_1, \phi_1, r_2, \phi_2) =
+	///     \frac{\ln(r_2 / r_1)}{\ln^2(r_2 / r_1) + (\phi_2 - \phi_1)^2}
+	///     \text{.}
+	/// \f]
+	static Number<Inexact> dAlphaDPhi1(const PolarPoint& p1, const PolarPoint& p2);
+	/// Computes \f$\frac{\partial\alpha}{\partial\phi_2}\f$ for the
+	/// \f$\alpha\f$ function (see \ref alpha).
+	///
+	/// This value can be computed as follows:
+	/// \f[
+	///     \frac{\partial\alpha}{\partial\phi_2}(r_1, \phi_1, r_2, \phi_2) =
+	///     \frac{-\ln(r_2 / r_1)}{\ln^2(r_2 / r_1) + (\phi_2 - \phi_1)^2}
+	///     \text{.}
+	/// \f]
+	static Number<Inexact> dAlphaDPhi2(const PolarPoint& p1, const PolarPoint& p2);
 
 	/// Returns the anchor of this spiral.
 	const PolarPoint& anchor() const;

From fde296bbbca7de62a29d4b8c4c6c8596da2cc486 Mon Sep 17 00:00:00 2001
From: Willem Sonke <willem@wimiso.nl>
Date: Thu, 7 Sep 2023 13:12:23 +0200
Subject: [PATCH 15/27] Add the smoothing factor to the cost computation

---
 cartocrow/flow_map/smooth_tree.cpp | 29 ++++++++++++++---------------
 cartocrow/flow_map/smooth_tree.h   |  8 +++++++-
 2 files changed, 21 insertions(+), 16 deletions(-)

diff --git a/cartocrow/flow_map/smooth_tree.cpp b/cartocrow/flow_map/smooth_tree.cpp
index 285632a5d..5fb700636 100644
--- a/cartocrow/flow_map/smooth_tree.cpp
+++ b/cartocrow/flow_map/smooth_tree.cpp
@@ -73,10 +73,10 @@ const std::vector<std::shared_ptr<Node>>& SmoothTree::nodes() const {
 	return m_nodes;
 }
 
-Number<Inexact> SmoothTree::computeSmoothingFunction(const std::shared_ptr<Node>& node) {
-	PolarPoint n = node->m_position;
-	PolarPoint p = node->m_parent->m_position;
-	PolarPoint c = node->m_children[0]->m_position;
+Number<Inexact> SmoothTree::computeSmoothingCost(int i, int iParent, int iChild) {
+	PolarPoint n = m_nodes[i]->m_position;
+	PolarPoint p = m_nodes[iParent]->m_parent->m_position;
+	PolarPoint c = m_nodes[iChild]->m_children[0]->m_position;
 	return std::pow(Spiral::alpha(p, n) - Spiral::alpha(n, c), 2);
 }
 
@@ -85,21 +85,20 @@ void SmoothTree::applySmoothingForce(int i, int iParent, int iChild) {
 	PolarPoint p = m_nodes[iParent]->m_position;
 	PolarPoint c = m_nodes[iChild]->m_position;
 
-	// chain rule: the derivative of (β_1 - β_2)² is 2(β_1 - β_2) * [β_1 - β_2]'
-	m_forces[i].r += 2 * (Spiral::alpha(p, n) - Spiral::alpha(n, c)) *
+	m_forces[i].r += 2 * m_smoothing_factor * (Spiral::alpha(p, n) - Spiral::alpha(n, c)) *
 	                 (Spiral::dAlphaDR2(p, n) - Spiral::dAlphaDR1(n, c));
-	m_forces[i].phi += 2 * (Spiral::alpha(p, n) - Spiral::alpha(n, c)) *
+	m_forces[i].phi += 2 * m_smoothing_factor * (Spiral::alpha(p, n) - Spiral::alpha(n, c)) *
 	                   (Spiral::dAlphaDPhi2(p, n) - Spiral::dAlphaDPhi1(n, c));
 
-	m_forces[iParent].r += 2 * (Spiral::alpha(p, n) - Spiral::alpha(n, c)) *
-	                 (Spiral::dAlphaDR1(p, n));
-	m_forces[iParent].phi += 2 * (Spiral::alpha(p, n) - Spiral::alpha(n, c)) *
-	                   (Spiral::dAlphaDPhi1(p, n));
+	m_forces[iParent].r += 2 * m_smoothing_factor * (Spiral::alpha(p, n) - Spiral::alpha(n, c)) *
+	                       (Spiral::dAlphaDR1(p, n));
+	m_forces[iParent].phi += 2 * m_smoothing_factor * (Spiral::alpha(p, n) - Spiral::alpha(n, c)) *
+	                         (Spiral::dAlphaDPhi1(p, n));
 
-	m_forces[iChild].r += 2 * (Spiral::alpha(p, n) - Spiral::alpha(n, c)) *
-	                 (Spiral::dAlphaDR2(n, c));
-	m_forces[iChild].phi += 2 * (Spiral::alpha(p, n) - Spiral::alpha(n, c)) *
-	                   (Spiral::dAlphaDPhi2(n, c));
+	m_forces[iChild].r += 2 * m_smoothing_factor * (Spiral::alpha(p, n) - Spiral::alpha(n, c)) *
+	                      (Spiral::dAlphaDR2(n, c));
+	m_forces[iChild].phi += 2 * m_smoothing_factor * (Spiral::alpha(p, n) - Spiral::alpha(n, c)) *
+	                        (Spiral::dAlphaDPhi2(n, c));
 }
 
 void SmoothTree::optimize() {
diff --git a/cartocrow/flow_map/smooth_tree.h b/cartocrow/flow_map/smooth_tree.h
index 8ca316372..ea8b84ddc 100644
--- a/cartocrow/flow_map/smooth_tree.h
+++ b/cartocrow/flow_map/smooth_tree.h
@@ -60,8 +60,14 @@ class SmoothTree {
 	};
 	std::vector<PolarForce> m_forces;
 
-	Number<Inexact> computeSmoothingFunction(const std::shared_ptr<Node>& node);
+	/// Computes the smoothing cost
+	Number<Inexact> computeSmoothingCost(int i, int iParent, int iChild);
 	void applySmoothingForce(int i, int iParent, int iChild);
+
+	Number<Inexact> m_obstacle_factor = 2.0;
+	Number<Inexact> m_smoothing_factor = 0.4;
+	Number<Inexact> m_straightening_factor = 0.4;
+	Number<Inexact> m_angle_restriction_factor = 0.077;
 };
 
 } // namespace cartocrow::flow_map

From fd4b6518dfe37d29da0bbd8e0e5f799452ff8eb2 Mon Sep 17 00:00:00 2001
From: Willem Sonke <willem@wimiso.nl>
Date: Thu, 7 Sep 2023 14:15:40 +0200
Subject: [PATCH 16/27] Show the current cost function value in the demo

Also, document the formulas used for the smoothing cost.
---
 cartocrow/flow_map/smooth_tree.cpp            | 45 ++++++++++++++-----
 cartocrow/flow_map/smooth_tree.h              | 43 +++++++++++++++++-
 .../optimization_demo/optimization_demo.cpp   | 14 +++++-
 .../optimization_demo/optimization_demo.h     |  2 +
 .../spiral_tree_demo/spiral_tree_demo.cpp     |  1 +
 5 files changed, 92 insertions(+), 13 deletions(-)

diff --git a/cartocrow/flow_map/smooth_tree.cpp b/cartocrow/flow_map/smooth_tree.cpp
index 5fb700636..8062eaf6c 100644
--- a/cartocrow/flow_map/smooth_tree.cpp
+++ b/cartocrow/flow_map/smooth_tree.cpp
@@ -33,6 +33,11 @@ SmoothTree::SmoothTree(const std::shared_ptr<SpiralTree> spiralTree) : m_tree(sp
 		}
 	}
 	constructSmoothTree(spiralTree->root(), rMax / 10); // TODO make 10 configurable
+
+	// set node IDs
+	for (int i = 0; i < m_nodes.size(); i++) {
+		m_nodes[i]->m_id = i;
+	}
 }
 
 std::shared_ptr<Node> SmoothTree::constructSmoothTree(const std::shared_ptr<Node>& node, Number<Inexact> maxRStep) {
@@ -75,9 +80,9 @@ const std::vector<std::shared_ptr<Node>>& SmoothTree::nodes() const {
 
 Number<Inexact> SmoothTree::computeSmoothingCost(int i, int iParent, int iChild) {
 	PolarPoint n = m_nodes[i]->m_position;
-	PolarPoint p = m_nodes[iParent]->m_parent->m_position;
-	PolarPoint c = m_nodes[iChild]->m_children[0]->m_position;
-	return std::pow(Spiral::alpha(p, n) - Spiral::alpha(n, c), 2);
+	PolarPoint p = m_nodes[iParent]->m_position;
+	PolarPoint c = m_nodes[iChild]->m_position;
+	return m_smoothing_factor * std::pow(Spiral::alpha(p, n) - Spiral::alpha(n, c), 2);
 }
 
 void SmoothTree::applySmoothingForce(int i, int iParent, int iChild) {
@@ -91,29 +96,47 @@ void SmoothTree::applySmoothingForce(int i, int iParent, int iChild) {
 	                   (Spiral::dAlphaDPhi2(p, n) - Spiral::dAlphaDPhi1(n, c));
 
 	m_forces[iParent].r += 2 * m_smoothing_factor * (Spiral::alpha(p, n) - Spiral::alpha(n, c)) *
-	                       (Spiral::dAlphaDR1(p, n));
+	                       Spiral::dAlphaDR1(p, n);
 	m_forces[iParent].phi += 2 * m_smoothing_factor * (Spiral::alpha(p, n) - Spiral::alpha(n, c)) *
-	                         (Spiral::dAlphaDPhi1(p, n));
+	                         Spiral::dAlphaDPhi1(p, n);
 
 	m_forces[iChild].r += 2 * m_smoothing_factor * (Spiral::alpha(p, n) - Spiral::alpha(n, c)) *
-	                      (Spiral::dAlphaDR2(n, c));
+	                      -Spiral::dAlphaDR2(n, c);
 	m_forces[iChild].phi += 2 * m_smoothing_factor * (Spiral::alpha(p, n) - Spiral::alpha(n, c)) *
-	                        (Spiral::dAlphaDPhi2(n, c));
+	                        -Spiral::dAlphaDPhi2(n, c);
 }
 
-void SmoothTree::optimize() {
-	m_forces = std::vector<PolarForce>(m_nodes.size(), PolarForce{});
+Number<Inexact> SmoothTree::computeAngleRestrictionCost(int i, int iParent, int iChild) {
+	return 0; // TODO
+}
+
+void SmoothTree::applyAngleRestrictionForce(int i, int iParent, int iChild) {
+	// TODO
+}
+
+Number<Inexact> SmoothTree::computeCost() {
+	Number<Inexact> cost = 0;
 	for (int i = 0; i < m_nodes.size(); i++) {
-		m_nodes[i]->m_id = i;
+		const auto& node = m_nodes[i];
+		if (node->getType() == Node::ConnectionType::kSubdivision) {
+			cost += computeSmoothingCost(i, m_nodes[i]->m_parent->m_id, m_nodes[i]->m_children[0]->m_id);
+			cost += computeAngleRestrictionCost(i, m_nodes[i]->m_parent->m_id, m_nodes[i]->m_children[0]->m_id);
+		}
 	}
+	return cost;
+}
+
+void SmoothTree::optimize() {
+	m_forces = std::vector<PolarForce>(m_nodes.size(), PolarForce{});
 	for (int i = 0; i < m_nodes.size(); i++) {
 		const auto& node = m_nodes[i];
 		if (node->getType() == Node::ConnectionType::kSubdivision) {
 			applySmoothingForce(i, m_nodes[i]->m_parent->m_id, m_nodes[i]->m_children[0]->m_id);
+			applyAngleRestrictionForce(i, m_nodes[i]->m_parent->m_id, m_nodes[i]->m_children[0]->m_id);
 		}
 	}
 	for (int i = 0; i < m_nodes.size(); i++) {
-		Number<Inexact> epsilon = 0.001;
+		Number<Inexact> epsilon = 0.001; // TODO
 		if (m_nodes[i]->getType() == Node::ConnectionType::kJoin) {
 			m_nodes[i]->m_position.setR(m_nodes[i]->m_position.r() + epsilon * m_forces[i].r);
 		}
diff --git a/cartocrow/flow_map/smooth_tree.h b/cartocrow/flow_map/smooth_tree.h
index ea8b84ddc..72acc3664 100644
--- a/cartocrow/flow_map/smooth_tree.h
+++ b/cartocrow/flow_map/smooth_tree.h
@@ -42,6 +42,9 @@ class SmoothTree {
 	/// Returns a list of the nodes in this smooth tree.
 	const std::vector<std::shared_ptr<Node>>& nodes() const;
 
+	/// Computes the total cost of the tree.
+	Number<Inexact> computeCost();
+
 	/// Performs one optimization step.
 	void optimize();
 
@@ -60,9 +63,47 @@ class SmoothTree {
 	};
 	std::vector<PolarForce> m_forces;
 
-	/// Computes the smoothing cost
+	/// Computes the smoothing cost for the subdivision node `i` at \f$(r,
+	/// \phi)\f$, with parent `iParent` at \f$(r_p, \phi_p)\f$ and child
+	/// `iChild` at \f$(r_c, \phi_c)\f$.
+	///
+	/// \f[
+	///     F_\text{smooth}(r, \phi, r_p, \phi_p, r_c, \phi_c) =
+	///     c_\text{smooth} \cdot
+	///     \big( \alpha(r_p, \phi_p, r, \phi)
+	///         - \alpha(r, \phi, r_c, \phi_c) \big)^2
+	///     \text{.}
+	/// \f]
 	Number<Inexact> computeSmoothingCost(int i, int iParent, int iChild);
+	/// Applies smoothing forces in \ref m_forces to the subdivision node `i`,
+	/// its parent `iParent`, and its child `iChild`.
+	///
+	/// The forces are defined by the partial derivatives of the smoothing cost
+	/// (see \ref computeSmoothingCost), which are:
+	///
+	/// \f{align*}{
+	///     \frac{\partial F_\text{smooth}}{\partial r}(r, \phi, r_p, \phi_p, r_c, \phi_c)\ &=
+	///     2c_\text{smooth} \cdot
+	///     \big( \alpha(r_p, \phi_p, r, \phi) - \alpha(r, \phi, r_c, \phi_c) \big) \cdot
+	///     \Big( \frac{\partial\alpha}{\partial r_2}(r_p, \phi_p, r, \phi)
+	///         - \frac{\partial\alpha}{\partial r_1}(r, \phi, r_c, \phi_c) \Big)
+	///     \text{;} \\
+	///     \frac{\partial F_\text{smooth}}{\partial r_p}(r, \phi, r_p, \phi_p, r_c, \phi_c)\ &=
+	///     2c_\text{smooth} \cdot
+	///     \big( \alpha(r_p, \phi_p, r, \phi) - \alpha(r, \phi, r_c, \phi_c) \big) \cdot
+	///     \frac{\partial\alpha}{\partial r_1}(r_p, \phi_p, r, \phi)
+	///     \text{;} \\
+	///     \frac{\partial F_\text{smooth}}{\partial r_c}(r, \phi, r_p, \phi_p, r_c, \phi_c)\ &=
+	///     2c_\text{smooth} \cdot
+	///     \big( \alpha(r_p, \phi_p, r, \phi) - \alpha(r, \phi, r_c, \phi_c) \big) \cdot
+	///     -\frac{\partial\alpha}{\partial r_2}(r, \phi, r_c, \phi_c)
+	///     \text{;} \\
+	/// \f}
+	/// et cetera.
 	void applySmoothingForce(int i, int iParent, int iChild);
+	
+	Number<Inexact> computeAngleRestrictionCost(int i, int iParent, int iChild);
+	void applyAngleRestrictionForce(int i, int iParent, int iChild);
 
 	Number<Inexact> m_obstacle_factor = 2.0;
 	Number<Inexact> m_smoothing_factor = 0.4;
diff --git a/demos/flow_map/optimization_demo/optimization_demo.cpp b/demos/flow_map/optimization_demo/optimization_demo.cpp
index fd5ab2492..20468f0cb 100644
--- a/demos/flow_map/optimization_demo/optimization_demo.cpp
+++ b/demos/flow_map/optimization_demo/optimization_demo.cpp
@@ -23,7 +23,7 @@ along with this program.  If not, see <https://www.gnu.org/licenses/>.
 
 #include <QApplication>
 #include <QCheckBox>
-#include <cstdlib>
+#include <QStatusBar>
 
 #include "cartocrow/core/core.h"
 #include "cartocrow/flow_map/painting.h"
@@ -80,12 +80,14 @@ OptimizationDemo::OptimizationDemo() {
 	connect(m_optimizeTimer, &QTimer::timeout, [&]() {
 		m_smoothTree->optimize();
 		m_renderer->update();
+		updateCostLabel();
 	});
 	m_optimizeOneStepButton = new QPushButton("Optimize one step");
 	toolBar->addWidget(m_optimizeOneStepButton);
 	connect(m_optimizeOneStepButton, &QPushButton::clicked, [&]() {
 		m_smoothTree->optimize();
 		m_renderer->update();
+		updateCostLabel();
 	});
 	toolBar->addSeparator();
 	toolBar->addWidget(new QLabel("α = "));
@@ -102,6 +104,10 @@ OptimizationDemo::OptimizationDemo() {
 		m_alphaLabel->setText(QString("%1π").arg(value / 1000.0, 0, 'f', 3));
 		recalculate();
 	});
+
+	m_costLabel = new QLabel();
+	statusBar()->addWidget(m_costLabel);
+
 	for (auto place : m_places) {
 		m_renderer->registerEditable(place);
 	}
@@ -134,8 +140,14 @@ void OptimizationDemo::recalculate() {
 	m_renderer->addPainting(painting, "Smooth tree");
 
 	m_renderer->update();
+	updateCostLabel();
 }
 
+void OptimizationDemo::updateCostLabel() {
+	m_costLabel->setText(QString("Cost function: %1").arg(m_smoothTree->computeCost()));
+}
+
+
 int main(int argc, char* argv[]) {
 	QApplication app(argc, argv);
 	OptimizationDemo demo;
diff --git a/demos/flow_map/optimization_demo/optimization_demo.h b/demos/flow_map/optimization_demo/optimization_demo.h
index e02b36f3e..533c2c26b 100644
--- a/demos/flow_map/optimization_demo/optimization_demo.h
+++ b/demos/flow_map/optimization_demo/optimization_demo.h
@@ -43,6 +43,7 @@ class OptimizationDemo : public QMainWindow {
 
   private:
 	void recalculate();
+	void updateCostLabel();
 	Number<Inexact> m_alpha = 25 * M_PI / 180;
 
 	std::vector<std::shared_ptr<Point<Inexact>>> m_places;
@@ -54,4 +55,5 @@ class OptimizationDemo : public QMainWindow {
 	QPushButton* m_optimizeButton;
 	QPushButton* m_optimizeOneStepButton;
 	QTimer* m_optimizeTimer;
+	QLabel* m_costLabel;
 };
diff --git a/demos/flow_map/spiral_tree_demo/spiral_tree_demo.cpp b/demos/flow_map/spiral_tree_demo/spiral_tree_demo.cpp
index e9b349040..6d8f57df1 100644
--- a/demos/flow_map/spiral_tree_demo/spiral_tree_demo.cpp
+++ b/demos/flow_map/spiral_tree_demo/spiral_tree_demo.cpp
@@ -23,6 +23,7 @@ along with this program.  If not, see <https://www.gnu.org/licenses/>.
 
 #include <QApplication>
 #include <QCheckBox>
+#include <QMainWindow>
 
 #include "cartocrow/core/core.h"
 #include "cartocrow/core/timer.h"

From 99b34bcc6381a8ed62d90a51b06047bacdd24d8c Mon Sep 17 00:00:00 2001
From: Willem Sonke <willem@wimiso.nl>
Date: Thu, 7 Sep 2023 14:38:34 +0200
Subject: [PATCH 17/27] Implement angle restriction cost

---
 cartocrow/flow_map/smooth_tree.cpp | 41 +++++++++++++++++++++++++-----
 cartocrow/flow_map/smooth_tree.h   | 36 ++++++++++++++++++++++++--
 2 files changed, 69 insertions(+), 8 deletions(-)

diff --git a/cartocrow/flow_map/smooth_tree.cpp b/cartocrow/flow_map/smooth_tree.cpp
index 8062eaf6c..48677aba4 100644
--- a/cartocrow/flow_map/smooth_tree.cpp
+++ b/cartocrow/flow_map/smooth_tree.cpp
@@ -106,12 +106,35 @@ void SmoothTree::applySmoothingForce(int i, int iParent, int iChild) {
 	                        -Spiral::dAlphaDPhi2(n, c);
 }
 
-Number<Inexact> SmoothTree::computeAngleRestrictionCost(int i, int iParent, int iChild) {
-	return 0; // TODO
+Number<Inexact> SmoothTree::computeAngleRestrictionCost(int i, int iChild1, int iChild2) {
+	PolarPoint n = m_nodes[i]->m_position;
+	PolarPoint c1 = m_nodes[iChild1]->m_position;
+	PolarPoint c2 = m_nodes[iChild2]->m_position;
+	return m_angle_restriction_factor * (std::log(1.0 / std::cos(Spiral::alpha(n, c1))) +
+	                                     std::log(1.0 / std::cos(Spiral::alpha(n, c2))));
 }
 
-void SmoothTree::applyAngleRestrictionForce(int i, int iParent, int iChild) {
-	// TODO
+void SmoothTree::applyAngleRestrictionForce(int i, int iChild1, int iChild2) {
+	PolarPoint n = m_nodes[i]->m_position;
+	PolarPoint c1 = m_nodes[iChild1]->m_position;
+	PolarPoint c2 = m_nodes[iChild2]->m_position;
+
+	m_forces[i].r +=
+	    m_angle_restriction_factor * (Spiral::dAlphaDR1(n, c1) * std::tan(Spiral::alpha(n, c1)) +
+	                                  Spiral::dAlphaDR1(n, c2) * std::tan(Spiral::alpha(n, c2)));
+	m_forces[i].phi +=
+	    m_angle_restriction_factor * (Spiral::dAlphaDPhi1(n, c1) * std::tan(Spiral::alpha(n, c1)) +
+	                                  Spiral::dAlphaDPhi1(n, c2) * std::tan(Spiral::alpha(n, c2)));
+
+	m_forces[iChild1].r +=
+	    m_angle_restriction_factor * Spiral::dAlphaDR2(n, c1) * std::tan(Spiral::alpha(n, c1));
+	m_forces[iChild1].phi +=
+	    m_angle_restriction_factor * Spiral::dAlphaDPhi2(n, c1) * std::tan(Spiral::alpha(n, c1));
+
+	m_forces[iChild2].r +=
+	    m_angle_restriction_factor * Spiral::dAlphaDR2(n, c2) * std::tan(Spiral::alpha(n, c2));
+	m_forces[iChild2].phi +=
+	    m_angle_restriction_factor * Spiral::dAlphaDPhi2(n, c2) * std::tan(Spiral::alpha(n, c2));
 }
 
 Number<Inexact> SmoothTree::computeCost() {
@@ -120,7 +143,10 @@ Number<Inexact> SmoothTree::computeCost() {
 		const auto& node = m_nodes[i];
 		if (node->getType() == Node::ConnectionType::kSubdivision) {
 			cost += computeSmoothingCost(i, m_nodes[i]->m_parent->m_id, m_nodes[i]->m_children[0]->m_id);
-			cost += computeAngleRestrictionCost(i, m_nodes[i]->m_parent->m_id, m_nodes[i]->m_children[0]->m_id);
+		} else if (node->getType() == Node::ConnectionType::kJoin) {
+			cost += computeAngleRestrictionCost(
+			    i, m_nodes[i]->m_children[0]->m_id,
+			    m_nodes[i]->m_children[m_nodes[i]->m_children.size() - 1]->m_id);
 		}
 	}
 	return cost;
@@ -132,7 +158,10 @@ void SmoothTree::optimize() {
 		const auto& node = m_nodes[i];
 		if (node->getType() == Node::ConnectionType::kSubdivision) {
 			applySmoothingForce(i, m_nodes[i]->m_parent->m_id, m_nodes[i]->m_children[0]->m_id);
-			applyAngleRestrictionForce(i, m_nodes[i]->m_parent->m_id, m_nodes[i]->m_children[0]->m_id);
+		} else if (node->getType() == Node::ConnectionType::kJoin) {
+			applyAngleRestrictionForce(
+			    i, m_nodes[i]->m_children[0]->m_id,
+			    m_nodes[i]->m_children[m_nodes[i]->m_children.size() - 1]->m_id);
 		}
 	}
 	for (int i = 0; i < m_nodes.size(); i++) {
diff --git a/cartocrow/flow_map/smooth_tree.h b/cartocrow/flow_map/smooth_tree.h
index 72acc3664..f9bc216ed 100644
--- a/cartocrow/flow_map/smooth_tree.h
+++ b/cartocrow/flow_map/smooth_tree.h
@@ -102,8 +102,40 @@ class SmoothTree {
 	/// et cetera.
 	void applySmoothingForce(int i, int iParent, int iChild);
 	
-	Number<Inexact> computeAngleRestrictionCost(int i, int iParent, int iChild);
-	void applyAngleRestrictionForce(int i, int iParent, int iChild);
+	/// Computes the angle restriction cost for the join node `i` at \f$(r,
+	/// \phi)\f$, with children `iChild1` at \f$(r_{c_1}, \phi_{c_1})\f$ and
+	/// `iChild2` at \f$(r_{c_2}, \phi_{c_2})\f$.
+	///
+	/// \f[
+	///     F_\text{AR}(r, \phi, r_{c_1}, \phi_{c_1}, r_{c_2}, \phi_{c_2}) =
+	///     c_\text{AR} \cdot
+	///     \big( \log(\sec \alpha(r, \phi, r_{c_1}, \phi_{c_1}))
+	///         + \log(\sec \alpha(r, \phi, r_{c_2}, \phi_{c_2})) \big)
+	///     \text{.}
+	/// \f]
+	Number<Inexact> computeAngleRestrictionCost(int i, int iChild1, int iChild2);
+	/// Applies angle restriction forces in \ref m_forces to the subdivision
+	/// node `i` and its children `iChild1` and `iChild2`.
+	///
+	/// The forces are defined by the partial derivatives of the angle
+	/// restriction cost (see \ref computeAngleRestrictionCost), which are:
+	///
+	/// \f{align*}{
+	///     \frac{\partial F_\text{AR}}{\partial r}(r, \phi, r_{c_1}, \phi_{c_1}, r_{c_2}, \phi_{c_2})\ &=
+	///     c_\text{AR} \cdot
+	///     \Big( \frac{\partial\alpha}{\partial r_1}(r, \phi, r_{c_1}, \phi_{c_1}) \cdot
+	///             \tan \alpha(r, \phi, r_{c_1}, \phi_{c_1}) +
+	///         \frac{\partial\alpha}{\partial r_1}(r, \phi, r_{c_2}, \phi_{c_2}) \cdot
+	///             \tan \alpha(r, \phi, r_{c_2}, \phi_{c_2}) \Big)
+	///     \text{;} \\
+	///     \frac{\partial F_\text{AR}}{\partial r_{c_1}}(r, \phi, r_{c_1}, \phi_{c_1}, r_{c_2}, \phi_{c_2})\ &=
+	///     c_\text{AR} \cdot
+	///     \frac{\partial\alpha}{\partial r_2}(r, \phi, r_{c_1}, \phi_{c_1}) \cdot
+	///         \tan \alpha(r, \phi, r_{c_1}, \phi_{c_1})
+	///     \text{;} \\
+	/// \f}
+	/// et cetera.
+	void applyAngleRestrictionForce(int i, int iChild1, int iChild2);
 
 	Number<Inexact> m_obstacle_factor = 2.0;
 	Number<Inexact> m_smoothing_factor = 0.4;

From dbf331c2f626623252db718a53195c9a93b01540 Mon Sep 17 00:00:00 2001
From: Willem Sonke <willem@wimiso.nl>
Date: Fri, 8 Sep 2023 19:36:19 +0200
Subject: [PATCH 18/27] Implement balancing cost

---
 cartocrow/flow_map/smooth_tree.cpp | 47 ++++++++++++++++++++++++++++--
 cartocrow/flow_map/smooth_tree.h   | 40 +++++++++++++++++++++++++
 2 files changed, 85 insertions(+), 2 deletions(-)

diff --git a/cartocrow/flow_map/smooth_tree.cpp b/cartocrow/flow_map/smooth_tree.cpp
index 48677aba4..86210fc56 100644
--- a/cartocrow/flow_map/smooth_tree.cpp
+++ b/cartocrow/flow_map/smooth_tree.cpp
@@ -25,7 +25,8 @@ along with this program.  If not, see <https://www.gnu.org/licenses/>.
 
 namespace cartocrow::flow_map {
 
-SmoothTree::SmoothTree(const std::shared_ptr<SpiralTree> spiralTree) : m_tree(spiralTree) {
+SmoothTree::SmoothTree(const std::shared_ptr<SpiralTree> spiralTree)
+    : m_tree(spiralTree), m_restrictingAngle(spiralTree->restrictingAngle()) {
 	Number<Inexact> rMax = 0;
 	for (const auto& node : m_tree->nodes()) {
 		if (node->m_position.r() > rMax) {
@@ -137,6 +138,42 @@ void SmoothTree::applyAngleRestrictionForce(int i, int iChild1, int iChild2) {
 	    m_angle_restriction_factor * Spiral::dAlphaDPhi2(n, c2) * std::tan(Spiral::alpha(n, c2));
 }
 
+Number<Inexact> SmoothTree::computeBalancingCost(int i, int iChild1, int iChild2) {
+	PolarPoint n = m_nodes[i]->m_position;
+	PolarPoint c1 = m_nodes[iChild1]->m_position;
+	PolarPoint c2 = m_nodes[iChild2]->m_position;
+	return m_angle_restriction_factor * 2 * std::pow(std::tan(m_restrictingAngle), 2) *
+	       std::log(1 / std::sin(0.5 * (Spiral::alpha(n, c1) - Spiral::alpha(n, c2))));
+}
+
+void SmoothTree::applyBalancingForce(int i, int iChild1, int iChild2) {
+	PolarPoint n = m_nodes[i]->m_position;
+	PolarPoint c1 = m_nodes[iChild1]->m_position;
+	PolarPoint c2 = m_nodes[iChild2]->m_position;
+
+	m_forces[i].r += m_angle_restriction_factor * -std::pow(std::tan(m_restrictingAngle), 2) *
+	                 (1 / std::tan(0.5 * (Spiral::alpha(n, c1) - Spiral::alpha(n, c2)))) *
+	                 (Spiral::dAlphaDR1(n, c1) - Spiral::dAlphaDR1(n, c2));
+	m_forces[i].phi += m_angle_restriction_factor * -std::pow(std::tan(m_restrictingAngle), 2) *
+	                 (1 / std::tan(0.5 * (Spiral::alpha(n, c1) - Spiral::alpha(n, c2)))) *
+	                 (Spiral::dAlphaDPhi1(n, c1) - Spiral::dAlphaDPhi1(n, c2));
+
+	m_forces[iChild1].r += m_angle_restriction_factor * -std::pow(std::tan(m_restrictingAngle), 2) *
+	                       (1 / std::tan(0.5 * (Spiral::alpha(n, c1) - Spiral::alpha(n, c2)))) *
+	                       Spiral::dAlphaDR2(n, c1);
+	m_forces[iChild1].phi += m_angle_restriction_factor * -std::pow(std::tan(m_restrictingAngle), 2) *
+	                         (1 / std::tan(0.5 * (Spiral::alpha(n, c1) - Spiral::alpha(n, c2)))) *
+	                         Spiral::dAlphaDPhi2(n, c1);
+
+	m_forces[iChild2].r += m_angle_restriction_factor * -std::pow(std::tan(m_restrictingAngle), 2) *
+	                       (1 / std::tan(0.5 * (Spiral::alpha(n, c1) - Spiral::alpha(n, c2)))) *
+	                       -Spiral::dAlphaDR2(n, c2);
+	m_forces[iChild2].phi += m_angle_restriction_factor * -std::pow(std::tan(m_restrictingAngle), 2) *
+	                         (1 / std::tan(0.5 * (Spiral::alpha(n, c1) - Spiral::alpha(n, c2)))) *
+	                         -Spiral::dAlphaDPhi2(n, c2);
+}
+
+
 Number<Inexact> SmoothTree::computeCost() {
 	Number<Inexact> cost = 0;
 	for (int i = 0; i < m_nodes.size(); i++) {
@@ -147,6 +184,9 @@ Number<Inexact> SmoothTree::computeCost() {
 			cost += computeAngleRestrictionCost(
 			    i, m_nodes[i]->m_children[0]->m_id,
 			    m_nodes[i]->m_children[m_nodes[i]->m_children.size() - 1]->m_id);
+			cost += computeBalancingCost(
+			    i, m_nodes[i]->m_children[0]->m_id,
+			    m_nodes[i]->m_children[m_nodes[i]->m_children.size() - 1]->m_id);
 		}
 	}
 	return cost;
@@ -162,12 +202,15 @@ void SmoothTree::optimize() {
 			applyAngleRestrictionForce(
 			    i, m_nodes[i]->m_children[0]->m_id,
 			    m_nodes[i]->m_children[m_nodes[i]->m_children.size() - 1]->m_id);
+			applyBalancingForce(
+			    i, m_nodes[i]->m_children[0]->m_id,
+			    m_nodes[i]->m_children[m_nodes[i]->m_children.size() - 1]->m_id);
 		}
 	}
 	for (int i = 0; i < m_nodes.size(); i++) {
 		Number<Inexact> epsilon = 0.001; // TODO
 		if (m_nodes[i]->getType() == Node::ConnectionType::kJoin) {
-			m_nodes[i]->m_position.setR(m_nodes[i]->m_position.r() + epsilon * m_forces[i].r);
+			//m_nodes[i]->m_position.setR(m_nodes[i]->m_position.r() + epsilon * m_forces[i].r);
 		}
 		if (m_nodes[i]->getType() == Node::ConnectionType::kJoin ||
 		    m_nodes[i]->getType() == Node::ConnectionType::kSubdivision) {
diff --git a/cartocrow/flow_map/smooth_tree.h b/cartocrow/flow_map/smooth_tree.h
index f9bc216ed..7472d6850 100644
--- a/cartocrow/flow_map/smooth_tree.h
+++ b/cartocrow/flow_map/smooth_tree.h
@@ -52,6 +52,8 @@ class SmoothTree {
 	/// The spiral tree underlying this smooth tree.
 	std::shared_ptr<SpiralTree> m_tree;
 
+	Number<Inexact> m_restrictingAngle;
+
 	/// List of nodes in this tree.
 	std::vector<std::shared_ptr<Node>> m_nodes;
 
@@ -136,6 +138,44 @@ class SmoothTree {
 	/// \f}
 	/// et cetera.
 	void applyAngleRestrictionForce(int i, int iChild1, int iChild2);
+	/// Computes the balancing cost for the join node `i` at \f$(r, \phi)\f$,
+	/// with children `iChild1` at \f$(r_{c_1}, \phi_{c_1})\f$ and `iChild2` at
+	/// \f$(r_{c_2}, \phi_{c_2})\f$.
+	///
+	/// \f[
+	///     F_\text{balance}(r, \phi, r_{c_1}, \phi_{c_1}, r_{c_2}, \phi_{c_2}) =
+	///     c_\text{AR} \cdot
+	///     2 \tan^2(\alpha) \log\Big( \csc \Big(
+	///         \frac{\alpha(r, \phi, r_{c_1}, \phi_{c_1}) - \alpha(r, \phi, r_{c_2}, \phi_{c_2})}{2}
+	///     \Big) \Big)
+	///     \text{.}
+	/// \f]
+	Number<Inexact> computeBalancingCost(int i, int iChild1, int iChild2);
+	/// Applies balancing forces in \ref m_forces to the subdivision node `i`
+	/// and its children `iChild1` and `iChild2`.
+	///
+	/// The forces are defined by the partial derivatives of the balancing cost
+	/// (see \ref computeBalancingCost), which are:
+	///
+	/// \f{multline}{
+	///     \frac{\partial F_\text{balance}}{\partial r}(r, \phi, r_{c_1}, \phi_{c_1}, r_{c_2}, \phi_{c_2})\ =
+	///     c_\text{AR} \cdot -\tan^2(\alpha) \cdot
+	///     \cot \Big(
+	///         \frac{\alpha(r, \phi, r_{c_1}, \phi_{c_1}) - \alpha(r, \phi, r_{c_2}, \phi_{c_2})}{2}
+	///     \Big) \\
+	///     \cdot \Big( \frac{\partial\alpha}{\partial r_1}(r, \phi, r_{c_1}, \phi_{c_1}) -
+	///         \frac{\partial\alpha}{\partial r_1}(r, \phi, r_{c_2}, \phi_{c_2}) \Big)
+	///     \text{;} \\
+	/// \f}
+	/// et cetera.
+	void applyBalancingForce(int i, int iChild1, int iChild2);
+
+	// [ws] TODO Somewhere in the optimization code there seems to be a sign
+	// issue, although the result seems correct. We are computing forces by
+	// computing the gradient, but the forces are supposed to point towards the
+	// negative gradient, not the positive. Yet, we are never inverting the
+	// sign. I'm confused as to why this seems to be working fine; maybe the
+	// sign of the α function is incorrect?
 
 	Number<Inexact> m_obstacle_factor = 2.0;
 	Number<Inexact> m_smoothing_factor = 0.4;

From 31b6650024631417ada84ab2f9eea6416980ce92 Mon Sep 17 00:00:00 2001
From: Willem Sonke <willem@wimiso.nl>
Date: Fri, 8 Sep 2023 21:30:14 +0200
Subject: [PATCH 19/27] Implement straightening cost

---
 cartocrow/flow_map/smooth_tree.cpp | 76 +++++++++++++++++++++++++++++-
 cartocrow/flow_map/smooth_tree.h   | 47 ++++++++++++++++--
 cartocrow/flow_map/spiral.cpp      | 12 +++--
 3 files changed, 127 insertions(+), 8 deletions(-)

diff --git a/cartocrow/flow_map/smooth_tree.cpp b/cartocrow/flow_map/smooth_tree.cpp
index 86210fc56..cb4c218bf 100644
--- a/cartocrow/flow_map/smooth_tree.cpp
+++ b/cartocrow/flow_map/smooth_tree.cpp
@@ -173,6 +173,74 @@ void SmoothTree::applyBalancingForce(int i, int iChild1, int iChild2) {
 	                         -Spiral::dAlphaDPhi2(n, c2);
 }
 
+Number<Inexact>
+SmoothTree::computeStraighteningCost(int i, int iParent,
+                                     const std::vector<std::shared_ptr<Node>>& children) {
+	PolarPoint n = m_nodes[i]->m_position;
+	PolarPoint p = m_nodes[iParent]->m_position;
+	Number<Inexact> maxFlow = 0;
+	for (const auto& child : children) {
+		if (child->m_flow > maxFlow) {
+			maxFlow = child->m_flow;
+		}
+	}
+	Number<Inexact> numerator = 0;
+	Number<Inexact> denominator = 0;
+	for (const auto& child : children) {
+		if (child->m_flow >= m_relevantFlowFactor * maxFlow) {
+			PolarPoint c = child->m_position;
+			numerator += child->m_flow * Spiral::alpha(n, c);
+			denominator += child->m_flow;
+		}
+	}
+	return m_straightening_factor * std::pow(Spiral::alpha(p, n) - numerator / denominator, 2);
+}
+
+void SmoothTree::applyStraighteningForce(int i, int iParent,
+                                         const std::vector<std::shared_ptr<Node>>& children) {
+	PolarPoint n = m_nodes[i]->m_position;
+	PolarPoint p = m_nodes[iParent]->m_position;
+	Number<Inexact> maxFlow = 0;
+	for (const auto& child : children) {
+		if (child->m_flow > maxFlow) {
+			maxFlow = child->m_flow;
+		}
+	}
+	Number<Inexact> numerator = 0;
+	Number<Inexact> numeratorDR1 = 0;
+	Number<Inexact> numeratorDPhi1 = 0;
+	Number<Inexact> denominator = 0;
+	for (const auto& child : children) {
+		if (child->m_flow >= m_relevantFlowFactor * maxFlow) {
+			PolarPoint c = child->m_position;
+			numerator += child->m_flow * Spiral::alpha(n, c);
+			numeratorDR1 += child->m_flow * Spiral::dAlphaDR1(n, c);
+			numeratorDPhi1 += child->m_flow * Spiral::dAlphaDPhi1(n, c);
+			denominator += child->m_flow;
+		}
+	}
+	m_forces[i].r += 2 * m_straightening_factor * (Spiral::alpha(p, n) - numerator / denominator) *
+	                 (Spiral::dAlphaDR2(p, n) - numeratorDR1 / denominator);
+	m_forces[i].phi += 2 * m_straightening_factor * (Spiral::alpha(p, n) - numerator / denominator) *
+	                   (Spiral::dAlphaDPhi2(p, n) - numeratorDPhi1 / denominator);
+
+	m_forces[iParent].r += 2 * m_straightening_factor * (Spiral::alpha(p, n) - numerator / denominator) *
+	                 Spiral::dAlphaDR1(p, n);
+	m_forces[iParent].phi += 2 * m_straightening_factor * (Spiral::alpha(p, n) - numerator / denominator) *
+	                   Spiral::dAlphaDPhi1(p, n);
+
+	for (const auto& child : children) {
+		if (child->m_flow > maxFlow) {
+			PolarPoint c = child->m_position;
+			m_forces[child->m_id].r += 2 * m_straightening_factor *
+			                           (Spiral::alpha(p, n) - numerator / denominator) *
+			                           -child->m_flow * Spiral::dAlphaDR2(n, c) / denominator;
+			m_forces[child->m_id].phi += 2 * m_straightening_factor *
+			                           (Spiral::alpha(p, n) - numerator / denominator) *
+			                           -child->m_flow * Spiral::dAlphaDPhi2(n, c) / denominator;
+		}
+	}
+}
 
 Number<Inexact> SmoothTree::computeCost() {
 	Number<Inexact> cost = 0;
@@ -187,6 +255,9 @@ Number<Inexact> SmoothTree::computeCost() {
 			cost += computeBalancingCost(
 			    i, m_nodes[i]->m_children[0]->m_id,
 			    m_nodes[i]->m_children[m_nodes[i]->m_children.size() - 1]->m_id);
+			cost += computeStraighteningCost(
+			    i, m_nodes[i]->m_parent->m_id,
+			    m_nodes[i]->m_children);
 		}
 	}
 	return cost;
@@ -205,10 +276,13 @@ void SmoothTree::optimize() {
 			applyBalancingForce(
 			    i, m_nodes[i]->m_children[0]->m_id,
 			    m_nodes[i]->m_children[m_nodes[i]->m_children.size() - 1]->m_id);
+			applyStraighteningForce(
+			    i, m_nodes[i]->m_parent->m_id,
+			    m_nodes[i]->m_children);
 		}
 	}
 	for (int i = 0; i < m_nodes.size(); i++) {
-		Number<Inexact> epsilon = 0.001; // TODO
+		Number<Inexact> epsilon = 0.0001; // TODO
 		if (m_nodes[i]->getType() == Node::ConnectionType::kJoin) {
 			//m_nodes[i]->m_position.setR(m_nodes[i]->m_position.r() + epsilon * m_forces[i].r);
 		}
diff --git a/cartocrow/flow_map/smooth_tree.h b/cartocrow/flow_map/smooth_tree.h
index 7472d6850..4ee4023c2 100644
--- a/cartocrow/flow_map/smooth_tree.h
+++ b/cartocrow/flow_map/smooth_tree.h
@@ -116,8 +116,8 @@ class SmoothTree {
 	///     \text{.}
 	/// \f]
 	Number<Inexact> computeAngleRestrictionCost(int i, int iChild1, int iChild2);
-	/// Applies angle restriction forces in \ref m_forces to the subdivision
-	/// node `i` and its children `iChild1` and `iChild2`.
+	/// Applies angle restriction forces in \ref m_forces to the join node `i`
+	/// and its children `iChild1` and `iChild2`.
 	///
 	/// The forces are defined by the partial derivatives of the angle
 	/// restriction cost (see \ref computeAngleRestrictionCost), which are:
@@ -151,8 +151,8 @@ class SmoothTree {
 	///     \text{.}
 	/// \f]
 	Number<Inexact> computeBalancingCost(int i, int iChild1, int iChild2);
-	/// Applies balancing forces in \ref m_forces to the subdivision node `i`
-	/// and its children `iChild1` and `iChild2`.
+	/// Applies balancing forces in \ref m_forces to the join node `i` and its
+	/// children `iChild1` and `iChild2`.
 	///
 	/// The forces are defined by the partial derivatives of the balancing cost
 	/// (see \ref computeBalancingCost), which are:
@@ -169,6 +169,43 @@ class SmoothTree {
 	/// \f}
 	/// et cetera.
 	void applyBalancingForce(int i, int iChild1, int iChild2);
+	/// Computes the straightening cost for the join node `i` at \f$(r,
+	/// \phi)\f$, with parent `iParent` at \f$(r_p, \phi_p)\f$ and children
+	/// `children` at \f$(r_{c_i}, \phi_{c_i})\f$.
+	///
+	/// \f[
+	///     F_\text{balance}(r, \phi, r_p, \phi_p, \{r_{c_i}\}, \{\phi_{c_i}\}) =
+	///     c_\text{straighten} \cdot
+	///     \Big( \alpha(r_p, \phi_p, r, \phi) -
+	///         \frac{\sum_{i \mid t_i \geq ct^*} t_i \alpha(r, \phi, r_{c_i}, \phi_{c_i})}
+	///             {\sum_{i \mid t_i \geq ct^*} t_i}
+	///     \Big)^2
+	///     \text{.}
+	/// \f]
+	Number<Inexact> computeStraighteningCost(int i, int iParent,
+	                                         const std::vector<std::shared_ptr<Node>>& children);
+	/// Applies straightening forces in \ref m_forces to the join node `i` at
+	/// \f$(r, \phi)\f$, with parent `iParent` at \f$(r_p, \phi_p)\f$ and
+	/// children `children` at \f$(r_{c_i}, \phi_{c_i})\f$.
+	///
+	/// The forces are defined by the partial derivatives of the straightening cost
+	/// (see \ref computeStraighteningCost), which are:
+	///
+	/// \f{multline}{
+	///     \frac{\partial F_\text{straighten}}{\partial r}(r, \phi, r_p, \phi_p, \{r_{c_i}\}, \{\phi_{c_i}\}) =
+	///     2c_\text{straighten} \cdot
+	///     \bigg( \alpha(r_p, \phi_p, r, \phi) -
+	///         \frac{\sum_{i \mid t_i \geq ct^*} t_i \alpha(r, \phi, r_{c_i}, \phi_{c_i})}
+	///             {\sum_{i \mid t_i \geq ct^*} t_i} \bigg) \\
+	///     \cdot \bigg(
+	///         \frac{\partial\alpha}{\partial r_2}(r_p, \phi_p, r, \phi) -
+	///         \frac{\sum_{i \mid t_i \geq ct^*} t_i \frac{\partial\alpha}{\partial r_1}(r, \phi, r_{c_i}, \phi_{c_i})}
+	///             {\sum_{i \mid t_i \geq ct^*} t_i} \bigg)
+	///     \text{;} \\
+	/// \f}
+	/// et cetera.
+	void applyStraighteningForce(int i, int iParent,
+	                             const std::vector<std::shared_ptr<Node>>& children);
 
 	// [ws] TODO Somewhere in the optimization code there seems to be a sign
 	// issue, although the result seems correct. We are computing forces by
@@ -181,6 +218,8 @@ class SmoothTree {
 	Number<Inexact> m_smoothing_factor = 0.4;
 	Number<Inexact> m_straightening_factor = 0.4;
 	Number<Inexact> m_angle_restriction_factor = 0.077;
+
+	Number<Inexact> m_relevantFlowFactor = 0.5;
 };
 
 } // namespace cartocrow::flow_map
diff --git a/cartocrow/flow_map/spiral.cpp b/cartocrow/flow_map/spiral.cpp
index a9c4a6ef9..aedb93013 100644
--- a/cartocrow/flow_map/spiral.cpp
+++ b/cartocrow/flow_map/spiral.cpp
@@ -70,7 +70,9 @@ Number<Inexact> Spiral::alpha(const PolarPoint& p1, const PolarPoint& p2) {
 Number<Inexact> Spiral::dAlphaDPhi1(const PolarPoint& p1, const PolarPoint& p2) {
 	const PolarPoint& source = p1.r() < p2.r() ? p2 : p1;
 	const PolarPoint& target = p1.r() < p2.r() ? p1 : p2;
-
+	if (target.r() == 0) {
+		return 0;
+	}
 	Number<Inexact> phiDiff = wrapAngle(target.phi() - source.phi(), -M_PI);
 	Number<Inexact> rDiffLog = std::log(target.r() / source.r());
 	return rDiffLog / (rDiffLog * rDiffLog + phiDiff * phiDiff);
@@ -83,7 +85,9 @@ Number<Inexact> Spiral::dAlphaDPhi2(const PolarPoint& p1, const PolarPoint& p2)
 Number<Inexact> Spiral::dAlphaDR1(const PolarPoint& p1, const PolarPoint& p2) {
 	const PolarPoint& source = p1.r() < p2.r() ? p2 : p1;
 	const PolarPoint& target = p1.r() < p2.r() ? p1 : p2;
-
+	if (target.r() == 0) {
+		return 0;
+	}
 	Number<Inexact> phiDiff = wrapAngle(target.phi() - source.phi(), -M_PI);
 	Number<Inexact> rDiffLog = std::log(target.r() / source.r());
 	return (-phiDiff / p1.r()) / (rDiffLog * rDiffLog + phiDiff * phiDiff);
@@ -92,7 +96,9 @@ Number<Inexact> Spiral::dAlphaDR1(const PolarPoint& p1, const PolarPoint& p2) {
 Number<Inexact> Spiral::dAlphaDR2(const PolarPoint& p1, const PolarPoint& p2) {
 	const PolarPoint& source = p1.r() < p2.r() ? p2 : p1;
 	const PolarPoint& target = p1.r() < p2.r() ? p1 : p2;
-
+	if (source.r() == 0) {
+		return 0;
+	}
 	Number<Inexact> phiDiff = wrapAngle(target.phi() - source.phi(), -M_PI);
 	Number<Inexact> rDiffLog = std::log(target.r() / source.r());
 	return (phiDiff / p2.r()) / (rDiffLog * rDiffLog + phiDiff * phiDiff);

From feb5a535c255f41df72cddc8eb401bffedd574dd Mon Sep 17 00:00:00 2001
From: Willem Sonke <willem@wimiso.nl>
Date: Tue, 12 Sep 2023 18:24:36 +0200
Subject: [PATCH 20/27] Show a cost graph in the optimization demo

---
 cartocrow/flow_map/smooth_tree.cpp            |  2 +-
 .../flow_map/optimization_demo/CMakeLists.txt |  1 +
 .../flow_map/optimization_demo/cost_graph.cpp | 79 +++++++++++++++++++
 demos/flow_map/optimization_demo/cost_graph.h | 64 +++++++++++++++
 .../optimization_demo/optimization_demo.cpp   | 23 +++++-
 .../optimization_demo/optimization_demo.h     |  4 +
 6 files changed, 168 insertions(+), 5 deletions(-)
 create mode 100644 demos/flow_map/optimization_demo/cost_graph.cpp
 create mode 100644 demos/flow_map/optimization_demo/cost_graph.h

diff --git a/cartocrow/flow_map/smooth_tree.cpp b/cartocrow/flow_map/smooth_tree.cpp
index cb4c218bf..c65587aaa 100644
--- a/cartocrow/flow_map/smooth_tree.cpp
+++ b/cartocrow/flow_map/smooth_tree.cpp
@@ -284,7 +284,7 @@ void SmoothTree::optimize() {
 	for (int i = 0; i < m_nodes.size(); i++) {
 		Number<Inexact> epsilon = 0.0001; // TODO
 		if (m_nodes[i]->getType() == Node::ConnectionType::kJoin) {
-			//m_nodes[i]->m_position.setR(m_nodes[i]->m_position.r() + epsilon * m_forces[i].r);
+			m_nodes[i]->m_position.setR(m_nodes[i]->m_position.r() + epsilon * m_forces[i].r);
 		}
 		if (m_nodes[i]->getType() == Node::ConnectionType::kJoin ||
 		    m_nodes[i]->getType() == Node::ConnectionType::kSubdivision) {
diff --git a/demos/flow_map/optimization_demo/CMakeLists.txt b/demos/flow_map/optimization_demo/CMakeLists.txt
index 1c9b3ba80..822d233cf 100644
--- a/demos/flow_map/optimization_demo/CMakeLists.txt
+++ b/demos/flow_map/optimization_demo/CMakeLists.txt
@@ -1,5 +1,6 @@
 set(SOURCES
     optimization_demo.cpp
+    cost_graph.cpp
 )
 
 add_executable(optimization_demo ${SOURCES})
diff --git a/demos/flow_map/optimization_demo/cost_graph.cpp b/demos/flow_map/optimization_demo/cost_graph.cpp
new file mode 100644
index 000000000..ca60b31d4
--- /dev/null
+++ b/demos/flow_map/optimization_demo/cost_graph.cpp
@@ -0,0 +1,79 @@
+/*
+The Necklace Map console application implements the algorithmic
+geo-visualization method by the same name, developed by
+Bettina Speckmann and Kevin Verbeek at TU Eindhoven
+(DOI: 10.1109/TVCG.2010.180 & 10.1142/S021819591550003X).
+Copyright (C) 2021  Netherlands eScience Center and TU Eindhoven
+
+This program is free software: you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation, either version 3 of the License, or
+(at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with this program.  If not, see <https://www.gnu.org/licenses/>.
+*/
+
+#include "cost_graph.h"
+
+#include <QPainter>
+#include <QPainterPath>
+
+using namespace cartocrow;
+
+CostGraph::CostGraph() {
+	setMinimumSize(QSize(300, 200));
+}
+
+void CostGraph::addStep(DataPoint costs) {
+	m_dataPoints.push_back(costs);
+	m_maxCost = std::max(m_maxCost, costs.m_obstacle_cost + costs.m_smoothing_cost +
+	                                    costs.m_angle_restriction_cost + costs.m_balancing_cost +
+	                                    costs.m_straightening_cost);
+	update();
+}
+
+void CostGraph::clear() {
+	m_dataPoints.clear();
+	m_maxCost = 0;
+	update();
+}
+
+void CostGraph::paintEvent(QPaintEvent* event) {
+	QPainter painter(this);
+	if (m_dataPoints.empty()) {
+		painter.drawText(QPoint(20, 20), "No data points");
+		return;
+	}
+
+	painter.setRenderHint(QPainter::Antialiasing);
+
+	int marginSize = 20;
+	int graphWidth = width() - 2 * marginSize;
+	int graphHeight = height() - 2 * marginSize;
+
+	painter.translate(QPoint(marginSize, graphHeight + marginSize));
+	painter.scale(1, -1);
+	painter.drawLine(0, 0, 0, graphHeight);
+	painter.drawLine(0, 0, graphWidth, 0);
+
+	QPainterPath dataLine;
+		dataLine.moveTo(0, graphHeight * m_dataPoints[0].m_obstacle_cost / m_maxCost);
+	for (int i = 1; i < m_dataPoints.size(); ++i) {
+		dataLine.lineTo(graphWidth * i / std::max(100.0f, static_cast<float>(m_dataPoints.size())),
+		                graphHeight * m_dataPoints[i].m_obstacle_cost / m_maxCost);
+	}
+	painter.setPen(QPen(QColor(240, 90, 40), 2));
+	painter.drawPath(dataLine);
+	painter.setBrush(QColor(240, 90, 40));
+	painter.drawEllipse(
+	    QPointF(graphWidth * (m_dataPoints.size() - 1) /
+	               std::max(100.0f, static_cast<float>(m_dataPoints.size())),
+	           graphHeight * m_dataPoints[m_dataPoints.size() - 1].m_obstacle_cost / m_maxCost),
+	    2, 2);
+}
diff --git a/demos/flow_map/optimization_demo/cost_graph.h b/demos/flow_map/optimization_demo/cost_graph.h
new file mode 100644
index 000000000..ee020c8c7
--- /dev/null
+++ b/demos/flow_map/optimization_demo/cost_graph.h
@@ -0,0 +1,64 @@
+/*
+The Necklace Map console application implements the algorithmic
+geo-visualization method by the same name, developed by
+Bettina Speckmann and Kevin Verbeek at TU Eindhoven
+(DOI: 10.1109/TVCG.2010.180 & 10.1142/S021819591550003X).
+Copyright (C) 2021  Netherlands eScience Center and TU Eindhoven
+
+This program is free software: you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation, either version 3 of the License, or
+(at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with this program.  If not, see <https://www.gnu.org/licenses/>.
+*/
+
+#ifndef CARTOCROW_DEMOS_FLOWMAP_OPTIMIZATIONDEMO_COSTGRAPH_H
+#define CARTOCROW_DEMOS_FLOWMAP_OPTIMIZATIONDEMO_COSTGRAPH_H
+
+#include <QWidget>
+
+#include <vector>
+
+#include "cartocrow/core/core.h"
+
+using namespace cartocrow;
+
+/// Simple widget that displays the cost of a \ref SmoothTree as it changes
+/// during the optimization procedure.
+class CostGraph : public QWidget {
+	Q_OBJECT;
+
+  public:
+	/// One data point containing all cost types.
+	struct DataPoint {
+		Number<Inexact> m_obstacle_cost;
+		Number<Inexact> m_smoothing_cost;
+		Number<Inexact> m_angle_restriction_cost;
+		Number<Inexact> m_balancing_cost;
+		Number<Inexact> m_straightening_cost;
+	};
+	/// Creates a cost graph without any cost data.
+	CostGraph();
+	/// Adds a new cost data point to the graph.
+	void addStep(DataPoint costs);
+	/// Clears all the cost data points.
+	void clear();
+
+  protected:
+	void paintEvent(QPaintEvent* event) override;
+
+  private:
+	/// The ordered list of cost data points.
+	std::vector<DataPoint> m_dataPoints;
+	/// The maximum total cost we've seen so far.
+	Number<Inexact> m_maxCost = 0;
+};
+
+#endif
diff --git a/demos/flow_map/optimization_demo/optimization_demo.cpp b/demos/flow_map/optimization_demo/optimization_demo.cpp
index 20468f0cb..56ce90587 100644
--- a/demos/flow_map/optimization_demo/optimization_demo.cpp
+++ b/demos/flow_map/optimization_demo/optimization_demo.cpp
@@ -23,6 +23,7 @@ along with this program.  If not, see <https://www.gnu.org/licenses/>.
 
 #include <QApplication>
 #include <QCheckBox>
+#include <QDockWidget>
 #include <QStatusBar>
 
 #include "cartocrow/core/core.h"
@@ -64,6 +65,12 @@ OptimizationDemo::OptimizationDemo() {
 	m_renderer->setGridMode(GeometryWidget::GridMode::POLAR);
 	setCentralWidget(m_renderer);
 
+	m_costGraph = new CostGraph();
+	QDockWidget* dockWidget = new QDockWidget();
+	dockWidget->setWindowTitle("Cost history");
+	dockWidget->setWidget(m_costGraph);
+	addDockWidget(Qt::RightDockWidgetArea, dockWidget);
+
 	QToolBar* toolBar = new QToolBar();
 	toolBar->addSeparator();
 	m_optimizeButton = new QPushButton("Run optimization");
@@ -78,13 +85,17 @@ OptimizationDemo::OptimizationDemo() {
 		}
 	});
 	connect(m_optimizeTimer, &QTimer::timeout, [&]() {
-		m_smoothTree->optimize();
+		for (int i = 0; i < 100; i++) {
+			m_iterationCount++;
+			m_smoothTree->optimize();
+			updateCostLabel();
+		}
 		m_renderer->update();
-		updateCostLabel();
 	});
 	m_optimizeOneStepButton = new QPushButton("Optimize one step");
 	toolBar->addWidget(m_optimizeOneStepButton);
 	connect(m_optimizeOneStepButton, &QPushButton::clicked, [&]() {
+		m_iterationCount++;
 		m_smoothTree->optimize();
 		m_renderer->update();
 		updateCostLabel();
@@ -119,6 +130,8 @@ OptimizationDemo::OptimizationDemo() {
 
 void OptimizationDemo::recalculate() {
 	m_renderer->clear();
+	m_costGraph->clear();
+	m_iterationCount = 0;
 	auto tree = std::make_shared<SpiralTree>(Point<Inexact>(0, 0), m_alpha);
 	for (const auto& place : m_places) {
 		tree->addPlace("", *place, 1);
@@ -144,10 +157,12 @@ void OptimizationDemo::recalculate() {
 }
 
 void OptimizationDemo::updateCostLabel() {
-	m_costLabel->setText(QString("Cost function: %1").arg(m_smoothTree->computeCost()));
+	m_costGraph->addStep({m_smoothTree->computeCost(), 0, 0, 0, 0});
+	m_costLabel->setText(QString("Iteration %1   |   Cost function: %2")
+	                         .arg(m_iterationCount)
+	                         .arg(m_smoothTree->computeCost()));
 }
 
-
 int main(int argc, char* argv[]) {
 	QApplication app(argc, argv);
 	OptimizationDemo demo;
diff --git a/demos/flow_map/optimization_demo/optimization_demo.h b/demos/flow_map/optimization_demo/optimization_demo.h
index 533c2c26b..9b1cd6f72 100644
--- a/demos/flow_map/optimization_demo/optimization_demo.h
+++ b/demos/flow_map/optimization_demo/optimization_demo.h
@@ -32,6 +32,8 @@ along with this program.  If not, see <https://www.gnu.org/licenses/>.
 #include "cartocrow/flow_map/smooth_tree.h"
 #include "cartocrow/renderer/geometry_widget.h"
 
+#include "cost_graph.h"
+
 using namespace cartocrow;
 using namespace cartocrow::renderer;
 
@@ -48,6 +50,7 @@ class OptimizationDemo : public QMainWindow {
 
 	std::vector<std::shared_ptr<Point<Inexact>>> m_places;
 	std::shared_ptr<cartocrow::flow_map::SmoothTree> m_smoothTree;
+	int m_iterationCount = 0;
 
 	GeometryWidget* m_renderer;
 	QSlider* m_alphaSlider;
@@ -56,4 +59,5 @@ class OptimizationDemo : public QMainWindow {
 	QPushButton* m_optimizeOneStepButton;
 	QTimer* m_optimizeTimer;
 	QLabel* m_costLabel;
+	CostGraph* m_costGraph;
 };

From bf4296d934f3d76d30d4ce3781c36d86c516f86f Mon Sep 17 00:00:00 2001
From: Willem Sonke <willem@wimiso.nl>
Date: Wed, 20 Sep 2023 13:56:01 +0200
Subject: [PATCH 21/27] Fix sign issues in the smooth tree optimization
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

The problem was that the α function was inverted (due to confusion over
which of the two input points is the one further from the origin).

Also renamed "force" to "gradient" to reduce sign confusion.
---
 cartocrow/flow_map/smooth_tree.cpp            | 153 +++++++++---------
 cartocrow/flow_map/smooth_tree.h              |  21 +--
 cartocrow/flow_map/spiral.cpp                 |  43 ++---
 cartocrow/flow_map/spiral.h                   |  19 ++-
 .../flow_map/optimization_demo/cost_graph.cpp |   2 +-
 .../optimization_demo/optimization_demo.cpp   |   4 +-
 6 files changed, 120 insertions(+), 122 deletions(-)

diff --git a/cartocrow/flow_map/smooth_tree.cpp b/cartocrow/flow_map/smooth_tree.cpp
index c65587aaa..41ef1915f 100644
--- a/cartocrow/flow_map/smooth_tree.cpp
+++ b/cartocrow/flow_map/smooth_tree.cpp
@@ -41,7 +41,8 @@ SmoothTree::SmoothTree(const std::shared_ptr<SpiralTree> spiralTree)
 	}
 }
 
-std::shared_ptr<Node> SmoothTree::constructSmoothTree(const std::shared_ptr<Node>& node, Number<Inexact> maxRStep) {
+std::shared_ptr<Node> SmoothTree::constructSmoothTree(const std::shared_ptr<Node>& node,
+                                                      Number<Inexact> maxRStep) {
 	auto smoothNode = std::make_shared<Node>(node->m_position);
 	m_nodes.push_back(smoothNode);
 	if (node->getType() == Node::ConnectionType::kLeaf) {
@@ -58,8 +59,8 @@ std::shared_ptr<Node> SmoothTree::constructSmoothTree(const std::shared_ptr<Node
 			int subdivisionCount = std::ceil(rRange / maxRStep);
 			auto previous = smoothNode;
 			for (double i = 1; i < subdivisionCount; ++i) {
-				const PolarPoint position =
-					spiral.evaluate(spiral.parameterForR(rMin + i * (rMax - rMin) / subdivisionCount));
+				const PolarPoint position = spiral.evaluate(
+				    spiral.parameterForR(rMin + i * (rMax - rMin) / subdivisionCount));
 				auto subdivision = std::make_shared<Node>(position);
 				subdivision->m_flow = smoothChild->m_flow;
 				m_nodes.push_back(subdivision);
@@ -86,25 +87,25 @@ Number<Inexact> SmoothTree::computeSmoothingCost(int i, int iParent, int iChild)
 	return m_smoothing_factor * std::pow(Spiral::alpha(p, n) - Spiral::alpha(n, c), 2);
 }
 
-void SmoothTree::applySmoothingForce(int i, int iParent, int iChild) {
+void SmoothTree::applySmoothingGradient(int i, int iParent, int iChild) {
 	PolarPoint n = m_nodes[i]->m_position;
 	PolarPoint p = m_nodes[iParent]->m_position;
 	PolarPoint c = m_nodes[iChild]->m_position;
 
-	m_forces[i].r += 2 * m_smoothing_factor * (Spiral::alpha(p, n) - Spiral::alpha(n, c)) *
-	                 (Spiral::dAlphaDR2(p, n) - Spiral::dAlphaDR1(n, c));
-	m_forces[i].phi += 2 * m_smoothing_factor * (Spiral::alpha(p, n) - Spiral::alpha(n, c)) *
-	                   (Spiral::dAlphaDPhi2(p, n) - Spiral::dAlphaDPhi1(n, c));
+	m_gradient[i].r += 2 * m_smoothing_factor * (Spiral::alpha(p, n) - Spiral::alpha(n, c)) *
+	                   (Spiral::dAlphaDR2(p, n) - Spiral::dAlphaDR1(n, c));
+	m_gradient[i].phi += 2 * m_smoothing_factor * (Spiral::alpha(p, n) - Spiral::alpha(n, c)) *
+	                     (Spiral::dAlphaDPhi2(p, n) - Spiral::dAlphaDPhi1(n, c));
 
-	m_forces[iParent].r += 2 * m_smoothing_factor * (Spiral::alpha(p, n) - Spiral::alpha(n, c)) *
-	                       Spiral::dAlphaDR1(p, n);
-	m_forces[iParent].phi += 2 * m_smoothing_factor * (Spiral::alpha(p, n) - Spiral::alpha(n, c)) *
-	                         Spiral::dAlphaDPhi1(p, n);
+	m_gradient[iParent].r += 2 * m_smoothing_factor * (Spiral::alpha(p, n) - Spiral::alpha(n, c)) *
+	                         Spiral::dAlphaDR1(p, n);
+	m_gradient[iParent].phi += 2 * m_smoothing_factor * (Spiral::alpha(p, n) - Spiral::alpha(n, c)) *
+	                           Spiral::dAlphaDPhi1(p, n);
 
-	m_forces[iChild].r += 2 * m_smoothing_factor * (Spiral::alpha(p, n) - Spiral::alpha(n, c)) *
-	                      -Spiral::dAlphaDR2(n, c);
-	m_forces[iChild].phi += 2 * m_smoothing_factor * (Spiral::alpha(p, n) - Spiral::alpha(n, c)) *
-	                        -Spiral::dAlphaDPhi2(n, c);
+	m_gradient[iChild].r += 2 * m_smoothing_factor * (Spiral::alpha(p, n) - Spiral::alpha(n, c)) *
+	                        -Spiral::dAlphaDR2(n, c);
+	m_gradient[iChild].phi += 2 * m_smoothing_factor * (Spiral::alpha(p, n) - Spiral::alpha(n, c)) *
+	                          -Spiral::dAlphaDPhi2(n, c);
 }
 
 Number<Inexact> SmoothTree::computeAngleRestrictionCost(int i, int iChild1, int iChild2) {
@@ -115,26 +116,26 @@ Number<Inexact> SmoothTree::computeAngleRestrictionCost(int i, int iChild1, int
 	                                     std::log(1.0 / std::cos(Spiral::alpha(n, c2))));
 }
 
-void SmoothTree::applyAngleRestrictionForce(int i, int iChild1, int iChild2) {
+void SmoothTree::applyAngleRestrictionGradient(int i, int iChild1, int iChild2) {
 	PolarPoint n = m_nodes[i]->m_position;
 	PolarPoint c1 = m_nodes[iChild1]->m_position;
 	PolarPoint c2 = m_nodes[iChild2]->m_position;
 
-	m_forces[i].r +=
+	m_gradient[i].r +=
 	    m_angle_restriction_factor * (Spiral::dAlphaDR1(n, c1) * std::tan(Spiral::alpha(n, c1)) +
 	                                  Spiral::dAlphaDR1(n, c2) * std::tan(Spiral::alpha(n, c2)));
-	m_forces[i].phi +=
+	m_gradient[i].phi +=
 	    m_angle_restriction_factor * (Spiral::dAlphaDPhi1(n, c1) * std::tan(Spiral::alpha(n, c1)) +
 	                                  Spiral::dAlphaDPhi1(n, c2) * std::tan(Spiral::alpha(n, c2)));
 
-	m_forces[iChild1].r +=
+	m_gradient[iChild1].r +=
 	    m_angle_restriction_factor * Spiral::dAlphaDR2(n, c1) * std::tan(Spiral::alpha(n, c1));
-	m_forces[iChild1].phi +=
+	m_gradient[iChild1].phi +=
 	    m_angle_restriction_factor * Spiral::dAlphaDPhi2(n, c1) * std::tan(Spiral::alpha(n, c1));
 
-	m_forces[iChild2].r +=
+	m_gradient[iChild2].r +=
 	    m_angle_restriction_factor * Spiral::dAlphaDR2(n, c2) * std::tan(Spiral::alpha(n, c2));
-	m_forces[iChild2].phi +=
+	m_gradient[iChild2].phi +=
 	    m_angle_restriction_factor * Spiral::dAlphaDPhi2(n, c2) * std::tan(Spiral::alpha(n, c2));
 }
 
@@ -146,31 +147,33 @@ Number<Inexact> SmoothTree::computeBalancingCost(int i, int iChild1, int iChild2
 	       std::log(1 / std::sin(0.5 * (Spiral::alpha(n, c1) - Spiral::alpha(n, c2))));
 }
 
-void SmoothTree::applyBalancingForce(int i, int iChild1, int iChild2) {
+void SmoothTree::applyBalancingGradient(int i, int iChild1, int iChild2) {
 	PolarPoint n = m_nodes[i]->m_position;
 	PolarPoint c1 = m_nodes[iChild1]->m_position;
 	PolarPoint c2 = m_nodes[iChild2]->m_position;
 
-	m_forces[i].r += m_angle_restriction_factor * -std::pow(std::tan(m_restrictingAngle), 2) *
-	                 (1 / std::tan(0.5 * (Spiral::alpha(n, c1) - Spiral::alpha(n, c2)))) *
-	                 (Spiral::dAlphaDR1(n, c1) - Spiral::dAlphaDR1(n, c2));
-	m_forces[i].phi += m_angle_restriction_factor * -std::pow(std::tan(m_restrictingAngle), 2) *
-	                 (1 / std::tan(0.5 * (Spiral::alpha(n, c1) - Spiral::alpha(n, c2)))) *
-	                 (Spiral::dAlphaDPhi1(n, c1) - Spiral::dAlphaDPhi1(n, c2));
-
-	m_forces[iChild1].r += m_angle_restriction_factor * -std::pow(std::tan(m_restrictingAngle), 2) *
-	                       (1 / std::tan(0.5 * (Spiral::alpha(n, c1) - Spiral::alpha(n, c2)))) *
-	                       Spiral::dAlphaDR2(n, c1);
-	m_forces[iChild1].phi += m_angle_restriction_factor * -std::pow(std::tan(m_restrictingAngle), 2) *
+	m_gradient[i].r += m_angle_restriction_factor * -std::pow(std::tan(m_restrictingAngle), 2) *
+	                   (1 / std::tan(0.5 * (Spiral::alpha(n, c1) - Spiral::alpha(n, c2)))) *
+	                   (Spiral::dAlphaDR1(n, c1) - Spiral::dAlphaDR1(n, c2));
+	m_gradient[i].phi += m_angle_restriction_factor * -std::pow(std::tan(m_restrictingAngle), 2) *
+	                     (1 / std::tan(0.5 * (Spiral::alpha(n, c1) - Spiral::alpha(n, c2)))) *
+	                     (Spiral::dAlphaDPhi1(n, c1) - Spiral::dAlphaDPhi1(n, c2));
+
+	m_gradient[iChild1].r += m_angle_restriction_factor * -std::pow(std::tan(m_restrictingAngle), 2) *
 	                         (1 / std::tan(0.5 * (Spiral::alpha(n, c1) - Spiral::alpha(n, c2)))) *
-	                         Spiral::dAlphaDPhi2(n, c1);
+	                         Spiral::dAlphaDR2(n, c1);
+	m_gradient[iChild1].phi += m_angle_restriction_factor *
+	                           -std::pow(std::tan(m_restrictingAngle), 2) *
+	                           (1 / std::tan(0.5 * (Spiral::alpha(n, c1) - Spiral::alpha(n, c2)))) *
+	                           Spiral::dAlphaDPhi2(n, c1);
 
-	m_forces[iChild2].r += m_angle_restriction_factor * -std::pow(std::tan(m_restrictingAngle), 2) *
-	                       (1 / std::tan(0.5 * (Spiral::alpha(n, c1) - Spiral::alpha(n, c2)))) *
-	                       -Spiral::dAlphaDR2(n, c2);
-	m_forces[iChild2].phi += m_angle_restriction_factor * -std::pow(std::tan(m_restrictingAngle), 2) *
+	m_gradient[iChild2].r += m_angle_restriction_factor * -std::pow(std::tan(m_restrictingAngle), 2) *
 	                         (1 / std::tan(0.5 * (Spiral::alpha(n, c1) - Spiral::alpha(n, c2)))) *
-	                         -Spiral::dAlphaDPhi2(n, c2);
+	                         -Spiral::dAlphaDR2(n, c2);
+	m_gradient[iChild2].phi += m_angle_restriction_factor *
+	                           -std::pow(std::tan(m_restrictingAngle), 2) *
+	                           (1 / std::tan(0.5 * (Spiral::alpha(n, c1) - Spiral::alpha(n, c2)))) *
+	                           -Spiral::dAlphaDPhi2(n, c2);
 }
 
 Number<Inexact>
@@ -196,8 +199,8 @@ SmoothTree::computeStraighteningCost(int i, int iParent,
 	return m_straightening_factor * std::pow(Spiral::alpha(p, n) - numerator / denominator, 2);
 }
 
-void SmoothTree::applyStraighteningForce(int i, int iParent,
-                                         const std::vector<std::shared_ptr<Node>>& children) {
+void SmoothTree::applyStraighteningGradient(int i, int iParent,
+                                            const std::vector<std::shared_ptr<Node>>& children) {
 	PolarPoint n = m_nodes[i]->m_position;
 	PolarPoint p = m_nodes[iParent]->m_position;
 	Number<Inexact> maxFlow = 0;
@@ -219,25 +222,28 @@ void SmoothTree::applyStraighteningForce(int i, int iParent,
 			denominator += child->m_flow;
 		}
 	}
-	m_forces[i].r += 2 * m_straightening_factor * (Spiral::alpha(p, n) - numerator / denominator) *
-	                 (Spiral::dAlphaDR2(p, n) - numeratorDR1 / denominator);
-	m_forces[i].phi += 2 * m_straightening_factor * (Spiral::alpha(p, n) - numerator / denominator) *
-	                   (Spiral::dAlphaDPhi2(p, n) - numeratorDPhi1 / denominator);
-
-	m_forces[iParent].r += 2 * m_straightening_factor * (Spiral::alpha(p, n) - numerator / denominator) *
-	                 Spiral::dAlphaDR1(p, n);
-	m_forces[iParent].phi += 2 * m_straightening_factor * (Spiral::alpha(p, n) - numerator / denominator) *
-	                   Spiral::dAlphaDPhi1(p, n);
+	m_gradient[i].r += 2 * m_straightening_factor * (Spiral::alpha(p, n) - numerator / denominator) *
+	                   (Spiral::dAlphaDR2(p, n) - numeratorDR1 / denominator);
+	m_gradient[i].phi += 2 * m_straightening_factor *
+	                     (Spiral::alpha(p, n) - numerator / denominator) *
+	                     (Spiral::dAlphaDPhi2(p, n) - numeratorDPhi1 / denominator);
+
+	m_gradient[iParent].r += 2 * m_straightening_factor *
+	                         (Spiral::alpha(p, n) - numerator / denominator) *
+	                         Spiral::dAlphaDR1(p, n);
+	m_gradient[iParent].phi += 2 * m_straightening_factor *
+	                           (Spiral::alpha(p, n) - numerator / denominator) *
+	                           Spiral::dAlphaDPhi1(p, n);
 
 	for (const auto& child : children) {
 		if (child->m_flow > maxFlow) {
 			PolarPoint c = child->m_position;
-			m_forces[child->m_id].r += 2 * m_straightening_factor *
-			                           (Spiral::alpha(p, n) - numerator / denominator) *
-			                           -child->m_flow * Spiral::dAlphaDR2(n, c) / denominator;
-			m_forces[child->m_id].phi += 2 * m_straightening_factor *
-			                           (Spiral::alpha(p, n) - numerator / denominator) *
-			                           -child->m_flow * Spiral::dAlphaDPhi2(n, c) / denominator;
+			m_gradient[child->m_id].r += 2 * m_straightening_factor *
+			                             (Spiral::alpha(p, n) - numerator / denominator) *
+			                             -child->m_flow * Spiral::dAlphaDR2(n, c) / denominator;
+			m_gradient[child->m_id].phi += 2 * m_straightening_factor *
+			                               (Spiral::alpha(p, n) - numerator / denominator) *
+			                               -child->m_flow * Spiral::dAlphaDPhi2(n, c) / denominator;
 		}
 	}
 }
@@ -247,48 +253,47 @@ Number<Inexact> SmoothTree::computeCost() {
 	for (int i = 0; i < m_nodes.size(); i++) {
 		const auto& node = m_nodes[i];
 		if (node->getType() == Node::ConnectionType::kSubdivision) {
-			cost += computeSmoothingCost(i, m_nodes[i]->m_parent->m_id, m_nodes[i]->m_children[0]->m_id);
+			cost +=
+			    computeSmoothingCost(i, m_nodes[i]->m_parent->m_id, m_nodes[i]->m_children[0]->m_id);
 		} else if (node->getType() == Node::ConnectionType::kJoin) {
-			cost += computeAngleRestrictionCost(
+			/*cost += computeAngleRestrictionCost(
 			    i, m_nodes[i]->m_children[0]->m_id,
 			    m_nodes[i]->m_children[m_nodes[i]->m_children.size() - 1]->m_id);
-			cost += computeBalancingCost(
-			    i, m_nodes[i]->m_children[0]->m_id,
-			    m_nodes[i]->m_children[m_nodes[i]->m_children.size() - 1]->m_id);
-			cost += computeStraighteningCost(
-			    i, m_nodes[i]->m_parent->m_id,
-			    m_nodes[i]->m_children);
+			cost +=
+			    computeBalancingCost(i, m_nodes[i]->m_children[0]->m_id,
+			                         m_nodes[i]->m_children[m_nodes[i]->m_children.size() - 1]->m_id);
+			cost += computeStraighteningCost(i, m_nodes[i]->m_parent->m_id, m_nodes[i]->m_children);*/
 		}
 	}
 	return cost;
 }
 
 void SmoothTree::optimize() {
-	m_forces = std::vector<PolarForce>(m_nodes.size(), PolarForce{});
+	m_gradient = std::vector<PolarGradient>(m_nodes.size(), PolarGradient{});
 	for (int i = 0; i < m_nodes.size(); i++) {
 		const auto& node = m_nodes[i];
 		if (node->getType() == Node::ConnectionType::kSubdivision) {
-			applySmoothingForce(i, m_nodes[i]->m_parent->m_id, m_nodes[i]->m_children[0]->m_id);
+			applySmoothingGradient(i, m_nodes[i]->m_parent->m_id, m_nodes[i]->m_children[0]->m_id);
 		} else if (node->getType() == Node::ConnectionType::kJoin) {
-			applyAngleRestrictionForce(
+			/*applyAngleRestrictionGradient(
 			    i, m_nodes[i]->m_children[0]->m_id,
 			    m_nodes[i]->m_children[m_nodes[i]->m_children.size() - 1]->m_id);
-			applyBalancingForce(
+			applyBalancingGradient(
 			    i, m_nodes[i]->m_children[0]->m_id,
 			    m_nodes[i]->m_children[m_nodes[i]->m_children.size() - 1]->m_id);
-			applyStraighteningForce(
+			applyStraighteningGradient(
 			    i, m_nodes[i]->m_parent->m_id,
-			    m_nodes[i]->m_children);
+			    m_nodes[i]->m_children);*/
 		}
 	}
 	for (int i = 0; i < m_nodes.size(); i++) {
 		Number<Inexact> epsilon = 0.0001; // TODO
 		if (m_nodes[i]->getType() == Node::ConnectionType::kJoin) {
-			m_nodes[i]->m_position.setR(m_nodes[i]->m_position.r() + epsilon * m_forces[i].r);
+			m_nodes[i]->m_position.setR(m_nodes[i]->m_position.r() - epsilon * m_gradient[i].r);
 		}
 		if (m_nodes[i]->getType() == Node::ConnectionType::kJoin ||
 		    m_nodes[i]->getType() == Node::ConnectionType::kSubdivision) {
-			m_nodes[i]->m_position.setPhi(m_nodes[i]->m_position.phi() + epsilon * m_forces[i].phi);
+			m_nodes[i]->m_position.setPhi(m_nodes[i]->m_position.phi() - epsilon * m_gradient[i].phi);
 		}
 	}
 }
diff --git a/cartocrow/flow_map/smooth_tree.h b/cartocrow/flow_map/smooth_tree.h
index 4ee4023c2..9490a5edb 100644
--- a/cartocrow/flow_map/smooth_tree.h
+++ b/cartocrow/flow_map/smooth_tree.h
@@ -59,11 +59,11 @@ class SmoothTree {
 
 	std::shared_ptr<Node> constructSmoothTree(const std::shared_ptr<Node>& node, Number<Inexact> maxRStep);
 
-	struct PolarForce {
+	struct PolarGradient {
 		Number<Inexact> r = 0;
 		Number<Inexact> phi = 0;
 	};
-	std::vector<PolarForce> m_forces;
+	std::vector<PolarGradient> m_gradient;
 
 	/// Computes the smoothing cost for the subdivision node `i` at \f$(r,
 	/// \phi)\f$, with parent `iParent` at \f$(r_p, \phi_p)\f$ and child
@@ -102,7 +102,7 @@ class SmoothTree {
 	///     \text{;} \\
 	/// \f}
 	/// et cetera.
-	void applySmoothingForce(int i, int iParent, int iChild);
+	void applySmoothingGradient(int i, int iParent, int iChild);
 	
 	/// Computes the angle restriction cost for the join node `i` at \f$(r,
 	/// \phi)\f$, with children `iChild1` at \f$(r_{c_1}, \phi_{c_1})\f$ and
@@ -137,7 +137,7 @@ class SmoothTree {
 	///     \text{;} \\
 	/// \f}
 	/// et cetera.
-	void applyAngleRestrictionForce(int i, int iChild1, int iChild2);
+	void applyAngleRestrictionGradient(int i, int iChild1, int iChild2);
 	/// Computes the balancing cost for the join node `i` at \f$(r, \phi)\f$,
 	/// with children `iChild1` at \f$(r_{c_1}, \phi_{c_1})\f$ and `iChild2` at
 	/// \f$(r_{c_2}, \phi_{c_2})\f$.
@@ -168,7 +168,7 @@ class SmoothTree {
 	///     \text{;} \\
 	/// \f}
 	/// et cetera.
-	void applyBalancingForce(int i, int iChild1, int iChild2);
+	void applyBalancingGradient(int i, int iChild1, int iChild2);
 	/// Computes the straightening cost for the join node `i` at \f$(r,
 	/// \phi)\f$, with parent `iParent` at \f$(r_p, \phi_p)\f$ and children
 	/// `children` at \f$(r_{c_i}, \phi_{c_i})\f$.
@@ -204,15 +204,8 @@ class SmoothTree {
 	///     \text{;} \\
 	/// \f}
 	/// et cetera.
-	void applyStraighteningForce(int i, int iParent,
-	                             const std::vector<std::shared_ptr<Node>>& children);
-
-	// [ws] TODO Somewhere in the optimization code there seems to be a sign
-	// issue, although the result seems correct. We are computing forces by
-	// computing the gradient, but the forces are supposed to point towards the
-	// negative gradient, not the positive. Yet, we are never inverting the
-	// sign. I'm confused as to why this seems to be working fine; maybe the
-	// sign of the α function is incorrect?
+	void applyStraighteningGradient(int i, int iParent,
+	                                const std::vector<std::shared_ptr<Node>>& children);
 
 	Number<Inexact> m_obstacle_factor = 2.0;
 	Number<Inexact> m_smoothing_factor = 0.4;
diff --git a/cartocrow/flow_map/spiral.cpp b/cartocrow/flow_map/spiral.cpp
index aedb93013..607d5adb7 100644
--- a/cartocrow/flow_map/spiral.cpp
+++ b/cartocrow/flow_map/spiral.cpp
@@ -22,7 +22,6 @@ Created by tvl (t.vanlankveld@esciencecenter.nl) on 16-10-2020
 #include "spiral.h"
 
 #include <cmath>
-
 #include <stdexcept>
 
 namespace cartocrow::flow_map {
@@ -38,39 +37,23 @@ Spiral::Spiral(const PolarPoint& p1, const PolarPoint& p2)
     : m_anchor(p1.r() < p2.r() ? p2 : p1), m_angle(alpha(p1, p2)) {}
 
 Number<Inexact> Spiral::alpha(const PolarPoint& p1, const PolarPoint& p2) {
-	const PolarPoint& source = p1.r() < p2.r() ? p2 : p1;
-	const PolarPoint& target = p1.r() < p2.r() ? p1 : p2;
-
+	const PolarPoint& source = p1.r() < p2.r() ? p1 : p2;
+	const PolarPoint& target = p1.r() < p2.r() ? p2 : p1;
 	if (source.r() == target.r()) {
 		throw std::runtime_error(
 		    "Cannot compute α for a spiral connecting two points equidistant to the root");
 	}
-
-	// Computing a spiral through two points (p, q), with unknown angle (alpha):
-	//
-	// p = (R_p, phi_p)  ->  R(0) = R_p; phi(0) = phi_p
-	// q = (R_q, phi_q)  ->  R(t) = R_p * e^{-t}; phi(t) = phi_p + tan(alpha) * t
-	//
-	// R(0) = R_p * e^0 = R_p
-	// phi(0) = phi_p + tan(alpha) * 0 = phi_p
-	// Assuming R_q < R_p and 0 < t:
-	//    R(t) = R_q = R_p * e^{-t}  =>  e^{-t} = R_q / R_p  =>  t = -ln(R_q / R_p)
-	//    phi(t) = phi_q = phi_p + tan(alpha) * t  =>  tan(alpha) = (phi_q - phi_p) / t  =>  alpha = tan^-1((phi_q - phi_p) / t)
-	//
-	//    =>  alpha = tan^-1((phi_q - phi_p) / -ln(R_q / R_p))
-
-	if (target.r() == 0) {
+	if (source.r() == 0) {
 		return 0;
-	} else {
-		Number<Inexact> diff_phi = wrapAngle(target.phi() - source.phi(), -M_PI);
-		return std::atan(diff_phi / -std::log(target.r() / source.r()));
 	}
+	Number<Inexact> diff_phi = wrapAngle(target.phi() - source.phi(), -M_PI);
+	return std::atan(diff_phi / -std::log(target.r() / source.r()));
 }
 
 Number<Inexact> Spiral::dAlphaDPhi1(const PolarPoint& p1, const PolarPoint& p2) {
-	const PolarPoint& source = p1.r() < p2.r() ? p2 : p1;
-	const PolarPoint& target = p1.r() < p2.r() ? p1 : p2;
-	if (target.r() == 0) {
+	const PolarPoint& source = p1.r() < p2.r() ? p1 : p2;
+	const PolarPoint& target = p1.r() < p2.r() ? p2 : p1;
+	if (source.r() == 0) {
 		return 0;
 	}
 	Number<Inexact> phiDiff = wrapAngle(target.phi() - source.phi(), -M_PI);
@@ -83,9 +66,9 @@ Number<Inexact> Spiral::dAlphaDPhi2(const PolarPoint& p1, const PolarPoint& p2)
 }
 
 Number<Inexact> Spiral::dAlphaDR1(const PolarPoint& p1, const PolarPoint& p2) {
-	const PolarPoint& source = p1.r() < p2.r() ? p2 : p1;
-	const PolarPoint& target = p1.r() < p2.r() ? p1 : p2;
-	if (target.r() == 0) {
+	const PolarPoint& source = p1.r() < p2.r() ? p1 : p2;
+	const PolarPoint& target = p1.r() < p2.r() ? p2 : p1;
+	if (source.r() == 0) {
 		return 0;
 	}
 	Number<Inexact> phiDiff = wrapAngle(target.phi() - source.phi(), -M_PI);
@@ -94,8 +77,8 @@ Number<Inexact> Spiral::dAlphaDR1(const PolarPoint& p1, const PolarPoint& p2) {
 }
 
 Number<Inexact> Spiral::dAlphaDR2(const PolarPoint& p1, const PolarPoint& p2) {
-	const PolarPoint& source = p1.r() < p2.r() ? p2 : p1;
-	const PolarPoint& target = p1.r() < p2.r() ? p1 : p2;
+	const PolarPoint& source = p1.r() < p2.r() ? p1 : p2;
+	const PolarPoint& target = p1.r() < p2.r() ? p2 : p1;
 	if (source.r() == 0) {
 		return 0;
 	}
diff --git a/cartocrow/flow_map/spiral.h b/cartocrow/flow_map/spiral.h
index 1d7a91e8c..9a9d707b1 100644
--- a/cartocrow/flow_map/spiral.h
+++ b/cartocrow/flow_map/spiral.h
@@ -71,12 +71,29 @@ class Spiral {
 	/// Computes the \f$\alpha\f$ of the shortest logarithmic spiral connecting
 	/// the two given points \f$(r_1, \phi_1)\f$ and \f$(r_2, \phi_2)\f$.
 	///
-	/// This value can be computed as follows:
+	/// Assuming that \f$0 < r_1 < r_2\f$, this value can be computed as
+	/// follows:
 	/// \f[
 	///     \alpha (r_1, \phi_1, r_2, \phi_2) = \arctan \left(
 	///         \frac{\phi_2 - \phi_1}{-\ln(r_2 / r_1)}
 	///     \right) \text{.}
 	/// \f]
+	///
+	/// We can derive this as follows. Let the spiral be \f$p(t) = (r(t),
+	/// \phi(t))\f$, where \f$r(t) = r \cdot e^{-t}\f$ and \f$\phi(t) = \phi +
+	/// \tan(\alpha) \cdot t\f$. The spiral passes through the given points, say
+	/// at \f$t = 0\f$ for \f$(r_1, \phi_1)\f$ and at \f$t = t' >
+	/// 0\f$ for \f$(r_2, \phi_1)\f$. So we have:
+	/// \f{align}{
+	///     r(0)\ &=\ r \cdot e^{-0}\ =\ r\ =\ r_1; \\
+	///     r(t')\ &=\ r \cdot e^{-t'}\ =\ r_1 \cdot e^{-t'}
+	///         \ = r_2 \quad \Rightarrow \quad t'\ =\ -\ln(r_2 / r_1); \\
+    ///     \phi(0)\ &=\ \phi + \tan(\alpha) \cdot 0\ =\ \phi\ =\ \phi_1; \\
+    ///     \phi(t')\ &=\ \phi + \tan(\alpha) \cdot t'
+	///         \ =\ \phi_1 + \tan(\alpha) \cdot -\ln(r_2 / r_1)
+	///         \ =\ \phi_2 \quad \Rightarrow \quad \alpha
+	///         \ =\ \arctan \left( \frac{\phi_2 - \phi_1}{-\ln(r_2 / r_1)} \right). \\
+	/// \f}
 	static Number<Inexact> alpha(const PolarPoint& p1, const PolarPoint& p2);
 	/// Computes \f$\frac{\partial\alpha}{\partial r_1}\f$ for the \f$\alpha\f$
 	/// function (see \ref alpha).
diff --git a/demos/flow_map/optimization_demo/cost_graph.cpp b/demos/flow_map/optimization_demo/cost_graph.cpp
index ca60b31d4..23615f23f 100644
--- a/demos/flow_map/optimization_demo/cost_graph.cpp
+++ b/demos/flow_map/optimization_demo/cost_graph.cpp
@@ -63,7 +63,7 @@ void CostGraph::paintEvent(QPaintEvent* event) {
 	painter.drawLine(0, 0, graphWidth, 0);
 
 	QPainterPath dataLine;
-		dataLine.moveTo(0, graphHeight * m_dataPoints[0].m_obstacle_cost / m_maxCost);
+	dataLine.moveTo(0, graphHeight * m_dataPoints[0].m_obstacle_cost / m_maxCost);
 	for (int i = 1; i < m_dataPoints.size(); ++i) {
 		dataLine.lineTo(graphWidth * i / std::max(100.0f, static_cast<float>(m_dataPoints.size())),
 		                graphHeight * m_dataPoints[i].m_obstacle_cost / m_maxCost);
diff --git a/demos/flow_map/optimization_demo/optimization_demo.cpp b/demos/flow_map/optimization_demo/optimization_demo.cpp
index 56ce90587..dd0c9c7b5 100644
--- a/demos/flow_map/optimization_demo/optimization_demo.cpp
+++ b/demos/flow_map/optimization_demo/optimization_demo.cpp
@@ -141,13 +141,13 @@ void OptimizationDemo::recalculate() {
 	ReachableRegionAlgorithm::ReachableRegion reachableRegion = ReachableRegionAlgorithm(tree).run();
 	SpiralTreeObstructedAlgorithm(tree, reachableRegion).run();
 	m_smoothTree = std::make_shared<SmoothTree>(tree);
-	/*for (auto& node : m_smoothTree->m_nodes) {
+	for (auto& node : m_smoothTree->m_nodes) {
 		if (node->getType() == Node::ConnectionType::kSubdivision) {
 			node->m_position.setPhi(wrapAngle(
 			    node->m_position.phi() + static_cast<float>(std::rand()) / RAND_MAX * 0.2 - 0.1,
 			    -M_PI));
 		}
-	}*/
+	}
 	SmoothTreePainting::Options options;
 	auto painting = std::make_shared<SmoothTreePainting>(m_smoothTree, options);
 	m_renderer->addPainting(painting, "Smooth tree");

From c2b216b7c1c558660d26446f1e29a24d363123c1 Mon Sep 17 00:00:00 2001
From: Willem Sonke <willem@wimiso.nl>
Date: Wed, 20 Sep 2023 22:07:53 +0200
Subject: [PATCH 22/27] Distinguish different types of costs in the GUI

---
 cartocrow/flow_map/smooth_tree.cpp            | 68 ++++++++++++++-----
 cartocrow/flow_map/smooth_tree.h              |  8 +++
 .../flow_map/optimization_demo/cost_graph.cpp | 66 ++++++++++++++----
 demos/flow_map/optimization_demo/cost_graph.h | 15 ++++
 .../optimization_demo/optimization_demo.cpp   |  4 +-
 5 files changed, 129 insertions(+), 32 deletions(-)

diff --git a/cartocrow/flow_map/smooth_tree.cpp b/cartocrow/flow_map/smooth_tree.cpp
index 41ef1915f..092b26ccf 100644
--- a/cartocrow/flow_map/smooth_tree.cpp
+++ b/cartocrow/flow_map/smooth_tree.cpp
@@ -80,6 +80,18 @@ const std::vector<std::shared_ptr<Node>>& SmoothTree::nodes() const {
 	return m_nodes;
 }
 
+Number<Inexact> SmoothTree::computeSmoothingCost() {
+	Number<Inexact> cost = 0;
+	for (int i = 0; i < m_nodes.size(); i++) {
+		const auto& node = m_nodes[i];
+		if (node->getType() == Node::ConnectionType::kSubdivision) {
+			cost +=
+			    computeSmoothingCost(i, m_nodes[i]->m_parent->m_id, m_nodes[i]->m_children[0]->m_id);
+		}
+	}
+	return cost;
+}
+
 Number<Inexact> SmoothTree::computeSmoothingCost(int i, int iParent, int iChild) {
 	PolarPoint n = m_nodes[i]->m_position;
 	PolarPoint p = m_nodes[iParent]->m_position;
@@ -108,6 +120,19 @@ void SmoothTree::applySmoothingGradient(int i, int iParent, int iChild) {
 	                          -Spiral::dAlphaDPhi2(n, c);
 }
 
+Number<Inexact> SmoothTree::computeAngleRestrictionCost() {
+	Number<Inexact> cost = 0;
+	for (int i = 0; i < m_nodes.size(); i++) {
+		const auto& node = m_nodes[i];
+		if (node->getType() == Node::ConnectionType::kJoin) {
+			cost += computeAngleRestrictionCost(
+			    i, m_nodes[i]->m_children[0]->m_id,
+			    m_nodes[i]->m_children[m_nodes[i]->m_children.size() - 1]->m_id);
+		}
+	}
+	return cost;
+}
+
 Number<Inexact> SmoothTree::computeAngleRestrictionCost(int i, int iChild1, int iChild2) {
 	PolarPoint n = m_nodes[i]->m_position;
 	PolarPoint c1 = m_nodes[iChild1]->m_position;
@@ -139,6 +164,19 @@ void SmoothTree::applyAngleRestrictionGradient(int i, int iChild1, int iChild2)
 	    m_angle_restriction_factor * Spiral::dAlphaDPhi2(n, c2) * std::tan(Spiral::alpha(n, c2));
 }
 
+Number<Inexact> SmoothTree::computeBalancingCost() {
+	Number<Inexact> cost = 0;
+	for (int i = 0; i < m_nodes.size(); i++) {
+		const auto& node = m_nodes[i];
+		if (node->getType() == Node::ConnectionType::kJoin) {
+			cost +=
+			    computeBalancingCost(i, m_nodes[i]->m_children[0]->m_id,
+			                         m_nodes[i]->m_children[m_nodes[i]->m_children.size() - 1]->m_id);
+		}
+	}
+	return cost;
+}
+
 Number<Inexact> SmoothTree::computeBalancingCost(int i, int iChild1, int iChild2) {
 	PolarPoint n = m_nodes[i]->m_position;
 	PolarPoint c1 = m_nodes[iChild1]->m_position;
@@ -176,6 +214,17 @@ void SmoothTree::applyBalancingGradient(int i, int iChild1, int iChild2) {
 	                           -Spiral::dAlphaDPhi2(n, c2);
 }
 
+Number<Inexact> SmoothTree::computeStraighteningCost() {
+	Number<Inexact> cost = 0;
+	for (int i = 0; i < m_nodes.size(); i++) {
+		const auto& node = m_nodes[i];
+		if (node->getType() == Node::ConnectionType::kJoin) {
+			cost += computeStraighteningCost(i, m_nodes[i]->m_parent->m_id, m_nodes[i]->m_children);
+		}
+	}
+	return cost;
+}
+
 Number<Inexact>
 SmoothTree::computeStraighteningCost(int i, int iParent,
                                      const std::vector<std::shared_ptr<Node>>& children) {
@@ -249,23 +298,8 @@ void SmoothTree::applyStraighteningGradient(int i, int iParent,
 }
 
 Number<Inexact> SmoothTree::computeCost() {
-	Number<Inexact> cost = 0;
-	for (int i = 0; i < m_nodes.size(); i++) {
-		const auto& node = m_nodes[i];
-		if (node->getType() == Node::ConnectionType::kSubdivision) {
-			cost +=
-			    computeSmoothingCost(i, m_nodes[i]->m_parent->m_id, m_nodes[i]->m_children[0]->m_id);
-		} else if (node->getType() == Node::ConnectionType::kJoin) {
-			/*cost += computeAngleRestrictionCost(
-			    i, m_nodes[i]->m_children[0]->m_id,
-			    m_nodes[i]->m_children[m_nodes[i]->m_children.size() - 1]->m_id);
-			cost +=
-			    computeBalancingCost(i, m_nodes[i]->m_children[0]->m_id,
-			                         m_nodes[i]->m_children[m_nodes[i]->m_children.size() - 1]->m_id);
-			cost += computeStraighteningCost(i, m_nodes[i]->m_parent->m_id, m_nodes[i]->m_children);*/
-		}
-	}
-	return cost;
+	return computeSmoothingCost() + computeAngleRestrictionCost() + computeBalancingCost() +
+	       computeStraighteningCost();
 }
 
 void SmoothTree::optimize() {
diff --git a/cartocrow/flow_map/smooth_tree.h b/cartocrow/flow_map/smooth_tree.h
index 9490a5edb..b40e18a7b 100644
--- a/cartocrow/flow_map/smooth_tree.h
+++ b/cartocrow/flow_map/smooth_tree.h
@@ -44,6 +44,14 @@ class SmoothTree {
 
 	/// Computes the total cost of the tree.
 	Number<Inexact> computeCost();
+	/// Computes the smoothing cost of the entire tree.
+	Number<Inexact> computeSmoothingCost();
+	/// Computes the angle restriction cost of the entire tree.
+	Number<Inexact> computeAngleRestrictionCost();
+	/// Computes the balancing cost of the entire tree.
+	Number<Inexact> computeBalancingCost();
+	/// Computes the straightening cost of the entire tree.
+	Number<Inexact> computeStraighteningCost();
 
 	/// Performs one optimization step.
 	void optimize();
diff --git a/demos/flow_map/optimization_demo/cost_graph.cpp b/demos/flow_map/optimization_demo/cost_graph.cpp
index 23615f23f..7369d4ca2 100644
--- a/demos/flow_map/optimization_demo/cost_graph.cpp
+++ b/demos/flow_map/optimization_demo/cost_graph.cpp
@@ -23,6 +23,8 @@ along with this program.  If not, see <https://www.gnu.org/licenses/>.
 
 #include <QPainter>
 #include <QPainterPath>
+#include <QPalette>
+#include <QStyle>
 
 using namespace cartocrow;
 
@@ -58,22 +60,58 @@ void CostGraph::paintEvent(QPaintEvent* event) {
 	int graphHeight = height() - 2 * marginSize;
 
 	painter.translate(QPoint(marginSize, graphHeight + marginSize));
-	painter.scale(1, -1);
-	painter.drawLine(0, 0, 0, graphHeight);
+
+	// axes
+	painter.drawLine(0, 0, 0, -graphHeight);
 	painter.drawLine(0, 0, graphWidth, 0);
 
-	QPainterPath dataLine;
-	dataLine.moveTo(0, graphHeight * m_dataPoints[0].m_obstacle_cost / m_maxCost);
-	for (int i = 1; i < m_dataPoints.size(); ++i) {
-		dataLine.lineTo(graphWidth * i / std::max(100.0f, static_cast<float>(m_dataPoints.size())),
-		                graphHeight * m_dataPoints[i].m_obstacle_cost / m_maxCost);
-	}
+	// graphs
+	QPainterPath obstaclePath = createDataPath(1, graphWidth, graphHeight);
+	QPainterPath smoothingPath = createDataPath(2, graphWidth, graphHeight);
+	QPainterPath angleRestrictionPath = createDataPath(3, graphWidth, graphHeight);
+	QPainterPath balancingPath = createDataPath(4, graphWidth, graphHeight);
+	QPainterPath straighteningPath = createDataPath(5, graphWidth, graphHeight);
+	painter.drawPath(obstaclePath);
+	painter.drawPath(smoothingPath);
+	painter.drawPath(angleRestrictionPath);
+	painter.drawPath(balancingPath);
 	painter.setPen(QPen(QColor(240, 90, 40), 2));
-	painter.drawPath(dataLine);
+	painter.drawPath(straighteningPath);
 	painter.setBrush(QColor(240, 90, 40));
-	painter.drawEllipse(
-	    QPointF(graphWidth * (m_dataPoints.size() - 1) /
-	               std::max(100.0f, static_cast<float>(m_dataPoints.size())),
-	           graphHeight * m_dataPoints[m_dataPoints.size() - 1].m_obstacle_cost / m_maxCost),
-	    2, 2);
+
+	// dot
+	DataPoint lastData = m_dataPoints[m_dataPoints.size() - 1];
+	painter.drawEllipse(QPointF(graphWidth * (m_dataPoints.size() - 1) /
+	                                std::max(1000.0f, static_cast<float>(m_dataPoints.size())),
+	                            -graphHeight * lastData.stackedCost(5) / m_maxCost),
+	                    2, 2);
+
+	// labels
+	style()->standardPalette();
+	painter.setPen(QPen(style()->standardPalette().windowText(), 1));
+	painter.drawText(0, -graphHeight * lastData.stackedCost(1) / m_maxCost, graphWidth,
+	                 graphHeight * lastData.m_obstacle_cost / m_maxCost,
+	                 Qt::AlignRight | Qt::AlignVCenter, "obs");
+	painter.drawText(0, -graphHeight * lastData.stackedCost(2) / m_maxCost, graphWidth,
+	                 graphHeight * lastData.m_smoothing_cost / m_maxCost,
+	                 Qt::AlignRight | Qt::AlignVCenter, "sm");
+	painter.drawText(0, -graphHeight * lastData.stackedCost(3) / m_maxCost, graphWidth,
+	                 graphHeight * lastData.m_angle_restriction_cost / m_maxCost,
+	                 Qt::AlignRight | Qt::AlignVCenter, "AR");
+	painter.drawText(0, -graphHeight * lastData.stackedCost(4) / m_maxCost, graphWidth,
+	                 graphHeight * lastData.m_balancing_cost / m_maxCost,
+	                 Qt::AlignRight | Qt::AlignVCenter, "bal");
+	painter.drawText(0, -graphHeight * lastData.stackedCost(5) / m_maxCost, graphWidth,
+	                 graphHeight * lastData.m_straightening_cost / m_maxCost,
+	                 Qt::AlignRight | Qt::AlignVCenter, "str");
+}
+
+QPainterPath CostGraph::createDataPath(int costIndex, int graphWidth, int graphHeight) const {
+	QPainterPath dataPath;
+	dataPath.moveTo(0, -graphHeight * m_dataPoints[0].stackedCost(costIndex) / m_maxCost);
+	for (int i = 1; i < m_dataPoints.size(); ++i) {
+		dataPath.lineTo(graphWidth * i / std::max(1000.0f, static_cast<float>(m_dataPoints.size())),
+		                -graphHeight * m_dataPoints[i].stackedCost(costIndex) / m_maxCost);
+	}
+	return dataPath;
 }
diff --git a/demos/flow_map/optimization_demo/cost_graph.h b/demos/flow_map/optimization_demo/cost_graph.h
index ee020c8c7..675bbfe3e 100644
--- a/demos/flow_map/optimization_demo/cost_graph.h
+++ b/demos/flow_map/optimization_demo/cost_graph.h
@@ -24,6 +24,7 @@ along with this program.  If not, see <https://www.gnu.org/licenses/>.
 
 #include <QWidget>
 
+#include <array>
 #include <vector>
 
 #include "cartocrow/core/core.h"
@@ -43,6 +44,18 @@ class CostGraph : public QWidget {
 		Number<Inexact> m_angle_restriction_cost;
 		Number<Inexact> m_balancing_cost;
 		Number<Inexact> m_straightening_cost;
+
+		inline Number<Inexact> stackedCost(int count) const {
+			std::array<Number<Inexact> DataPoint::*, 5> stackingOrder = {
+			    &DataPoint::m_obstacle_cost, &DataPoint::m_smoothing_cost,
+			    &DataPoint::m_angle_restriction_cost, &DataPoint::m_balancing_cost,
+			    &DataPoint::m_straightening_cost};
+			Number<Inexact> sum = 0;
+			for (int i = 0; i < count; ++i) {
+				sum += this->*stackingOrder[i];
+			}
+			return sum;
+		}
 	};
 	/// Creates a cost graph without any cost data.
 	CostGraph();
@@ -59,6 +72,8 @@ class CostGraph : public QWidget {
 	std::vector<DataPoint> m_dataPoints;
 	/// The maximum total cost we've seen so far.
 	Number<Inexact> m_maxCost = 0;
+
+	QPainterPath createDataPath(int costIndex, int graphWidth, int graphHeight) const;
 };
 
 #endif
diff --git a/demos/flow_map/optimization_demo/optimization_demo.cpp b/demos/flow_map/optimization_demo/optimization_demo.cpp
index dd0c9c7b5..a8e05f7dc 100644
--- a/demos/flow_map/optimization_demo/optimization_demo.cpp
+++ b/demos/flow_map/optimization_demo/optimization_demo.cpp
@@ -157,7 +157,9 @@ void OptimizationDemo::recalculate() {
 }
 
 void OptimizationDemo::updateCostLabel() {
-	m_costGraph->addStep({m_smoothTree->computeCost(), 0, 0, 0, 0});
+	m_costGraph->addStep(
+	    {0, m_smoothTree->computeSmoothingCost(), m_smoothTree->computeAngleRestrictionCost(),
+	     m_smoothTree->computeBalancingCost(), m_smoothTree->computeStraighteningCost()});
 	m_costLabel->setText(QString("Iteration %1   |   Cost function: %2")
 	                         .arg(m_iterationCount)
 	                         .arg(m_smoothTree->computeCost()));

From ce8bf02b7a6581e703df9c1b9c5779a580b84d6d Mon Sep 17 00:00:00 2001
From: Willem Sonke <willem@wimiso.nl>
Date: Thu, 21 Sep 2023 11:46:45 +0200
Subject: [PATCH 23/27] Add colors to the cost graph in the demo

---
 .../flow_map/optimization_demo/cost_graph.cpp | 48 +++++++++++++++++--
 demos/flow_map/optimization_demo/cost_graph.h |  1 +
 2 files changed, 45 insertions(+), 4 deletions(-)

diff --git a/demos/flow_map/optimization_demo/cost_graph.cpp b/demos/flow_map/optimization_demo/cost_graph.cpp
index 7369d4ca2..2fc49896f 100644
--- a/demos/flow_map/optimization_demo/cost_graph.cpp
+++ b/demos/flow_map/optimization_demo/cost_graph.cpp
@@ -65,12 +65,41 @@ void CostGraph::paintEvent(QPaintEvent* event) {
 	painter.drawLine(0, 0, 0, -graphHeight);
 	painter.drawLine(0, 0, graphWidth, 0);
 
-	// graphs
+	QColor green(52, 140, 80);
+	QColor orange(255, 120, 0);
+	QColor purple(110, 60, 190);
+	QColor blue(33, 142, 252);
+	QColor red(213, 0, 74);
+
+	QPainterPath zeroPath;
+	zeroPath.moveTo(0, 0);
+	zeroPath.lineTo(graphWidth * m_dataPoints.size() /
+	                    std::max(100.0f, static_cast<float>(m_dataPoints.size())),
+	                0);
 	QPainterPath obstaclePath = createDataPath(1, graphWidth, graphHeight);
 	QPainterPath smoothingPath = createDataPath(2, graphWidth, graphHeight);
 	QPainterPath angleRestrictionPath = createDataPath(3, graphWidth, graphHeight);
 	QPainterPath balancingPath = createDataPath(4, graphWidth, graphHeight);
 	QPainterPath straighteningPath = createDataPath(5, graphWidth, graphHeight);
+
+	// shading
+	painter.setPen(Qt::NoPen);
+	painter.setOpacity(0.2);
+	painter.setBrush(green);
+	painter.drawPath(createRegionBetween(zeroPath, obstaclePath));
+	painter.setBrush(orange);
+	painter.drawPath(createRegionBetween(obstaclePath, smoothingPath));
+	painter.setBrush(purple);
+	painter.drawPath(createRegionBetween(smoothingPath, angleRestrictionPath));
+	painter.setBrush(blue);
+	painter.drawPath(createRegionBetween(angleRestrictionPath, balancingPath));
+	painter.setBrush(red);
+	painter.drawPath(createRegionBetween(balancingPath, straighteningPath));
+	painter.setOpacity(1);
+
+	// graphs
+	painter.setPen(QPen(QColor(0, 0, 0), 1));
+	painter.setBrush(Qt::NoBrush);
 	painter.drawPath(obstaclePath);
 	painter.drawPath(smoothingPath);
 	painter.drawPath(angleRestrictionPath);
@@ -82,25 +111,29 @@ void CostGraph::paintEvent(QPaintEvent* event) {
 	// dot
 	DataPoint lastData = m_dataPoints[m_dataPoints.size() - 1];
 	painter.drawEllipse(QPointF(graphWidth * (m_dataPoints.size() - 1) /
-	                                std::max(1000.0f, static_cast<float>(m_dataPoints.size())),
+	                                std::max(100.0f, static_cast<float>(m_dataPoints.size())),
 	                            -graphHeight * lastData.stackedCost(5) / m_maxCost),
 	                    2, 2);
 
 	// labels
 	style()->standardPalette();
-	painter.setPen(QPen(style()->standardPalette().windowText(), 1));
+	painter.setPen(green);
 	painter.drawText(0, -graphHeight * lastData.stackedCost(1) / m_maxCost, graphWidth,
 	                 graphHeight * lastData.m_obstacle_cost / m_maxCost,
 	                 Qt::AlignRight | Qt::AlignVCenter, "obs");
+	painter.setPen(orange);
 	painter.drawText(0, -graphHeight * lastData.stackedCost(2) / m_maxCost, graphWidth,
 	                 graphHeight * lastData.m_smoothing_cost / m_maxCost,
 	                 Qt::AlignRight | Qt::AlignVCenter, "sm");
+	painter.setPen(purple);
 	painter.drawText(0, -graphHeight * lastData.stackedCost(3) / m_maxCost, graphWidth,
 	                 graphHeight * lastData.m_angle_restriction_cost / m_maxCost,
 	                 Qt::AlignRight | Qt::AlignVCenter, "AR");
+	painter.setPen(blue);
 	painter.drawText(0, -graphHeight * lastData.stackedCost(4) / m_maxCost, graphWidth,
 	                 graphHeight * lastData.m_balancing_cost / m_maxCost,
 	                 Qt::AlignRight | Qt::AlignVCenter, "bal");
+	painter.setPen(red);
 	painter.drawText(0, -graphHeight * lastData.stackedCost(5) / m_maxCost, graphWidth,
 	                 graphHeight * lastData.m_straightening_cost / m_maxCost,
 	                 Qt::AlignRight | Qt::AlignVCenter, "str");
@@ -110,8 +143,15 @@ QPainterPath CostGraph::createDataPath(int costIndex, int graphWidth, int graphH
 	QPainterPath dataPath;
 	dataPath.moveTo(0, -graphHeight * m_dataPoints[0].stackedCost(costIndex) / m_maxCost);
 	for (int i = 1; i < m_dataPoints.size(); ++i) {
-		dataPath.lineTo(graphWidth * i / std::max(1000.0f, static_cast<float>(m_dataPoints.size())),
+		dataPath.lineTo(graphWidth * i / std::max(100.0f, static_cast<float>(m_dataPoints.size())),
 		                -graphHeight * m_dataPoints[i].stackedCost(costIndex) / m_maxCost);
 	}
 	return dataPath;
 }
+
+QPainterPath CostGraph::createRegionBetween(const QPainterPath& first,
+                                            const QPainterPath& second) const {
+	QPainterPath region = first;
+	region.connectPath(second.toReversed());
+	return region;
+}
diff --git a/demos/flow_map/optimization_demo/cost_graph.h b/demos/flow_map/optimization_demo/cost_graph.h
index 675bbfe3e..2cd39e2a2 100644
--- a/demos/flow_map/optimization_demo/cost_graph.h
+++ b/demos/flow_map/optimization_demo/cost_graph.h
@@ -74,6 +74,7 @@ class CostGraph : public QWidget {
 	Number<Inexact> m_maxCost = 0;
 
 	QPainterPath createDataPath(int costIndex, int graphWidth, int graphHeight) const;
+	QPainterPath createRegionBetween(const QPainterPath& first, const QPainterPath& second) const;
 };
 
 #endif

From 261e9ce71d3fb3e4d3f825ef1870b9b5f25c2e93 Mon Sep 17 00:00:00 2001
From: Willem Sonke <willem@wimiso.nl>
Date: Thu, 21 Sep 2023 14:26:27 +0200
Subject: [PATCH 24/27] Stop optimizing when encountering a nan

---
 cartocrow/flow_map/smooth_tree.cpp                  |  4 ++--
 demos/flow_map/optimization_demo/cost_graph.cpp     |  3 ++-
 demos/flow_map/optimization_demo/cost_graph.h       |  4 +++-
 .../optimization_demo/optimization_demo.cpp         | 13 ++++++++++---
 4 files changed, 17 insertions(+), 7 deletions(-)

diff --git a/cartocrow/flow_map/smooth_tree.cpp b/cartocrow/flow_map/smooth_tree.cpp
index 092b26ccf..945482638 100644
--- a/cartocrow/flow_map/smooth_tree.cpp
+++ b/cartocrow/flow_map/smooth_tree.cpp
@@ -309,7 +309,7 @@ void SmoothTree::optimize() {
 		if (node->getType() == Node::ConnectionType::kSubdivision) {
 			applySmoothingGradient(i, m_nodes[i]->m_parent->m_id, m_nodes[i]->m_children[0]->m_id);
 		} else if (node->getType() == Node::ConnectionType::kJoin) {
-			/*applyAngleRestrictionGradient(
+			applyAngleRestrictionGradient(
 			    i, m_nodes[i]->m_children[0]->m_id,
 			    m_nodes[i]->m_children[m_nodes[i]->m_children.size() - 1]->m_id);
 			applyBalancingGradient(
@@ -317,7 +317,7 @@ void SmoothTree::optimize() {
 			    m_nodes[i]->m_children[m_nodes[i]->m_children.size() - 1]->m_id);
 			applyStraighteningGradient(
 			    i, m_nodes[i]->m_parent->m_id,
-			    m_nodes[i]->m_children);*/
+			    m_nodes[i]->m_children);
 		}
 	}
 	for (int i = 0; i < m_nodes.size(); i++) {
diff --git a/demos/flow_map/optimization_demo/cost_graph.cpp b/demos/flow_map/optimization_demo/cost_graph.cpp
index 2fc49896f..71aa36d6a 100644
--- a/demos/flow_map/optimization_demo/cost_graph.cpp
+++ b/demos/flow_map/optimization_demo/cost_graph.cpp
@@ -25,6 +25,7 @@ along with this program.  If not, see <https://www.gnu.org/licenses/>.
 #include <QPainterPath>
 #include <QPalette>
 #include <QStyle>
+#include <cmath>
 
 using namespace cartocrow;
 
@@ -142,7 +143,7 @@ void CostGraph::paintEvent(QPaintEvent* event) {
 QPainterPath CostGraph::createDataPath(int costIndex, int graphWidth, int graphHeight) const {
 	QPainterPath dataPath;
 	dataPath.moveTo(0, -graphHeight * m_dataPoints[0].stackedCost(costIndex) / m_maxCost);
-	for (int i = 1; i < m_dataPoints.size(); ++i) {
+	for (int i = 1; i < m_dataPoints.size(); i += 1 + i / width()) {
 		dataPath.lineTo(graphWidth * i / std::max(100.0f, static_cast<float>(m_dataPoints.size())),
 		                -graphHeight * m_dataPoints[i].stackedCost(costIndex) / m_maxCost);
 	}
diff --git a/demos/flow_map/optimization_demo/cost_graph.h b/demos/flow_map/optimization_demo/cost_graph.h
index 2cd39e2a2..bb0b9e6fd 100644
--- a/demos/flow_map/optimization_demo/cost_graph.h
+++ b/demos/flow_map/optimization_demo/cost_graph.h
@@ -52,7 +52,9 @@ class CostGraph : public QWidget {
 			    &DataPoint::m_straightening_cost};
 			Number<Inexact> sum = 0;
 			for (int i = 0; i < count; ++i) {
-				sum += this->*stackingOrder[i];
+				if (!std::isnan(this->*stackingOrder[i])) {
+					sum += this->*stackingOrder[i];
+				}
 			}
 			return sum;
 		}
diff --git a/demos/flow_map/optimization_demo/optimization_demo.cpp b/demos/flow_map/optimization_demo/optimization_demo.cpp
index a8e05f7dc..963cf037a 100644
--- a/demos/flow_map/optimization_demo/optimization_demo.cpp
+++ b/demos/flow_map/optimization_demo/optimization_demo.cpp
@@ -89,6 +89,11 @@ OptimizationDemo::OptimizationDemo() {
 			m_iterationCount++;
 			m_smoothTree->optimize();
 			updateCostLabel();
+			if (std::isnan(m_smoothTree->computeCost())){
+				m_optimizeButton->setChecked(false);
+				m_optimizeTimer->stop();
+				break;
+			}
 		}
 		m_renderer->update();
 	});
@@ -160,9 +165,11 @@ void OptimizationDemo::updateCostLabel() {
 	m_costGraph->addStep(
 	    {0, m_smoothTree->computeSmoothingCost(), m_smoothTree->computeAngleRestrictionCost(),
 	     m_smoothTree->computeBalancingCost(), m_smoothTree->computeStraighteningCost()});
-	m_costLabel->setText(QString("Iteration %1   |   Cost function: %2")
-	                         .arg(m_iterationCount)
-	                         .arg(m_smoothTree->computeCost()));
+	Number<Inexact> cost = m_smoothTree->computeCost();
+	m_costLabel->setText(
+	    QString("Iteration %1 | Cost function: %2")
+	        .arg(m_iterationCount)
+	        .arg(std::isnan(cost) ? "<b><font color=\"#d5004a\">nan</font></b>" : QString::number(cost)));
 }
 
 int main(int argc, char* argv[]) {

From 1026480ac3b3913286e740b6d0d7265c12e247f0 Mon Sep 17 00:00:00 2001
From: Willem Sonke <willem@wimiso.nl>
Date: Thu, 21 Sep 2023 14:52:03 +0200
Subject: [PATCH 25/27] Add a toggle for the "stop on nan" behavior

---
 demos/flow_map/optimization_demo/optimization_demo.cpp | 5 ++++-
 demos/flow_map/optimization_demo/optimization_demo.h   | 2 ++
 2 files changed, 6 insertions(+), 1 deletion(-)

diff --git a/demos/flow_map/optimization_demo/optimization_demo.cpp b/demos/flow_map/optimization_demo/optimization_demo.cpp
index 963cf037a..7ba401ed9 100644
--- a/demos/flow_map/optimization_demo/optimization_demo.cpp
+++ b/demos/flow_map/optimization_demo/optimization_demo.cpp
@@ -25,6 +25,7 @@ along with this program.  If not, see <https://www.gnu.org/licenses/>.
 #include <QCheckBox>
 #include <QDockWidget>
 #include <QStatusBar>
+#include <qcheckbox.h>
 
 #include "cartocrow/core/core.h"
 #include "cartocrow/flow_map/painting.h"
@@ -89,7 +90,7 @@ OptimizationDemo::OptimizationDemo() {
 			m_iterationCount++;
 			m_smoothTree->optimize();
 			updateCostLabel();
-			if (std::isnan(m_smoothTree->computeCost())){
+			if (m_stopOnNanCheckbox->isChecked() && std::isnan(m_smoothTree->computeCost())) {
 				m_optimizeButton->setChecked(false);
 				m_optimizeTimer->stop();
 				break;
@@ -97,6 +98,8 @@ OptimizationDemo::OptimizationDemo() {
 		}
 		m_renderer->update();
 	});
+	m_stopOnNanCheckbox = new QCheckBox("Stop on nan");
+	toolBar->addWidget(m_stopOnNanCheckbox);
 	m_optimizeOneStepButton = new QPushButton("Optimize one step");
 	toolBar->addWidget(m_optimizeOneStepButton);
 	connect(m_optimizeOneStepButton, &QPushButton::clicked, [&]() {
diff --git a/demos/flow_map/optimization_demo/optimization_demo.h b/demos/flow_map/optimization_demo/optimization_demo.h
index 9b1cd6f72..37ddb8805 100644
--- a/demos/flow_map/optimization_demo/optimization_demo.h
+++ b/demos/flow_map/optimization_demo/optimization_demo.h
@@ -19,6 +19,7 @@ You should have received a copy of the GNU General Public License
 along with this program.  If not, see <https://www.gnu.org/licenses/>.
 */
 
+#include <QCheckBox>
 #include <QLabel>
 #include <QMainWindow>
 #include <QPushButton>
@@ -56,6 +57,7 @@ class OptimizationDemo : public QMainWindow {
 	QSlider* m_alphaSlider;
 	QLabel* m_alphaLabel;
 	QPushButton* m_optimizeButton;
+	QCheckBox* m_stopOnNanCheckbox;
 	QPushButton* m_optimizeOneStepButton;
 	QTimer* m_optimizeTimer;
 	QLabel* m_costLabel;

From ab5d0c00e132f90d7745210972ec932aa0a9021d Mon Sep 17 00:00:00 2001
From: Willem Sonke <willem@wimiso.nl>
Date: Fri, 22 Sep 2023 16:23:32 +0200
Subject: [PATCH 26/27] =?UTF-8?q?Disable=20=CF=86=20randomization=20in=20t?=
 =?UTF-8?q?he=20demo?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

---
 demos/flow_map/optimization_demo/optimization_demo.cpp | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/demos/flow_map/optimization_demo/optimization_demo.cpp b/demos/flow_map/optimization_demo/optimization_demo.cpp
index 7ba401ed9..136f27454 100644
--- a/demos/flow_map/optimization_demo/optimization_demo.cpp
+++ b/demos/flow_map/optimization_demo/optimization_demo.cpp
@@ -149,13 +149,13 @@ void OptimizationDemo::recalculate() {
 	ReachableRegionAlgorithm::ReachableRegion reachableRegion = ReachableRegionAlgorithm(tree).run();
 	SpiralTreeObstructedAlgorithm(tree, reachableRegion).run();
 	m_smoothTree = std::make_shared<SmoothTree>(tree);
-	for (auto& node : m_smoothTree->m_nodes) {
+	/*for (auto& node : m_smoothTree->m_nodes) {
 		if (node->getType() == Node::ConnectionType::kSubdivision) {
 			node->m_position.setPhi(wrapAngle(
 			    node->m_position.phi() + static_cast<float>(std::rand()) / RAND_MAX * 0.2 - 0.1,
 			    -M_PI));
 		}
-	}
+	}*/
 	SmoothTreePainting::Options options;
 	auto painting = std::make_shared<SmoothTreePainting>(m_smoothTree, options);
 	m_renderer->addPainting(painting, "Smooth tree");

From 1eb7eabcc606ae9e22be9cd7bc02a003b9db19a8 Mon Sep 17 00:00:00 2001
From: Willem Sonke <willem@wimiso.nl>
Date: Fri, 22 Sep 2023 16:33:59 +0200
Subject: [PATCH 27/27] Fix some small code style issues

---
 cartocrow/flow_map/smooth_tree.cpp            | 56 +++++++++----------
 cartocrow/flow_map/smooth_tree.h              | 10 ++--
 cartocrow/flow_map/smooth_tree_painting.cpp   |  3 +-
 cartocrow/renderer/geometry_widget.cpp        |  1 -
 .../optimization_demo/optimization_demo.cpp   | 11 ----
 5 files changed, 35 insertions(+), 46 deletions(-)

diff --git a/cartocrow/flow_map/smooth_tree.cpp b/cartocrow/flow_map/smooth_tree.cpp
index 945482638..4ddd29300 100644
--- a/cartocrow/flow_map/smooth_tree.cpp
+++ b/cartocrow/flow_map/smooth_tree.cpp
@@ -96,7 +96,7 @@ Number<Inexact> SmoothTree::computeSmoothingCost(int i, int iParent, int iChild)
 	PolarPoint n = m_nodes[i]->m_position;
 	PolarPoint p = m_nodes[iParent]->m_position;
 	PolarPoint c = m_nodes[iChild]->m_position;
-	return m_smoothing_factor * std::pow(Spiral::alpha(p, n) - Spiral::alpha(n, c), 2);
+	return m_smoothingFactor * std::pow(Spiral::alpha(p, n) - Spiral::alpha(n, c), 2);
 }
 
 void SmoothTree::applySmoothingGradient(int i, int iParent, int iChild) {
@@ -104,19 +104,19 @@ void SmoothTree::applySmoothingGradient(int i, int iParent, int iChild) {
 	PolarPoint p = m_nodes[iParent]->m_position;
 	PolarPoint c = m_nodes[iChild]->m_position;
 
-	m_gradient[i].r += 2 * m_smoothing_factor * (Spiral::alpha(p, n) - Spiral::alpha(n, c)) *
+	m_gradient[i].r += 2 * m_smoothingFactor * (Spiral::alpha(p, n) - Spiral::alpha(n, c)) *
 	                   (Spiral::dAlphaDR2(p, n) - Spiral::dAlphaDR1(n, c));
-	m_gradient[i].phi += 2 * m_smoothing_factor * (Spiral::alpha(p, n) - Spiral::alpha(n, c)) *
+	m_gradient[i].phi += 2 * m_smoothingFactor * (Spiral::alpha(p, n) - Spiral::alpha(n, c)) *
 	                     (Spiral::dAlphaDPhi2(p, n) - Spiral::dAlphaDPhi1(n, c));
 
-	m_gradient[iParent].r += 2 * m_smoothing_factor * (Spiral::alpha(p, n) - Spiral::alpha(n, c)) *
+	m_gradient[iParent].r += 2 * m_smoothingFactor * (Spiral::alpha(p, n) - Spiral::alpha(n, c)) *
 	                         Spiral::dAlphaDR1(p, n);
-	m_gradient[iParent].phi += 2 * m_smoothing_factor * (Spiral::alpha(p, n) - Spiral::alpha(n, c)) *
+	m_gradient[iParent].phi += 2 * m_smoothingFactor * (Spiral::alpha(p, n) - Spiral::alpha(n, c)) *
 	                           Spiral::dAlphaDPhi1(p, n);
 
-	m_gradient[iChild].r += 2 * m_smoothing_factor * (Spiral::alpha(p, n) - Spiral::alpha(n, c)) *
+	m_gradient[iChild].r += 2 * m_smoothingFactor * (Spiral::alpha(p, n) - Spiral::alpha(n, c)) *
 	                        -Spiral::dAlphaDR2(n, c);
-	m_gradient[iChild].phi += 2 * m_smoothing_factor * (Spiral::alpha(p, n) - Spiral::alpha(n, c)) *
+	m_gradient[iChild].phi += 2 * m_smoothingFactor * (Spiral::alpha(p, n) - Spiral::alpha(n, c)) *
 	                          -Spiral::dAlphaDPhi2(n, c);
 }
 
@@ -137,7 +137,7 @@ Number<Inexact> SmoothTree::computeAngleRestrictionCost(int i, int iChild1, int
 	PolarPoint n = m_nodes[i]->m_position;
 	PolarPoint c1 = m_nodes[iChild1]->m_position;
 	PolarPoint c2 = m_nodes[iChild2]->m_position;
-	return m_angle_restriction_factor * (std::log(1.0 / std::cos(Spiral::alpha(n, c1))) +
+	return m_angle_restrictionFactor * (std::log(1.0 / std::cos(Spiral::alpha(n, c1))) +
 	                                     std::log(1.0 / std::cos(Spiral::alpha(n, c2))));
 }
 
@@ -147,21 +147,21 @@ void SmoothTree::applyAngleRestrictionGradient(int i, int iChild1, int iChild2)
 	PolarPoint c2 = m_nodes[iChild2]->m_position;
 
 	m_gradient[i].r +=
-	    m_angle_restriction_factor * (Spiral::dAlphaDR1(n, c1) * std::tan(Spiral::alpha(n, c1)) +
+	    m_angle_restrictionFactor * (Spiral::dAlphaDR1(n, c1) * std::tan(Spiral::alpha(n, c1)) +
 	                                  Spiral::dAlphaDR1(n, c2) * std::tan(Spiral::alpha(n, c2)));
 	m_gradient[i].phi +=
-	    m_angle_restriction_factor * (Spiral::dAlphaDPhi1(n, c1) * std::tan(Spiral::alpha(n, c1)) +
+	    m_angle_restrictionFactor * (Spiral::dAlphaDPhi1(n, c1) * std::tan(Spiral::alpha(n, c1)) +
 	                                  Spiral::dAlphaDPhi1(n, c2) * std::tan(Spiral::alpha(n, c2)));
 
 	m_gradient[iChild1].r +=
-	    m_angle_restriction_factor * Spiral::dAlphaDR2(n, c1) * std::tan(Spiral::alpha(n, c1));
+	    m_angle_restrictionFactor * Spiral::dAlphaDR2(n, c1) * std::tan(Spiral::alpha(n, c1));
 	m_gradient[iChild1].phi +=
-	    m_angle_restriction_factor * Spiral::dAlphaDPhi2(n, c1) * std::tan(Spiral::alpha(n, c1));
+	    m_angle_restrictionFactor * Spiral::dAlphaDPhi2(n, c1) * std::tan(Spiral::alpha(n, c1));
 
 	m_gradient[iChild2].r +=
-	    m_angle_restriction_factor * Spiral::dAlphaDR2(n, c2) * std::tan(Spiral::alpha(n, c2));
+	    m_angle_restrictionFactor * Spiral::dAlphaDR2(n, c2) * std::tan(Spiral::alpha(n, c2));
 	m_gradient[iChild2].phi +=
-	    m_angle_restriction_factor * Spiral::dAlphaDPhi2(n, c2) * std::tan(Spiral::alpha(n, c2));
+	    m_angle_restrictionFactor * Spiral::dAlphaDPhi2(n, c2) * std::tan(Spiral::alpha(n, c2));
 }
 
 Number<Inexact> SmoothTree::computeBalancingCost() {
@@ -181,7 +181,7 @@ Number<Inexact> SmoothTree::computeBalancingCost(int i, int iChild1, int iChild2
 	PolarPoint n = m_nodes[i]->m_position;
 	PolarPoint c1 = m_nodes[iChild1]->m_position;
 	PolarPoint c2 = m_nodes[iChild2]->m_position;
-	return m_angle_restriction_factor * 2 * std::pow(std::tan(m_restrictingAngle), 2) *
+	return m_angle_restrictionFactor * 2 * std::pow(std::tan(m_restrictingAngle), 2) *
 	       std::log(1 / std::sin(0.5 * (Spiral::alpha(n, c1) - Spiral::alpha(n, c2))));
 }
 
@@ -190,25 +190,25 @@ void SmoothTree::applyBalancingGradient(int i, int iChild1, int iChild2) {
 	PolarPoint c1 = m_nodes[iChild1]->m_position;
 	PolarPoint c2 = m_nodes[iChild2]->m_position;
 
-	m_gradient[i].r += m_angle_restriction_factor * -std::pow(std::tan(m_restrictingAngle), 2) *
+	m_gradient[i].r += m_angle_restrictionFactor * -std::pow(std::tan(m_restrictingAngle), 2) *
 	                   (1 / std::tan(0.5 * (Spiral::alpha(n, c1) - Spiral::alpha(n, c2)))) *
 	                   (Spiral::dAlphaDR1(n, c1) - Spiral::dAlphaDR1(n, c2));
-	m_gradient[i].phi += m_angle_restriction_factor * -std::pow(std::tan(m_restrictingAngle), 2) *
+	m_gradient[i].phi += m_angle_restrictionFactor * -std::pow(std::tan(m_restrictingAngle), 2) *
 	                     (1 / std::tan(0.5 * (Spiral::alpha(n, c1) - Spiral::alpha(n, c2)))) *
 	                     (Spiral::dAlphaDPhi1(n, c1) - Spiral::dAlphaDPhi1(n, c2));
 
-	m_gradient[iChild1].r += m_angle_restriction_factor * -std::pow(std::tan(m_restrictingAngle), 2) *
+	m_gradient[iChild1].r += m_angle_restrictionFactor * -std::pow(std::tan(m_restrictingAngle), 2) *
 	                         (1 / std::tan(0.5 * (Spiral::alpha(n, c1) - Spiral::alpha(n, c2)))) *
 	                         Spiral::dAlphaDR2(n, c1);
-	m_gradient[iChild1].phi += m_angle_restriction_factor *
+	m_gradient[iChild1].phi += m_angle_restrictionFactor *
 	                           -std::pow(std::tan(m_restrictingAngle), 2) *
 	                           (1 / std::tan(0.5 * (Spiral::alpha(n, c1) - Spiral::alpha(n, c2)))) *
 	                           Spiral::dAlphaDPhi2(n, c1);
 
-	m_gradient[iChild2].r += m_angle_restriction_factor * -std::pow(std::tan(m_restrictingAngle), 2) *
+	m_gradient[iChild2].r += m_angle_restrictionFactor * -std::pow(std::tan(m_restrictingAngle), 2) *
 	                         (1 / std::tan(0.5 * (Spiral::alpha(n, c1) - Spiral::alpha(n, c2)))) *
 	                         -Spiral::dAlphaDR2(n, c2);
-	m_gradient[iChild2].phi += m_angle_restriction_factor *
+	m_gradient[iChild2].phi += m_angle_restrictionFactor *
 	                           -std::pow(std::tan(m_restrictingAngle), 2) *
 	                           (1 / std::tan(0.5 * (Spiral::alpha(n, c1) - Spiral::alpha(n, c2)))) *
 	                           -Spiral::dAlphaDPhi2(n, c2);
@@ -245,7 +245,7 @@ SmoothTree::computeStraighteningCost(int i, int iParent,
 			denominator += child->m_flow;
 		}
 	}
-	return m_straightening_factor * std::pow(Spiral::alpha(p, n) - numerator / denominator, 2);
+	return m_straighteningFactor * std::pow(Spiral::alpha(p, n) - numerator / denominator, 2);
 }
 
 void SmoothTree::applyStraighteningGradient(int i, int iParent,
@@ -271,26 +271,26 @@ void SmoothTree::applyStraighteningGradient(int i, int iParent,
 			denominator += child->m_flow;
 		}
 	}
-	m_gradient[i].r += 2 * m_straightening_factor * (Spiral::alpha(p, n) - numerator / denominator) *
+	m_gradient[i].r += 2 * m_straighteningFactor * (Spiral::alpha(p, n) - numerator / denominator) *
 	                   (Spiral::dAlphaDR2(p, n) - numeratorDR1 / denominator);
-	m_gradient[i].phi += 2 * m_straightening_factor *
+	m_gradient[i].phi += 2 * m_straighteningFactor *
 	                     (Spiral::alpha(p, n) - numerator / denominator) *
 	                     (Spiral::dAlphaDPhi2(p, n) - numeratorDPhi1 / denominator);
 
-	m_gradient[iParent].r += 2 * m_straightening_factor *
+	m_gradient[iParent].r += 2 * m_straighteningFactor *
 	                         (Spiral::alpha(p, n) - numerator / denominator) *
 	                         Spiral::dAlphaDR1(p, n);
-	m_gradient[iParent].phi += 2 * m_straightening_factor *
+	m_gradient[iParent].phi += 2 * m_straighteningFactor *
 	                           (Spiral::alpha(p, n) - numerator / denominator) *
 	                           Spiral::dAlphaDPhi1(p, n);
 
 	for (const auto& child : children) {
 		if (child->m_flow > maxFlow) {
 			PolarPoint c = child->m_position;
-			m_gradient[child->m_id].r += 2 * m_straightening_factor *
+			m_gradient[child->m_id].r += 2 * m_straighteningFactor *
 			                             (Spiral::alpha(p, n) - numerator / denominator) *
 			                             -child->m_flow * Spiral::dAlphaDR2(n, c) / denominator;
-			m_gradient[child->m_id].phi += 2 * m_straightening_factor *
+			m_gradient[child->m_id].phi += 2 * m_straighteningFactor *
 			                               (Spiral::alpha(p, n) - numerator / denominator) *
 			                               -child->m_flow * Spiral::dAlphaDPhi2(n, c) / denominator;
 		}
diff --git a/cartocrow/flow_map/smooth_tree.h b/cartocrow/flow_map/smooth_tree.h
index b40e18a7b..fce2be64b 100644
--- a/cartocrow/flow_map/smooth_tree.h
+++ b/cartocrow/flow_map/smooth_tree.h
@@ -56,7 +56,7 @@ class SmoothTree {
 	/// Performs one optimization step.
 	void optimize();
 
-  //private:  // TODO temporary
+  private:
 	/// The spiral tree underlying this smooth tree.
 	std::shared_ptr<SpiralTree> m_tree;
 
@@ -215,10 +215,10 @@ class SmoothTree {
 	void applyStraighteningGradient(int i, int iParent,
 	                                const std::vector<std::shared_ptr<Node>>& children);
 
-	Number<Inexact> m_obstacle_factor = 2.0;
-	Number<Inexact> m_smoothing_factor = 0.4;
-	Number<Inexact> m_straightening_factor = 0.4;
-	Number<Inexact> m_angle_restriction_factor = 0.077;
+	Number<Inexact> m_obstacleFactor = 2.0;
+	Number<Inexact> m_smoothingFactor = 0.4;
+	Number<Inexact> m_straighteningFactor = 0.4;
+	Number<Inexact> m_angle_restrictionFactor = 0.077;
 
 	Number<Inexact> m_relevantFlowFactor = 0.5;
 };
diff --git a/cartocrow/flow_map/smooth_tree_painting.cpp b/cartocrow/flow_map/smooth_tree_painting.cpp
index e1fd06f4e..61b7573cc 100644
--- a/cartocrow/flow_map/smooth_tree_painting.cpp
+++ b/cartocrow/flow_map/smooth_tree_painting.cpp
@@ -1,8 +1,9 @@
 #include "smooth_tree_painting.h"
 
+#include <CGAL/Origin.h>
+
 #include "../core/core.h"
 #include "../renderer/geometry_renderer.h"
-#include <CGAL/Origin.h>
 
 namespace cartocrow::flow_map {
 
diff --git a/cartocrow/renderer/geometry_widget.cpp b/cartocrow/renderer/geometry_widget.cpp
index d482f255f..86f9abf6f 100644
--- a/cartocrow/renderer/geometry_widget.cpp
+++ b/cartocrow/renderer/geometry_widget.cpp
@@ -31,7 +31,6 @@ along with this program.  If not, see <https://www.gnu.org/licenses/>.
 
 #include <cmath>
 #include <limits>
-#include <qnamespace.h>
 
 namespace cartocrow::renderer {
 
diff --git a/demos/flow_map/optimization_demo/optimization_demo.cpp b/demos/flow_map/optimization_demo/optimization_demo.cpp
index 136f27454..20ffc86e0 100644
--- a/demos/flow_map/optimization_demo/optimization_demo.cpp
+++ b/demos/flow_map/optimization_demo/optimization_demo.cpp
@@ -22,10 +22,8 @@ along with this program.  If not, see <https://www.gnu.org/licenses/>.
 #include "optimization_demo.h"
 
 #include <QApplication>
-#include <QCheckBox>
 #include <QDockWidget>
 #include <QStatusBar>
-#include <qcheckbox.h>
 
 #include "cartocrow/core/core.h"
 #include "cartocrow/flow_map/painting.h"
@@ -46,15 +44,6 @@ using namespace cartocrow::renderer;
 OptimizationDemo::OptimizationDemo() {
 	setWindowTitle("CartoCrow – Optimization demo");
 
-	/*m_places.push_back(std::make_shared<Point<Inexact>>(11.2121212, 17.0707070));
-	m_places.push_back(std::make_shared<Point<Inexact>>(13.9393939, -14.1414141));
-	m_places.push_back(std::make_shared<Point<Inexact>>(-4.5454545, -18.9898989));
-	m_places.push_back(std::make_shared<Point<Inexact>>(16.6666666, 6.1616161));
-	m_places.push_back(std::make_shared<Point<Inexact>>(-9.8989898, 13.9393939));
-	m_places.push_back(std::make_shared<Point<Inexact>>(-16.1616161, -2.6262626));*/
-	/*m_places.push_back(std::make_shared<Point<Inexact>>(2, 1));
-	m_places.push_back(std::make_shared<Point<Inexact>>(2, -1));*/
-
 	for (int i = 0; i < 40; i++) {
 		m_places.push_back(std::make_shared<Point<Inexact>>(
 		    static_cast<float>(std::rand()) / RAND_MAX * 200 - 100,