actsvg/display/datamodel.hpp

0001 // This file is part of the actsvg package.
0002 //
0003 // Copyright (C) 2022 CERN for the benefit of the ACTS project
0004 //
0005 // This Source Code Form is subject to the terms of the Mozilla Public
0006 // License, v. 2.0. If a copy of the MPL was not distributed with this
0007 // file, You can obtain one at http://mozilla.org/MPL/2.0/.
0008
0009 #pragma once
0010
0011 #include <algorithm>
0012 #include <string>
0013 #include <utility>
0014 #include <vector>
0015
0016 #include "actsvg/core.hpp"
0017 #include "actsvg/proto/cluster.hpp"
0018 #include "actsvg/proto/track.hpp"
0019
0020 namespace actsvg {
0021
0022 using namespace defaults;
0023
0024 namespace display {
0025
0026 /** Make a cluster  in 2D
0027  *
0028  * @param grid_ is the grid in which this cluster sits
0029  * @param id_ is the identification tag
0030  * @param cluster_ is the cluster to be drawn
0031  * @param fill_low_ is the low color band
0032  * @param fill_high is the high color band
0033  * @param fill_m_ is the measurement fill color
0034  * @param expand_ is the focus expansion
0035  *
0036  * @return a cluster display and the focussed grid area
0037  *
0038  **/
0039 template <size_t DIM>
0040 std::pair<svg::object, svg::object> cluster(
0041     const svg::object& grid_, const std::string& id_,
0042     const proto::cluster<DIM>& cluster_,
0043     const style::fill& fill_low_ = style::fill{style::color{{255, 255, 0}}},
0044     const style::fill& fill_high_ = style::fill{style::color{{250, 0, 0}}},
0045     const style::fill& fill_m_ = style::fill{style::color{{0, 0, 255}}},
0046     const std::array<unsigned int, 2>& expand_ = {2, 2}) {
0047
0048     svg::object cluster_group = svg::object::create_group(id_);
0049
0050     // Grid indices
0051     unsigned int i_min = std::numeric_limits<unsigned int>::max();
0052     unsigned int i_max = std::numeric_limits<unsigned int>::min();
0053     unsigned int j_min = std::numeric_limits<unsigned int>::max();
0054     unsigned int j_max = std::numeric_limits<unsigned int>::min();
0055
0056     // Cluster low/high data
0057     scalar low_data = std::numeric_limits<scalar>::max();
0058     scalar high_data = std::numeric_limits<scalar>::lowest();
0059
0060     for (const auto& c : cluster_._channels) {
0061         low_data = std::min(low_data, c._data);
0062         high_data = std::max(high_data, c._data);
0063     }
0064     // Cluster data scaling
0065     scalar delta_data = high_data - low_data;
0066
0067     int low_r = fill_low_._fc._rgb[0];
0068     int low_g = fill_low_._fc._rgb[1];
0069     int low_b = fill_low_._fc._rgb[2];
0070     int delta_r = (fill_high_._fc._rgb[0] - low_r);
0071     int delta_g = (fill_high_._fc._rgb[1] - low_g);
0072     int delta_b = (fill_high_._fc._rgb[2] - low_b);
0073
0074     // Measurement, covariance, truth
0075     scalar t_x = 0.;
0076     scalar t_y = 0.;
0077     scalar t_r = 0.;
0078     scalar t_phi = 0.;
0079
0080     scalar m_x = 0.;
0081     scalar m_y = 0.;
0082     scalar m_r = 0.;
0083     scalar m_phi = 0.;
0084
0085     scalar d_r = 0.;
0086     scalar d_phi = 0.;
0087
0088     scalar c_x = 0.;
0089     scalar c_y = 0.;
0090     scalar c_r = 0.;
0091
0092     std::array<scalar, 2> x_range = {std::numeric_limits<scalar>::max(),
0093                                      std::numeric_limits<scalar>::lowest()};
0094     std::array<scalar, 2> y_range = x_range;
0095     std::array<scalar, 2> r_range = x_range;
0096     std::array<scalar, 2> phi_range = x_range;
0097
0098     /// Helper function to generate a fill color
0099     ///
0100     /// @param channel_ the channel in question
0101     auto generate_fill =
0102         [&](const proto::channel<DIM>& channel_) -> style::fill {
0103         // Scale the colors for the filling
0104         scalar rel_delta_data = (channel_._data - low_data) / delta_data;
0105
0106         int ch_r = low_r + int(rel_delta_data * delta_r);
0107         int ch_g = low_g + int(rel_delta_data * delta_g);
0108         int ch_b = low_b + int(rel_delta_data * delta_b);
0109
0110         return style::fill{style::color{{ch_r, ch_g, ch_b}}};
0111     };
0112
0113     if constexpr (DIM == 2) {
0114
0115         if (cluster_._type == proto::cluster<DIM>::e_polar) {
0116             // Polar case, measurement and truth
0117             m_r = cluster_._measurement[DIM - 2];
0118             m_phi = cluster_._measurement[DIM - 1];
0119
0120             d_r = cluster_._variance[DIM - 2];
0121             d_phi = cluster_._variance[DIM - 1];
0122
0123             t_r = cluster_._measurement[DIM - 2];
0124             t_phi = cluster_._measurement[DIM - 1];
0125
0126         } else {
0127             // Cartesian/Fan case
0128             m_x = cluster_._measurement[DIM - 2];
0129             m_y = cluster_._measurement[DIM - 1];
0130             c_x = cluster_._variance[DIM - 2];
0131             c_y = cluster_._variance[DIM - 1];
0132             t_x = cluster_._truth[DIM - 2];
0133             t_y = cluster_._truth[DIM - 1];
0134
0135             // Correlation between -1 (-45 deg)and 1 (45 deg)
0136             scalar corr = cluster_._correlation;
0137             corr = corr < -1._scalar ? -1._scalar
0138                                      : (corr > 1._scalar ? 1._scalar : corr);
0139             c_r = -corr * 45;
0140         }
0141     }
0142
0143     // Loop over the cluster channels and draw them
0144     for (const auto& c : cluster_._channels) {
0145         // Conventional indices
0146         size_t i = 0;
0147         size_t j = 0;
0148         // 2-Dimensional clusters: trivial
0149         if constexpr (DIM == 2) {
0150             i = c._cid[DIM - 2];
0151             j = c._cid[DIM - 1];
0152         }
0153         // 1-Dimensional cluster
0154         if constexpr (DIM == 1) {
0155             if (cluster_._coords[DIM - 1] == proto::cluster<DIM>::e_r or
0156                 cluster_._coords[DIM - 1] == proto::cluster<DIM>::e_x) {
0157                 i = c._cid[DIM - 1];
0158             } else {
0159                 j = c._cid[DIM - 1];
0160             }
0161         }
0162         // Remember for the focus point
0163         i_min = i < i_min ? i : i_min;
0164         i_max = i_max < i ? i : i_max;
0165         j_min = j < j_min ? j : j_min;
0166         j_max = j_max < j ? j : j_max;
0167
0168         // Get the corresponding grid tile (will be a copy already if
0169         // successful)
0170         std::string tile_id =
0171             grid_._id + "_" + std::to_string(i) + "_" + std::to_string(j);
0172         std::optional<svg::object> grid_tile = grid_.find_object(tile_id);
0173
0174         if (grid_tile.has_value()) {
0175             // Get the tile
0176             svg::object tile = grid_tile.value();
0177             // Record the x/y/r/phi area
0178             x_range[0] = std::min(x_range[0], tile._x_range[0]);
0179             x_range[1] = std::max(x_range[1], tile._x_range[1]);
0180             y_range[0] = std::min(y_range[0], tile._y_range[0]);
0181             y_range[1] = std::max(y_range[1], tile._y_range[1]);
0182             r_range[0] = std::min(r_range[0], tile._r_range[0]);
0183             r_range[1] = std::max(r_range[1], tile._r_range[1]);
0184             phi_range[0] = std::min(phi_range[0], tile._phi_range[0]);
0185             phi_range[1] = std::max(phi_range[1], tile._phi_range[1]);
0186
0187             auto points = tile._attribute_map.find("points");
0188             if (points != tile._attribute_map.end()) {
0189                 svg::object cell;
0190                 cell._tag = "polygon";
0191                 cell._id = id_ + "_cell_" + std::to_string(i) + "_" +
0192                            std::to_string(j);
0193                 cell._fill = generate_fill(c);
0194                 cell._stroke = style::stroke();
0195                 cell._attribute_map["points"] = points->second;
0196                 // Add the cell object
0197                 cluster_group.add_object(cell);
0198             }
0199         }
0200     }
0201
0202     // Need to fill the remaining parameters for measurement & truth
0203     if constexpr (DIM == 1) {
0204         if (cluster_._coords[DIM - 1] == proto::cluster<DIM>::e_r) {
0205             m_r = cluster_._measurement[DIM - 1];
0206             t_r = cluster_._truth[DIM - 1];
0207
0208             d_r = cluster_._variance[DIM - 1];
0209             d_phi = std::abs(phi_range[1] - phi_range[0]) /
0210                     static_cast<scalar>(std::sqrt(12.));
0211
0212             m_phi = 0.5_scalar * (phi_range[0] + phi_range[1]);
0213             t_phi = m_phi;
0214         } else if (cluster_._coords[DIM - 1] == proto::cluster<DIM>::e_phi) {
0215             m_phi = cluster_._measurement[DIM - 1];
0216             t_phi = cluster_._truth[DIM - 1];
0217
0218             d_r = std::abs(r_range[1] - r_range[0]) /
0219                   static_cast<scalar>(std::sqrt(12.));
0220             d_phi = cluster_._variance[DIM - 1];
0221
0222             m_r = 0.5_scalar * (r_range[0] + r_range[1]);
0223             t_r = m_r;
0224         } else if (cluster_._coords[DIM - 1] == proto::cluster<DIM>::e_x) {
0225             m_x = cluster_._measurement[DIM - 1];
0226             t_x = cluster_._truth[DIM - 1];
0227             c_x = cluster_._variance[DIM - 1];
0228             c_y = std::abs(y_range[1] - y_range[0]) /
0229                   static_cast<scalar>(std::sqrt(12.));
0230             m_y = 0.5_scalar * (y_range[0] + y_range[1]);
0231             t_y = m_y;
0232
0233             // Correlation between -1 (-45 deg)and 1 (45 deg)
0234             scalar corr = cluster_._correlation;
0235             corr = corr < -1._scalar ? -1._scalar
0236                                      : (corr > 1._scalar ? 1._scalar : corr);
0237             c_r = -corr * 45;
0238
0239         } else {
0240             m_y = cluster_._measurement[DIM - 1];
0241             t_y = cluster_._truth[DIM - 1];
0242             c_y = cluster_._variance[DIM - 1];
0243             c_x = std::abs(x_range[1] - x_range[0]) /
0244                   static_cast<scalar>(std::sqrt(12.));
0245             m_x = 0.5_scalar * (x_range[0] + x_range[1]);
0246             t_x = m_x;
0247
0248             // Correlation between -1 (-45 deg)and 1 (45 deg)
0249             scalar corr = cluster_._correlation;
0250             corr = corr < -1._scalar ? -1._scalar
0251                                      : (corr > 1._scalar ? 1._scalar : corr);
0252             c_r = corr * 45;
0253         }
0254     }
0255
0256     // Finally, if the cluster is polar, transform
0257     if (cluster_._type == proto::cluster<DIM>::e_polar) {
0258
0259         m_x = m_r * std::cos(m_phi);
0260         m_y = m_r * std::sin(m_phi);
0261
0262         scalar cos_phi = std::cos(m_phi);
0263         scalar sin_phi = std::sin(m_phi);
0264
0265         c_x = std::abs(cos_phi * d_r - m_r * sin_phi * d_phi);
0266         c_y = std::abs(sin_phi * d_r + m_r * cos_phi * d_phi);
0267
0268         c_r = m_phi * 180 / pi;
0269
0270         t_x = t_r * std::cos(t_phi);
0271         t_y = t_r * std::sin(t_phi);
0272     }
0273
0274     // Add an error ellipse
0275     style::fill fill_c = fill_m_;
0276     fill_c._fc._opacity = 0.5;
0277     style::stroke stroke_c;
0278     style::transform t_c;
0279     t_c._rot = {c_r, m_x, m_y};
0280     cluster_group.add_object(
0281         draw::ellipse("c", {m_x, m_y}, {c_x, c_y}, fill_c, stroke_c, t_c));
0282
0283     // Add the marker
0284     style::marker marker{"o"};
0285     marker._fill = fill_m_;
0286     cluster_group.add_object(draw::marker("m", {m_x, m_y}, marker));
0287
0288     if (cluster_._mc) {
0289         style::marker truth_marker{"o"};
0290         cluster_group.add_object(draw::marker("t", {t_x, t_y}, truth_marker));
0291     }
0292
0293     // Get the focussed grid
0294     svg::object grid_area = svg::object::create_group(id_ + "_focussed_grid");
0295
0296     // Expand
0297     i_min = (i_min >= expand_[0]) ? i_min - expand_[0] : 0;
0298     i_max += expand_[0];
0299     j_min = (j_min >= expand_[1]) ? j_min - expand_[1] : 0;
0300     j_max += expand_[1];
0301
0302     for (unsigned int i = i_min; i <= i_max; ++i) {
0303         for (unsigned int j = j_min; j < j_max; ++j) {
0304             std::string tile_id =
0305                 grid_._id + "_" + std::to_string(i) + "_" + std::to_string(j);
0306             std::optional<svg::object> grid_tile = grid_.find_object(tile_id);
0307             if (grid_tile.has_value()) {
0308                 grid_area.add_object(grid_tile.value());
0309             }
0310         }
0311     }
0312
0313     return {cluster_group, grid_area};
0314 }
0315
0316 /** Draw a trajectory object
0317  *
0318  * @tparam trajectory_type is the trajectory type
0319  * @tparam view_type is the view type
0320  *
0321  * @param id_ is the identification tag
0322  * @param trajectory_ is the trajectory object
0323  * @param view_ is the view object
0324  * @param bezier_ is the bezier curve flag
0325  */
0326 template <typename trajectory_type, typename view_type>
0327 svg::object trajectory(const std::string& id,
0328                        const trajectory_type& trajectory_,
0329                        const view_type& view_, bool bezier_ = false) {
0330     svg::object trajectory_group = svg::object::create_group(id);
0331
0332     std::vector<point2> points;
0333     points.reserve(trajectory_._path.size());
0334     std::vector<point2> directions;
0335     directions.reserve(trajectory_._path.size());
0336     for (const auto& [pos, dir] : trajectory_._path) {
0337         points.push_back(view_.point(pos));
0338         if (dir.has_value()) {
0339             point2 dir2 = view_.point(dir.value());
0340             // normalize the direction
0341             scalar norm = std::sqrt(dir2[0] * dir2[0] + dir2[1] * dir2[1]);
0342             dir2[0] /= norm;
0343             dir2[1] /= norm;
0344             directions.push_back(view_.point(dir2));
0345         }
0346     }
0347
0348     // Draw the line as polyline when there are no directions
0349     if (directions.empty() or !bezier_) {
0350         trajectory_group.add_object(
0351             draw::polyline(id + "_path", points, trajectory_._path_stroke));
0352     } else {
0353         // Draw the line as a path with directions as bezier curve
0354     }
0355
0356     // Add origin if configured
0357     if (trajectory_._origin_size > 0.) {
0358         trajectory_group.add_object(
0359             draw::circle(id + "_origin", view_.point(trajectory_._origin),
0360                          trajectory_._origin_size, trajectory_._origin_fill,
0361                          trajectory_._origin_stroke));
0362     }
0363
0364     // Add arrow head if configured
0365     if (trajectory_._path_arrow._size > 0. and not directions.empty()) {
0366         point2 arrow_pos = points.back();
0367         point2 arrow_dir = directions.back();
0368         scalar arrow_phi = std::atan2(arrow_dir[1], arrow_dir[0]);
0369         trajectory_group.add_object(draw::marker(
0370             id + "_arrow", arrow_pos, trajectory_._path_arrow, arrow_phi));
0371     }
0372
0373     // Return the trajectory group
0374     return trajectory_group;
0375 }
0376
0377 /** Draw a seed object
0378  *
0379  * @tparam seed_type is the seed type
0380  * @tparam view_type is the view type
0381  *
0382  * @param id_ is the identification tag
0383  * @param seed_ is the seed object
0384  * @param view_ is the view object
0385  *
0386  */
0387 template <typename seed_type, typename view_type>
0388 svg::object seed(const std::string& id_, const seed_type& seed_,
0389                  const view_type& view_) {
0390     svg::object seed_group = svg::object::create_group(id_);
0391
0392     // Draw the tranjectory if it exist
0393     if (seed_._trajectory.has_value()) {
0394         auto traj = trajectory(id_ + "_traj", seed_._trajectory.value(), view_);
0395         seed_group.add_object(traj);
0396     }
0397
0398     // Draw the space points
0399     for (auto [is, sp] : utils::enumerate(seed_._space_points)) {
0400         auto pos = view_.point(sp);
0401         seed_group.add_object(draw::circle(
0402             id_ + "_sp_" + std::to_string(is), pos, seed_._space_point_size,
0403             seed_._space_point_fill, seed_._space_point_stroke));
0404     }
0405
0406     return seed_group;
0407 }
0408
0409 }  // namespace display
0410 }  // namespace actsvg