BioDynaMo  v1.03.61-fe34ced0
adaptor.h
Go to the documentation of this file.
1 // -----------------------------------------------------------------------------
2 //
3 // Copyright (C) 2021 CERN & University of Surrey for the benefit of the
4 // BioDynaMo collaboration. All Rights Reserved.
5 //
6 // Licensed under the Apache License, Version 2.0 (the "License");
7 // you may not use this file except in compliance with the License.
8 //
9 // See the LICENSE file distributed with this work for details.
10 // See the NOTICE file distributed with this work for additional information
11 // regarding copyright ownership.
12 //
13 // -----------------------------------------------------------------------------
14 #ifndef CORE_VISUALIZATION_ROOT_ADAPTOR_H_
15 #define CORE_VISUALIZATION_ROOT_ADAPTOR_H_
16 
17 #include <string>
18 
19 #include <TCanvas.h>
20 #include <TGeoManager.h>
21 #include <TMath.h>
22 #include <TSystem.h>
23 #include <TVector3.h>
24 #include "Math/AxisAngle.h"
25 #include "Math/EulerAngles.h"
26 
27 #include "core/param/param.h"
28 #include "core/resource_manager.h"
29 #include "core/scheduler.h"
30 #include "core/shape.h"
31 #include "core/simulation.h"
32 #include "core/util/log.h"
33 #include "core/util/math.h"
35 
37 using ROOT::Math::AxisAngle;
38 using ROOT::Math::EulerAngles;
39 
40 namespace bdm {
41 
43 class RootAdaptor {
44  public:
46 
48  void Visualize(uint64_t total_steps) {
49  if (!initialized_) {
50  Initialize();
51  initialized_ = true;
52  }
53 
54  auto *param = Simulation::GetActive()->GetParam();
55  if (total_steps % param->visualization_interval != 0) {
56  return;
57  }
58 
59  top_->ClearNodes();
60 
62 
63  rm->ForEachAgent([&](Agent *agent) {
64  auto container = new TGeoVolumeAssembly("A");
65  this->AddBranch(agent, container);
66  top_->AddNode(container, top_->GetNdaughters());
67  });
68 
69  gSystem->ProcessEvents();
70  gGeoManager->Export(outfile_.c_str(), "biodynamo");
71  }
72 
73  void DrawInCanvas(size_t w = 300, size_t h = 300, std::string opt = "") {
74  canvas_->GetListOfPrimitives()->Clear();
75  if (opt.empty()) {
76  canvas_->GetListOfPrimitives()->Add(gGeoManager->GetTopVolume(), "all");
77  } else {
78  opt = "all;" + opt;
79  canvas_->GetListOfPrimitives()->Add(gGeoManager->GetTopVolume(),
80  opt.c_str());
81  }
82  canvas_->SetCanvasSize(static_cast<UInt_t>(w), static_cast<UInt_t>(h));
83  canvas_->Update();
84  canvas_->Draw();
85  }
86 
87  private:
88  void Initialize() {
89  gGeoManager = new TGeoManager("Visualization", "BioDynaMo");
90  canvas_ = new TCanvas("BioDynaMo Canvas", "For ROOT Notebooks ", 300, 300);
91 
93 
94  // Set number of segments for approximating circles in drawing.
95  // Keep it low for better performance.
96  gGeoManager->SetNsegments(15);
97 
98  mat_empty_space_ = new TGeoMaterial("EmptySpace", 0, 0, 0);
99  mat_solid_ = new TGeoMaterial("Solid", .938, 1., 10000.);
100  med_empty_space_ = new TGeoMedium("Empty", 1, mat_empty_space_);
101  med_solid_ = new TGeoMedium("Solid", 1, mat_solid_);
102 
103  // Another way to make top volume for world. In this way it will be
104  // unbounded.
105  top_ = new TGeoVolumeAssembly("WORLD");
106 
107  gGeoManager->SetTopVolume(top_);
108  gGeoManager->SetVisLevel(4);
109 
110  // Number of visualized nodes inside one volume. If you exceed this number,
111  // ROOT will draw nothing.
112  gGeoManager->SetMaxVisNodes(max_viz_nodes_);
113 
114  // Assign the gGeoManager to our canvas to enable resizing of the output
115  // display. The "all" option is necessary to prevent ROOT from hiding
116  // objects (which it does by default for performance reasons).
117  canvas_->GetListOfPrimitives()->Add(gGeoManager->GetTopVolume(), "all");
118 
119  initialized_ = true;
120  }
121 
123  void AddBranch(const Agent *agent, TGeoVolume *container) {
124  switch (agent->GetShape()) {
125  case Shape::kSphere:
126  AddSphere(agent, container);
127  break;
128  case Shape::kCylinder:
129  AddCylinder(agent, container);
130  break;
131  default:
132  Log::Error("RootAdaptor",
133  "Tried to add a shape to the Root visualization that's not "
134  "one of the supported types : ",
135  agent->GetShape());
136  }
137  // to be extended for other object
138  }
139 
141  void AddSphere(const Agent *agent, TGeoVolume *container) {
142  std::string name = agent->GetTypeName() + std::to_string(agent->GetUid());
143  auto radius = agent->GetDiameter() / 2;
144  auto massLocation = agent->GetPosition();
145  auto x = massLocation[0];
146  auto y = massLocation[1];
147  auto z = massLocation[2];
148  auto position = new TGeoTranslation(x, y, z);
149  auto volume = gGeoManager->MakeSphere(name.c_str(), med_solid_, 0, radius);
150  volume->SetLineColor(kBlue);
151  container->AddNode(volume, container->GetNdaughters(), position);
152  }
153 
155  void AddCylinder(const Agent *agent, TGeoVolume *container) {
156  if (auto neurite = dynamic_cast<const NeuriteElement *>(agent)) {
157  std::string name = agent->GetTypeName() + std::to_string(agent->GetUid());
158  auto radius = neurite->GetDiameter() / 2;
159  auto half_length = neurite->GetLength() / 2;
160  auto massLocation = neurite->GetPosition();
161  auto x = massLocation[0];
162  auto y = massLocation[1];
163  auto z = massLocation[2];
164  auto trans = new TGeoTranslation(x, y, z);
165 
166  // The initial orientation of the cylinder
167  TVector3 orig(0, 0, 1);
168 
169  // The vector that we want to orient the cylinder to (symmetry axis
170  // aligns with this vector)
171  Real3 d = neurite->GetSpringAxis();
172  TVector3 dir(d[0], d[1], d[2]);
173  dir = dir.Unit();
174 
175  // Compute the Axis-Angle rotation representation
176  auto dot_product = dir.Dot(orig);
177  auto angle = std::acos(dot_product);
178  // TODO(ahmad): make sure it is `z x dir, and not `dir x z`
179  TVector3 n = dir.Cross(orig);
180  n = n.Unit();
181  auto axis = AxisAngle::AxisVector(n[0], n[1], n[2]);
182  AxisAngle aa(axis, angle);
183  EulerAngles ea(aa);
184  TGeoRotation *rot = new TGeoRotation("rot", Math::ToDegree(ea.Phi()),
185  Math::ToDegree(ea.Theta()),
186  Math::ToDegree(ea.Psi()));
187  TGeoCombiTrans *transrot = new TGeoCombiTrans(*trans, *rot);
188  auto volume = gGeoManager->MakeTube(name.c_str(), med_solid_, 0, radius,
189  half_length);
190  volume->SetLineColor(kBlue);
191  container->AddNode(volume, container->GetNdaughters(), transrot);
192  } else {
193  Log::Error("RootAdaptor", "This is not a cylindrical shaped object!");
194  }
195  }
196 
197  std::string outfile_;
198 
199  TCanvas *canvas_ = nullptr;
200 
202  TGeoVolume *top_;
203 
205  TGeoMaterial *mat_empty_space_;
206  TGeoMaterial *mat_solid_;
207  TGeoMedium *med_empty_space_;
208  TGeoMedium *med_solid_;
209 
211 
214 };
215 
216 } // namespace bdm
217 
218 #endif // CORE_VISUALIZATION_ROOT_ADAPTOR_H_
shape.h
bdm::RootAdaptor::outfile_
std::string outfile_
Definition: adaptor.h:197
bdm
Definition: agent.cc:39
bdm::RootAdaptor::med_solid_
TGeoMedium * med_solid_
Definition: adaptor.h:208
bdm::RootAdaptor::AddSphere
void AddSphere(const Agent *agent, TGeoVolume *container)
Adds a sphere object to the volume.
Definition: adaptor.h:141
bdm::Simulation::GetUniqueName
const std::string & GetUniqueName() const
Returns the name of the simulation.
Definition: simulation.cc:299
bdm::RootAdaptor::max_viz_nodes_
int max_viz_nodes_
Max visualized shapes per volume.
Definition: adaptor.h:213
scheduler.h
bdm::Agent::GetUid
const AgentUid & GetUid() const
Definition: agent.cc:123
bdm::RootAdaptor::AddBranch
void AddBranch(const Agent *agent, TGeoVolume *container)
Recursively adds sphere and its daughters to the container.
Definition: adaptor.h:123
bdm::Agent::GetDiameter
virtual real_t GetDiameter() const =0
bdm::Agent::GetTypeName
virtual const char * GetTypeName() const
Definition: agent.h:134
neurite_element.h
bdm::Agent::GetShape
virtual Shape GetShape() const =0
bdm::RootAdaptor::top_
TGeoVolume * top_
Top volumes for TGeo and TEve (world)
Definition: adaptor.h:202
math.h
bdm::ResourceManager::ForEachAgent
virtual void ForEachAgent(const std::function< void(Agent *)> &function, Functor< bool, Agent * > *filter=nullptr)
Definition: resource_manager.h:206
bdm::Agent
Contains code required by all agents.
Definition: agent.h:79
bdm::RootAdaptor
The class that bridges the simulation code with ROOT Visualization.
Definition: adaptor.h:43
bdm::RootAdaptor::Initialize
void Initialize()
Definition: adaptor.h:88
bdm::RootAdaptor::canvas_
TCanvas * canvas_
Definition: adaptor.h:199
bdm::Simulation::GetResourceManager
ResourceManager * GetResourceManager()
Returns the ResourceManager instance.
Definition: simulation.cc:261
bdm::RootAdaptor::mat_solid_
TGeoMaterial * mat_solid_
Definition: adaptor.h:206
bdm::RootAdaptor::RootAdaptor
RootAdaptor()
Definition: adaptor.h:45
bdm::kSphere
@ kSphere
Definition: shape.h:20
bdm::RootAdaptor::initialized_
bool initialized_
Definition: adaptor.h:210
bdm::Log::Error
static void Error(const std::string &location, const Args &... parts)
Prints error message.
Definition: log.h:79
log.h
bdm::RootAdaptor::mat_empty_space_
TGeoMaterial * mat_empty_space_
Eve materials and medium.
Definition: adaptor.h:205
bdm::Math::ToDegree
static real_t ToDegree(real_t rad)
Definition: math.h:37
simulation.h
bdm::Simulation::GetParam
const Param * GetParam() const
Returns the simulation parameters.
Definition: simulation.cc:271
bdm::MathArray< real_t, 3 >
bdm::Agent::GetPosition
virtual const Real3 & GetPosition() const =0
bdm::kCylinder
@ kCylinder
Definition: shape.h:20
resource_manager.h
param.h
bdm::Simulation::GetActive
static Simulation * GetActive()
This function returns the currently active Simulation simulation.
Definition: simulation.cc:67
bdm::RootAdaptor::DrawInCanvas
void DrawInCanvas(size_t w=300, size_t h=300, std::string opt="")
Definition: adaptor.h:73
bdm::RootAdaptor::Visualize
void Visualize(uint64_t total_steps)
Visualize one timestep.
Definition: adaptor.h:48
bdm::RootAdaptor::med_empty_space_
TGeoMedium * med_empty_space_
Definition: adaptor.h:207
bdm::RootAdaptor::AddCylinder
void AddCylinder(const Agent *agent, TGeoVolume *container)
Adds a cylinder object to the volume.
Definition: adaptor.h:155
bdm::neuroscience::NeuriteElement
Definition: neurite_element.h:56