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/node.h b/cartocrow/flow_map/node.h index a5c407f4e..9381d6a2b 100644 --- a/cartocrow/flow_map/node.h +++ b/cartocrow/flow_map/node.h @@ -57,6 +57,12 @@ 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> m_children; + + /// The amount of flow through this node. + Number m_flow; + + /// The index of this node in its tree. + int m_id = -1; }; } // namespace cartocrow::flow_map 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 PolarPoint::phi() const { return m_phi; } +void PolarPoint::setR(Number r) { + m_r = r; +} + +void PolarPoint::setPhi(Number phi) { + m_phi = phi; +} + Point PolarPoint::toCartesian() const { const Vector d = Vector(std::cos(phi()), std::sin(phi())); return Point(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& R, const Number& 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& t); - /// Constructs a polar point from a point in Cartesian coordinates. explicit PolarPoint(const Point& p); - PolarPoint(const Point& p, const Vector& 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 phi() const; + /// Sets the distance \f$r\f$ from the origin. + void setR(Number r); + /// Sets the distance \f$\phi\f$ relative to the origin. + void setPhi(Number phi); + /// Returns the point in Cartesian coordinates corresponding to this polar /// point. Point toCartesian() const; diff --git a/cartocrow/flow_map/smooth_tree.cpp b/cartocrow/flow_map/smooth_tree.cpp new file mode 100644 index 000000000..4ddd29300 --- /dev/null +++ b/cartocrow/flow_map/smooth_tree.cpp @@ -0,0 +1,335 @@ +/* +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 . +*/ + +#include "smooth_tree.h" +#include "cartocrow/flow_map/smooth_tree_painting.h" +#include + +namespace cartocrow::flow_map { + +SmoothTree::SmoothTree(const std::shared_ptr spiralTree) + : m_tree(spiralTree), m_restrictingAngle(spiralTree->restrictingAngle()) { + Number 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 + + // set node IDs + for (int i = 0; i < m_nodes.size(); i++) { + m_nodes[i]->m_id = i; + } +} + +std::shared_ptr SmoothTree::constructSmoothTree(const std::shared_ptr& node, + Number maxRStep) { + auto smoothNode = std::make_shared(node->m_position); + m_nodes.push_back(smoothNode); + 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 rMin = node->m_position.r(); + Number rMax = child->m_position.r(); + Number 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(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); + } + } + + return smoothNode; +} + +const std::vector>& SmoothTree::nodes() const { + return m_nodes; +} + +Number SmoothTree::computeSmoothingCost() { + Number 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 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_smoothingFactor * std::pow(Spiral::alpha(p, n) - Spiral::alpha(n, c), 2); +} + +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_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_smoothingFactor * (Spiral::alpha(p, n) - Spiral::alpha(n, c)) * + (Spiral::dAlphaDPhi2(p, n) - Spiral::dAlphaDPhi1(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_smoothingFactor * (Spiral::alpha(p, n) - Spiral::alpha(n, c)) * + Spiral::dAlphaDPhi1(p, n); + + 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_smoothingFactor * (Spiral::alpha(p, n) - Spiral::alpha(n, c)) * + -Spiral::dAlphaDPhi2(n, c); +} + +Number SmoothTree::computeAngleRestrictionCost() { + Number 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 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_restrictionFactor * (std::log(1.0 / std::cos(Spiral::alpha(n, c1))) + + std::log(1.0 / std::cos(Spiral::alpha(n, c2)))); +} + +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_gradient[i].r += + 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_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_restrictionFactor * Spiral::dAlphaDR2(n, c1) * std::tan(Spiral::alpha(n, c1)); + m_gradient[iChild1].phi += + m_angle_restrictionFactor * Spiral::dAlphaDPhi2(n, c1) * std::tan(Spiral::alpha(n, c1)); + + m_gradient[iChild2].r += + m_angle_restrictionFactor * Spiral::dAlphaDR2(n, c2) * std::tan(Spiral::alpha(n, c2)); + m_gradient[iChild2].phi += + m_angle_restrictionFactor * Spiral::dAlphaDPhi2(n, c2) * std::tan(Spiral::alpha(n, c2)); +} + +Number SmoothTree::computeBalancingCost() { + Number 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 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_restrictionFactor * 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::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_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_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_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_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_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_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); +} + +Number SmoothTree::computeStraighteningCost() { + Number 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 +SmoothTree::computeStraighteningCost(int i, int iParent, + const std::vector>& children) { + PolarPoint n = m_nodes[i]->m_position; + PolarPoint p = m_nodes[iParent]->m_position; + Number maxFlow = 0; + for (const auto& child : children) { + if (child->m_flow > maxFlow) { + maxFlow = child->m_flow; + } + } + Number numerator = 0; + Number 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_straighteningFactor * std::pow(Spiral::alpha(p, n) - numerator / denominator, 2); +} + +void SmoothTree::applyStraighteningGradient(int i, int iParent, + const std::vector>& children) { + PolarPoint n = m_nodes[i]->m_position; + PolarPoint p = m_nodes[iParent]->m_position; + Number maxFlow = 0; + for (const auto& child : children) { + if (child->m_flow > maxFlow) { + maxFlow = child->m_flow; + } + } + Number numerator = 0; + Number numeratorDR1 = 0; + Number numeratorDPhi1 = 0; + Number 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_gradient[i].r += 2 * m_straighteningFactor * (Spiral::alpha(p, n) - numerator / denominator) * + (Spiral::dAlphaDR2(p, n) - numeratorDR1 / denominator); + 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_straighteningFactor * + (Spiral::alpha(p, n) - numerator / denominator) * + Spiral::dAlphaDR1(p, n); + 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_straighteningFactor * + (Spiral::alpha(p, n) - numerator / denominator) * + -child->m_flow * Spiral::dAlphaDR2(n, c) / denominator; + m_gradient[child->m_id].phi += 2 * m_straighteningFactor * + (Spiral::alpha(p, n) - numerator / denominator) * + -child->m_flow * Spiral::dAlphaDPhi2(n, c) / denominator; + } + } +} + +Number SmoothTree::computeCost() { + return computeSmoothingCost() + computeAngleRestrictionCost() + computeBalancingCost() + + computeStraighteningCost(); +} + +void SmoothTree::optimize() { + m_gradient = std::vector(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) { + 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( + i, m_nodes[i]->m_children[0]->m_id, + m_nodes[i]->m_children[m_nodes[i]->m_children.size() - 1]->m_id); + applyBalancingGradient( + i, m_nodes[i]->m_children[0]->m_id, + 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); + } + } + for (int i = 0; i < m_nodes.size(); i++) { + Number 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_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_gradient[i].phi); + } + } +} + +} // 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..fce2be64b --- /dev/null +++ b/cartocrow/flow_map/smooth_tree.h @@ -0,0 +1,228 @@ +/* +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 . +*/ + +#ifndef CARTOCROW_FLOW_MAP_SMOOTH_TREE_H +#define CARTOCROW_FLOW_MAP_SMOOTH_TREE_H + +#include +#include +#include +#include +#include + +#include "spiral_tree.h" + +namespace cartocrow::flow_map { + +/// A smoothed tree. +class SmoothTree { + + public: + /// Constructs a smooth tree. + SmoothTree(const std::shared_ptr spiralTree); + + /// Returns a list of the nodes in this smooth tree. + const std::vector>& nodes() const; + + /// Computes the total cost of the tree. + Number computeCost(); + /// Computes the smoothing cost of the entire tree. + Number computeSmoothingCost(); + /// Computes the angle restriction cost of the entire tree. + Number computeAngleRestrictionCost(); + /// Computes the balancing cost of the entire tree. + Number computeBalancingCost(); + /// Computes the straightening cost of the entire tree. + Number computeStraighteningCost(); + + /// Performs one optimization step. + void optimize(); + + private: + /// The spiral tree underlying this smooth tree. + std::shared_ptr m_tree; + + Number m_restrictingAngle; + + /// List of nodes in this tree. + std::vector> m_nodes; + + std::shared_ptr constructSmoothTree(const std::shared_ptr& node, Number maxRStep); + + struct PolarGradient { + Number r = 0; + Number phi = 0; + }; + std::vector 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 + /// `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 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 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 + /// `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 computeAngleRestrictionCost(int i, int iChild1, int 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: + /// + /// \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 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$. + /// + /// \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 computeBalancingCost(int i, int iChild1, int 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: + /// + /// \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 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$. + /// + /// \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 computeStraighteningCost(int i, int iParent, + const std::vector>& 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 applyStraighteningGradient(int i, int iParent, + const std::vector>& children); + + Number m_obstacleFactor = 2.0; + Number m_smoothingFactor = 0.4; + Number m_straighteningFactor = 0.4; + Number m_angle_restrictionFactor = 0.077; + + Number m_relevantFlowFactor = 0.5; +}; + +} // 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..61b7573cc --- /dev/null +++ b/cartocrow/flow_map/smooth_tree_painting.cpp @@ -0,0 +1,57 @@ +#include "smooth_tree_painting.h" + +#include + +#include "../core/core.h" +#include "../renderer/geometry_renderer.h" + +namespace cartocrow::flow_map { + +SmoothTreePainting::SmoothTreePainting(std::shared_ptr 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); + + // 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(node->m_parent->m_position.toCartesian(), + node->m_position.toCartesian())); + } + + // interior + 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(node->m_parent->m_position.toCartesian(), + node->m_position.toCartesian())); + } +} + +void SmoothTreePainting::paintNodes(renderer::GeometryRenderer& renderer) const { + renderer.setMode(renderer::GeometryRenderer::vertices | renderer::GeometryRenderer::stroke); + 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::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()); + } + } +} + +} // 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 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 m_tree; + Options m_options; +}; + +} // namespace cartocrow::flow_map + +#endif //CARTOCROW_FLOW_MAP_PAINTING diff --git a/cartocrow/flow_map/spiral.cpp b/cartocrow/flow_map/spiral.cpp index 6b6d107f0..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 - #include namespace cartocrow::flow_map { @@ -34,36 +33,58 @@ Spiral::Spiral(const PolarPoint& anchor, const Number& angle) } } -Spiral::Spiral(const PolarPoint& p1, const PolarPoint& p2) { - const PolarPoint& source = p1.r() < p2.r() ? p2 : p1; - const PolarPoint& target = p1.r() < p2.r() ? p1 : p2; +Spiral::Spiral(const PolarPoint& p1, const PolarPoint& p2) + : m_anchor(p1.r() < p2.r() ? p2 : p1), m_angle(alpha(p1, p2)) {} +Number Spiral::alpha(const PolarPoint& p1, const PolarPoint& 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( - "Tried to construct a spiral containing two points equidistant to the root"); + "Cannot compute α for a spiral connecting two points equidistant to the root"); + } + if (source.r() == 0) { + return 0; + } + Number diff_phi = wrapAngle(target.phi() - source.phi(), -M_PI); + return std::atan(diff_phi / -std::log(target.r() / source.r())); +} + +Number Spiral::dAlphaDPhi1(const PolarPoint& p1, const PolarPoint& 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; } + Number phiDiff = wrapAngle(target.phi() - source.phi(), -M_PI); + Number rDiffLog = std::log(target.r() / source.r()); + return rDiffLog / (rDiffLog * rDiffLog + phiDiff * phiDiff); +} + +Number Spiral::dAlphaDPhi2(const PolarPoint& p1, const PolarPoint& p2) { + return -dAlphaDPhi1(p1, p2); +} + +Number Spiral::dAlphaDR1(const PolarPoint& p1, const PolarPoint& 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; + } + Number phiDiff = wrapAngle(target.phi() - source.phi(), -M_PI); + Number rDiffLog = std::log(target.r() / source.r()); + return (-phiDiff / p1.r()) / (rDiffLog * rDiffLog + phiDiff * phiDiff); +} - 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 - // 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) { - m_angle = 0; - } else { - Number diff_phi = wrapAngle(target.phi() - source.phi(), -M_PI); - m_angle = std::atan(diff_phi / -std::log(target.r() / source.r())); +Number Spiral::dAlphaDR2(const PolarPoint& p1, const PolarPoint& 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; } + Number phiDiff = wrapAngle(target.phi() - source.phi(), -M_PI); + Number 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 af8b5f497..9a9d707b1 100644 --- a/cartocrow/flow_map/spiral.h +++ b/cartocrow/flow_map/spiral.h @@ -68,6 +68,74 @@ 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 \f$(r_1, \phi_1)\f$ and \f$(r_2, \phi_2)\f$. + /// + /// 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 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 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 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 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 dAlphaDPhi2(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/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::obstacles() { return m_obstacles; } +std::shared_ptr SpiralTree::root() { + return m_root; +} + void SpiralTree::addPlace(const std::string& name, const Point& position, Number flow) { auto newPlace = std::make_shared(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>& nodes() const; /// Returns a list of the obstacles in this spiral tree. std::vector& obstacles(); + /// Returns the root of this spiral tree. + std::shared_ptr root(); /// Adds a node to the spiral tree. //void addNode(std::shared_ptr node); diff --git a/cartocrow/renderer/geometry_widget.cpp b/cartocrow/renderer/geometry_widget.cpp index 58d956813..86f9abf6f 100644 --- a/cartocrow/renderer/geometry_widget.cpp +++ b/cartocrow/renderer/geometry_widget.cpp @@ -475,8 +475,7 @@ void GeometryWidget::draw(const Point& p) { } void GeometryWidget::draw(const Segment& 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 +545,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); } 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..822d233cf --- /dev/null +++ b/demos/flow_map/optimization_demo/CMakeLists.txt @@ -0,0 +1,20 @@ +set(SOURCES + optimization_demo.cpp + cost_graph.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/cost_graph.cpp b/demos/flow_map/optimization_demo/cost_graph.cpp new file mode 100644 index 000000000..71aa36d6a --- /dev/null +++ b/demos/flow_map/optimization_demo/cost_graph.cpp @@ -0,0 +1,158 @@ +/* +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 . +*/ + +#include "cost_graph.h" + +#include +#include +#include +#include +#include + +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)); + + // axes + painter.drawLine(0, 0, 0, -graphHeight); + painter.drawLine(0, 0, graphWidth, 0); + + 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(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); + painter.drawPath(balancingPath); + painter.setPen(QPen(QColor(240, 90, 40), 2)); + painter.drawPath(straighteningPath); + painter.setBrush(QColor(240, 90, 40)); + + // dot + DataPoint lastData = m_dataPoints[m_dataPoints.size() - 1]; + painter.drawEllipse(QPointF(graphWidth * (m_dataPoints.size() - 1) / + std::max(100.0f, static_cast(m_dataPoints.size())), + -graphHeight * lastData.stackedCost(5) / m_maxCost), + 2, 2); + + // labels + style()->standardPalette(); + 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"); +} + +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 += 1 + i / width()) { + dataPath.lineTo(graphWidth * i / std::max(100.0f, static_cast(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 new file mode 100644 index 000000000..bb0b9e6fd --- /dev/null +++ b/demos/flow_map/optimization_demo/cost_graph.h @@ -0,0 +1,82 @@ +/* +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 . +*/ + +#ifndef CARTOCROW_DEMOS_FLOWMAP_OPTIMIZATIONDEMO_COSTGRAPH_H +#define CARTOCROW_DEMOS_FLOWMAP_OPTIMIZATIONDEMO_COSTGRAPH_H + +#include + +#include +#include + +#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 m_obstacle_cost; + Number m_smoothing_cost; + Number m_angle_restriction_cost; + Number m_balancing_cost; + Number m_straightening_cost; + + inline Number stackedCost(int count) const { + std::array 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 sum = 0; + for (int i = 0; i < count; ++i) { + if (!std::isnan(this->*stackingOrder[i])) { + sum += this->*stackingOrder[i]; + } + } + return sum; + } + }; + /// 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 m_dataPoints; + /// The maximum total cost we've seen so far. + Number m_maxCost = 0; + + QPainterPath createDataPath(int costIndex, int graphWidth, int graphHeight) const; + QPainterPath createRegionBetween(const QPainterPath& first, const QPainterPath& second) const; +}; + +#endif 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..20ffc86e0 --- /dev/null +++ b/demos/flow_map/optimization_demo/optimization_demo.cpp @@ -0,0 +1,172 @@ +/* +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 . +*/ + +#include "optimization_demo.h" + +#include +#include +#include + +#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/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" +#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"); + + for (int i = 0; i < 40; i++) { + m_places.push_back(std::make_shared>( + static_cast(std::rand()) / RAND_MAX * 200 - 100, + static_cast(std::rand()) / RAND_MAX * 200 - 100)); + } + + m_renderer = new GeometryWidget(); + m_renderer->setMaxZoom(10000); + 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"); + 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, [&]() { + for (int i = 0; i < 100; i++) { + m_iterationCount++; + m_smoothTree->optimize(); + updateCostLabel(); + if (m_stopOnNanCheckbox->isChecked() && std::isnan(m_smoothTree->computeCost())) { + m_optimizeButton->setChecked(false); + m_optimizeTimer->stop(); + break; + } + } + 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, [&]() { + m_iterationCount++; + m_smoothTree->optimize(); + m_renderer->update(); + updateCostLabel(); + }); + 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(); + }); + + m_costLabel = new QLabel(); + statusBar()->addWidget(m_costLabel); + + for (auto place : m_places) { + m_renderer->registerEditable(place); + } + connect(m_renderer, &GeometryWidget::edited, [&]() { + recalculate(); + }); + recalculate(); +} + +void OptimizationDemo::recalculate() { + m_renderer->clear(); + m_costGraph->clear(); + m_iterationCount = 0; + auto tree = std::make_shared(Point(0, 0), m_alpha); + for (const auto& place : m_places) { + tree->addPlace("", *place, 1); + } + //tree->addShields(); + + ReachableRegionAlgorithm::ReachableRegion reachableRegion = ReachableRegionAlgorithm(tree).run(); + SpiralTreeObstructedAlgorithm(tree, reachableRegion).run(); + m_smoothTree = std::make_shared(tree); + /*for (auto& node : m_smoothTree->m_nodes) { + if (node->getType() == Node::ConnectionType::kSubdivision) { + node->m_position.setPhi(wrapAngle( + node->m_position.phi() + static_cast(std::rand()) / RAND_MAX * 0.2 - 0.1, + -M_PI)); + } + }*/ + SmoothTreePainting::Options options; + auto painting = std::make_shared(m_smoothTree, options); + m_renderer->addPainting(painting, "Smooth tree"); + + m_renderer->update(); + updateCostLabel(); +} + +void OptimizationDemo::updateCostLabel() { + m_costGraph->addStep( + {0, m_smoothTree->computeSmoothingCost(), m_smoothTree->computeAngleRestrictionCost(), + m_smoothTree->computeBalancingCost(), m_smoothTree->computeStraighteningCost()}); + Number cost = m_smoothTree->computeCost(); + m_costLabel->setText( + QString("Iteration %1 | Cost function: %2") + .arg(m_iterationCount) + .arg(std::isnan(cost) ? "nan" : QString::number(cost))); +} + +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..37ddb8805 --- /dev/null +++ b/demos/flow_map/optimization_demo/optimization_demo.h @@ -0,0 +1,65 @@ +/* +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 . +*/ + +#include +#include +#include +#include +#include +#include + +#include +#include + +#include "cartocrow/core/core.h" +#include "cartocrow/flow_map/smooth_tree.h" +#include "cartocrow/renderer/geometry_widget.h" + +#include "cost_graph.h" + +using namespace cartocrow; +using namespace cartocrow::renderer; + +class OptimizationDemo : public QMainWindow { + Q_OBJECT; + + public: + OptimizationDemo(); + + private: + void recalculate(); + void updateCostLabel(); + Number m_alpha = 25 * M_PI / 180; + + std::vector>> m_places; + std::shared_ptr m_smoothTree; + int m_iterationCount = 0; + + GeometryWidget* m_renderer; + QSlider* m_alphaSlider; + QLabel* m_alphaLabel; + QPushButton* m_optimizeButton; + QCheckBox* m_stopOnNanCheckbox; + QPushButton* m_optimizeOneStepButton; + QTimer* m_optimizeTimer; + QLabel* m_costLabel; + CostGraph* m_costGraph; +}; 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..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 . #include #include +#include #include "cartocrow/core/core.h" #include "cartocrow/core/timer.h" @@ -194,26 +195,12 @@ void SpiralTreeDemo::recalculate() { t.output(); Painting::Options options; - Painting p(nullptr, tree, options); auto painting = std::make_shared(nullptr, tree, options); m_renderer->addPainting(painting, "Spiral tree"); m_renderer->update(); } -/*Point* SpiralTreeDemo::findClosestPoint(Point p, Number radius) { - Point* closest = nullptr; - Number minSquaredDistance = radius * radius; - for (auto& vertex : m_obstacle) { - Number 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;