Line data Source code
1 : /*
2 : * Copyright (c) 2016: G-CSC, Goethe University Frankfurt
3 : * Author: Markus Breit
4 : *
5 : * This file is part of UG4.
6 : *
7 : * UG4 is free software: you can redistribute it and/or modify it under the
8 : * terms of the GNU Lesser General Public License version 3 (as published by the
9 : * Free Software Foundation) with the following additional attribution
10 : * requirements (according to LGPL/GPL v3 §7):
11 : *
12 : * (1) The following notice must be displayed in the Appropriate Legal Notices
13 : * of covered and combined works: "Based on UG4 (www.ug4.org/license)".
14 : *
15 : * (2) The following notice must be displayed at a prominent place in the
16 : * terminal output of covered works: "Based on UG4 (www.ug4.org/license)".
17 : *
18 : * (3) The following bibliography is recommended for citation and must be
19 : * preserved in all covered files:
20 : * "Reiter, S., Vogel, A., Heppner, I., Rupp, M., and Wittum, G. A massively
21 : * parallel geometric multigrid solver on hierarchically distributed grids.
22 : * Computing and visualization in science 16, 4 (2013), 151-164"
23 : * "Vogel, A., Reiter, S., Rupp, M., Nägel, A., and Wittum, G. UG4 -- a novel
24 : * flexible software system for simulating pde based models on high performance
25 : * computers. Computing and visualization in science 16, 4 (2013), 165-179"
26 : *
27 : * This program is distributed in the hope that it will be useful,
28 : * but WITHOUT ANY WARRANTY; without even the implied warranty of
29 : * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
30 : * GNU Lesser General Public License for more details.
31 : *
32 : * Created on: 2016-12-19
33 : */
34 :
35 : #ifndef UG__LIB_GRID__REFINEMENT__PROJECTORS__NEURITE_PROJECTOR_H
36 : #define UG__LIB_GRID__REFINEMENT__PROJECTORS__NEURITE_PROJECTOR_H
37 :
38 : #include "common/types.h"
39 : #include "lib_grid/refinement/projectors/refinement_projector.h"
40 :
41 : #include <boost/serialization/split_member.hpp> // for separate load/save methods
42 :
43 : namespace ug {
44 :
45 : class NeuriteProjector
46 : : public RefinementProjector
47 : {
48 : public:
49 : NeuriteProjector();
50 : NeuriteProjector(SPIGeometry3d geometry);
51 :
52 : virtual ~NeuriteProjector();
53 :
54 : virtual void set_geometry(SPIGeometry3d geometry);
55 :
56 : /// called when a new vertex was created from an old vertex
57 : virtual number new_vertex(Vertex* vrt, Vertex* parent);
58 :
59 : /// called when a new vertex was created from an old edge
60 : virtual number new_vertex(Vertex* vrt, Edge* parent);
61 :
62 : /// called when a new vertex was created from an old face
63 : virtual number new_vertex(Vertex* vrt, Face* parent);
64 :
65 : /// called when a new vertex was created from an old volume
66 : virtual number new_vertex(Vertex* vrt, Volume* parent);
67 :
68 : /// project a vertex to its model position
69 : void project(Vertex* vrt);
70 :
71 : /// spline direction at some grid object
72 : void direction_at_grid_object(vector3& dirOut, GridObject* o) const;
73 :
74 :
75 : protected:
76 : void attach_surf_params();
77 :
78 : public:
79 : struct Section
80 : {
81 0 : Section() : endParam(0) {} // constructor for serialization
82 : Section(number _endParam) // constructor for search with CompareSections
83 0 : : endParam(_endParam) {}
84 :
85 : number endParam;
86 :
87 : number splineParamsX[4];
88 : number splineParamsY[4];
89 : number splineParamsZ[4];
90 : number splineParamsR[4];
91 :
92 : template <class Archive>
93 0 : void serialize(Archive& ar, const unsigned int version)
94 : {
95 0 : ar & endParam;
96 0 : for (size_t i = 0; i < 4; ++i)
97 : {
98 0 : ar & splineParamsX[i];
99 0 : ar & splineParamsY[i];
100 0 : ar & splineParamsZ[i];
101 0 : ar & splineParamsR[i];
102 : }
103 0 : }
104 : };
105 :
106 : struct BranchingPoint;
107 0 : struct BranchingRegion
108 : {
109 0 : BranchingRegion() : t(0.0), bp(SPNULL) {} // constructor for serialization
110 : BranchingRegion(number _t) // constructor for search with CompareBranchingPointEnds
111 0 : : t(_t) {}
112 :
113 : /// the axial parameter where other neurite(s) branch off
114 : number t;
115 :
116 : /// pointer to whole branching point
117 : SmartPtr<BranchingPoint> bp;
118 :
119 : template<class Archive>
120 0 : void save(Archive & ar, const unsigned int version) const
121 : {
122 0 : ar << t;
123 :
124 0 : bool owns_bp = bp->vRegions[0] == this;
125 : ar << owns_bp;
126 0 : if (owns_bp)
127 : ar << *bp;
128 0 : }
129 :
130 : template<class Archive>
131 0 : void load(Archive & ar, const unsigned int version)
132 : {
133 : // invoke serialization of the base class
134 0 : ar >> t;
135 :
136 : bool owns_bp;
137 : ar >> owns_bp;
138 0 : if (owns_bp)
139 : {
140 0 : bp = make_sp(new BranchingPoint());
141 : ar >> *bp;
142 : }
143 0 : }
144 :
145 : BOOST_SERIALIZATION_SPLIT_MEMBER()
146 : };
147 :
148 :
149 0 : struct BranchingPoint
150 : {
151 : std::vector<uint32_t> vNid;
152 : std::vector<BranchingRegion*> vRegions;
153 :
154 : template <class Archive>
155 0 : void serialize(Archive& ar, const unsigned int version)
156 : {
157 : // please note: We could simply use
158 : // ar & vNid
159 : // but this somehow takes up a huge number of template recursions
160 0 : size_t sz = vNid.size(); // this does not hurt in the load-case
161 : ar & sz;
162 0 : vNid.resize(sz); // this does not hurt in the save-case
163 0 : for (size_t i = 0; i < sz; ++i)
164 : ar & vNid[i];
165 0 : }
166 : };
167 :
168 : struct SomaPoint {
169 : vector3 soma;
170 : number radius;
171 : SomaPoint(const vector3& soma, number radius) : soma(soma), radius(radius) {}
172 : };
173 :
174 0 : struct SomaBranchingRegion {
175 : SmartPtr<SomaPoint> somaPt;
176 : number radius;
177 : vector3 center;
178 : number t;
179 : vector3 bp;
180 :
181 : template <class Archive>
182 : void serialize(Archive& ar, const unsigned int version) {
183 0 : ar & radius;
184 0 : ar & t;
185 : }
186 :
187 : template<class Archive>
188 : void load(Archive & ar, const unsigned int version)
189 : {
190 : ar >> radius;
191 : ar >> t;
192 : }
193 :
194 : SomaBranchingRegion(const vector3& center, number radius, number t)
195 : : radius(radius),
196 : t(t),
197 : bp(center)
198 : {}
199 :
200 : SomaBranchingRegion(number t)
201 0 : : radius(-1),
202 0 : t(t)
203 : {}
204 :
205 : SomaBranchingRegion()
206 0 : : radius(-1),
207 0 : t(-1)
208 : {}
209 : };
210 :
211 : /*!
212 : * \brief Mapping of model to surface vertices
213 : */
214 : struct Mapping {
215 : /// TODO: Mapping should store unique indices of 1d vertices from SWC file
216 : vector3 v1; /// start vertex of edge
217 : vector3 v2; /// end vertex of edge
218 : number lambda; /// projection parameter lambda
219 : };
220 :
221 : struct Neurite
222 : {
223 : vector3 refDir;
224 : std::vector<Section> vSec;
225 : std::vector<BranchingRegion> vBR;
226 : std::vector<Section> vSomaSec;
227 : std::vector<SomaBranchingRegion> vSBR;
228 :
229 : template <class Archive>
230 0 : void serialize(Archive& ar, const unsigned int version)
231 : {
232 : /// neurite information
233 0 : ar & refDir;
234 :
235 0 : size_t sz = vSec.size();
236 : ar & sz;
237 0 : vSec.resize(sz);
238 0 : for (size_t i = 0; i < sz; ++i)
239 : ar & vSec[i];
240 :
241 0 : sz = vBR.size();
242 : ar & sz;
243 0 : vBR.resize(sz);
244 0 : for (size_t i = 0; i < sz; ++i)
245 : ar & vBR[i];
246 :
247 : /// soma information
248 0 : sz = vSomaSec.size();
249 : ar & sz;
250 0 : vSomaSec.resize(sz);
251 0 : for (size_t i = 0; i < sz; ++i)
252 : ar & vSomaSec[i];
253 :
254 0 : sz = vSBR.size();
255 : ar & sz;
256 0 : vSBR.resize(sz);
257 0 : for (size_t i = 0; i < sz; ++i) {
258 : ar & vSBR[i];
259 : }
260 0 : }
261 : };
262 :
263 : struct SurfaceParams
264 : {
265 : /**
266 : * @brief Neurite ID a vertex belongs to
267 : *
268 : * Only the low 20 bits are used for ordinary vertices, i.e.,
269 : * vertices that do not belong to two neurites at a branching point.
270 : * In branching points, the high 4 bits encode which of the child
271 : * neurites share the vertex: bit 28 for child 0, bit 29 for child 1,
272 : * bit 30 for child 2 and bit 31 for child 3.
273 : * Bits 20-27 are used to encode the branching region on the parent
274 : * neurite.
275 : * There cannot be more than 256 branching regions per neurite.
276 : * There cannot be more than 4 children per branching point.
277 : */
278 : uint32 neuriteID;
279 : float axial;
280 : float angular;
281 : float radial;
282 : };
283 :
284 :
285 : public:
286 : struct CompareSections
287 : {
288 : bool operator()(const Section& a, const Section& b)
289 0 : {return a.endParam < b.endParam;}
290 : };
291 : struct CompareBranchingRegionEnds
292 : {
293 : bool operator()(const BranchingRegion& a, const BranchingRegion& b)
294 0 : {return a.t < b.t;}
295 : };
296 :
297 : // helper struct for branching point calculations
298 : struct BPProjectionHelper
299 : {
300 : number start;
301 : number end;
302 : std::vector<Section>::const_iterator sec_start;
303 : const SomaBranchingRegion* sr;
304 : };
305 :
306 : void debug_neurites() const;
307 :
308 : struct CompareSomaBranchingRegionsEnd {
309 : bool operator()(const SomaBranchingRegion& a, const SomaBranchingRegion& b) {
310 0 : return a.t < b.t;
311 : }
312 : };
313 :
314 : public:
315 : std::vector<Neurite>& neurites();
316 : std::vector<SomaBranchingRegion>& somata();
317 : const Neurite& neurite(uint32_t nid) const;
318 : Grid::VertexAttachmentAccessor<Attachment<SurfaceParams> >& surface_params_accessor();
319 : const Grid::VertexAttachmentAccessor<Attachment<SurfaceParams> >& surface_params_accessor() const;
320 : const std::vector<std::pair<number, number> >& quadrature_points() const;
321 : void average_pos_from_parent(vector3& posOut, const IVertexGroup* const parent) const;
322 : vector3 position(Vertex* vrt) const;
323 :
324 : number axial_range_around_branching_region
325 : (
326 : uint32_t nid,
327 : size_t brInd,
328 : number numberOfRadii = 5.0
329 : ) const;
330 :
331 : number axial_range_around_soma_region
332 : (
333 : const SomaBranchingRegion& sr,
334 : size_t numRadii,
335 : size_t nid,
336 : Vertex* vrt
337 : ) const;
338 :
339 : bool is_in_axial_range_around_soma_region
340 : (
341 : const SomaBranchingRegion& sr,
342 : size_t nid,
343 : Vertex* vrt
344 : ) const;
345 :
346 : void print_surface_params(const Vertex* const v) const;
347 :
348 : protected:
349 : std::vector<Section>::const_iterator get_section_iterator(uint32_t nid, float t) const;
350 : std::vector<NeuriteProjector::Section>::const_iterator get_soma_section_iterator(uint32_t nid, float t) const;
351 :
352 : void prepare_quadrature();
353 :
354 : void average_params
355 : (
356 : uint32_t& neuriteID,
357 : float& t,
358 : float& angle,
359 : float& radius,
360 : const IVertexGroup* parent
361 : ) const;
362 :
363 : number push_into_place(Vertex* vrt, const IVertexGroup* parent);
364 :
365 : private:
366 : Attachment<SurfaceParams> m_aSurfParams;
367 : Grid::VertexAttachmentAccessor<Attachment<SurfaceParams> > m_aaSurfParams;
368 :
369 : /**
370 : * @brief storage for cubic splines:
371 : * Vector of sorted vectors of sections; each inner vector represents a single neurite.
372 : * Any complete neurite is parameterized by t in [0,1]. Each section in the neurite vector
373 : * consists of the parameter t at which the section ends and the sixteen coefficients
374 : * describing the spline in each dimension (monomial basis {(t_(i+1) - t)^i}_i).
375 : **/
376 : std::vector<Neurite> m_vNeurites; //!< spline information for neurites
377 :
378 : std::vector<std::pair<number, number> > m_qPoints; //!< for quadrature when projecting within branching points
379 :
380 : /*!
381 : * \brief check if an element is inside a sphere
382 : * \param[in] elem
383 : * \param[in] center
384 : * \param[in] radius
385 : * \return \c true if element is inside or on sphere else false
386 : */
387 : template <class TElem>
388 0 : bool IsElementInsideSphere
389 : (
390 : const TElem* elem,
391 : const ug::vector3& center,
392 : const number radius
393 : ) const
394 : {
395 0 : Grid::VertexAttachmentAccessor<APosition> aaPos(this->geometry()->grid(), aPosition);
396 : vector3 c = CalculateCenter(elem, aaPos);
397 : vector3 diff;
398 : VecSubtract(diff, c, center);
399 : VecPow(diff, diff, 2.0);
400 0 : number s = diff.x() + diff.y() + diff.z();
401 0 : return (s <= (sq(radius) + SMALL));
402 : }
403 :
404 : friend class boost::serialization::access;
405 : template<class Archive>
406 0 : void save(Archive & ar, const unsigned int version) const
407 : {
408 : UG_EMPTY_BASE_CLASS_SERIALIZATION(NeuriteProjector, RefinementProjector);
409 :
410 : // only write if data is to be written
411 : if(ArchiveInfo<Archive>::TYPE == AT_DATA)
412 : {
413 : // neurites
414 0 : size_t sz = m_vNeurites.size();
415 : ar << sz;
416 0 : for (size_t i = 0; i < sz; ++i) {
417 : ar << m_vNeurites[i];
418 : }
419 : }
420 :
421 : // do not do anything otherwise
422 : else if (ArchiveInfo<Archive>::TYPE == AT_GUI)
423 : {}
424 0 : }
425 :
426 : template<class Archive>
427 0 : void load(Archive & ar, const unsigned int version)
428 : {
429 : UG_EMPTY_BASE_CLASS_SERIALIZATION(NeuriteProjector, RefinementProjector);
430 :
431 : // neurites
432 : size_t sz;
433 : ar >> sz;
434 0 : m_vNeurites.resize(sz);
435 0 : for (size_t i = 0; i < sz; ++i) {
436 : ar >> m_vNeurites[i];
437 : }
438 :
439 : // reconstruct uninitialized pointers in branching points/ranges
440 : size_t nNeurites = m_vNeurites.size();
441 0 : for (size_t n = 0; n < nNeurites; ++n)
442 : {
443 : Neurite& neurite = m_vNeurites[n];
444 : size_t nBR = neurite.vBR.size();
445 0 : for (size_t b = 0; b < nBR; ++b)
446 : {
447 : BranchingRegion& br = neurite.vBR[b];
448 : SmartPtr<BranchingPoint> bp = br.bp;
449 0 : if (bp.valid() && !bp->vRegions.size())
450 : {
451 : size_t nReg = bp->vNid.size();
452 0 : bp->vRegions.resize(nReg);
453 0 : bp->vRegions[0] = &br;
454 0 : for (size_t r = 1; r < nReg; ++r)
455 : {
456 : UG_ASSERT(m_vNeurites[bp->vNid[r]].vBR.size(),
457 : "No branching region in neurite where there should be one.")
458 0 : bp->vRegions[r] = &m_vNeurites[bp->vNid[r]].vBR[0];
459 0 : bp->vRegions[r]->bp = bp;
460 : }
461 : }
462 : }
463 : }
464 : // debug_neurites();
465 0 : }
466 : BOOST_SERIALIZATION_SPLIT_MEMBER()
467 : };
468 :
469 : // DO NOT CHANGE LINES BELOW! Needed for serialization! //
470 : std::ostream& operator<<(std::ostream& os, const NeuriteProjector::SurfaceParams& surfParams);
471 : std::istream& operator>>(std::istream& in, NeuriteProjector::SurfaceParams& surfParams);
472 : DECLARE_ATTACHMENT_INFO_TRAITS(Attachment<NeuriteProjector::SurfaceParams>, "NeuriteProjectorSurfaceParams");
473 : std::ostream& operator<<(std::ostream& os, const NeuriteProjector::Mapping& mapping);
474 : std::istream& operator>>(std::istream& in, NeuriteProjector::Mapping& mapping);
475 : DECLARE_ATTACHMENT_INFO_TRAITS(Attachment<NeuriteProjector::Mapping>, "Mapping");
476 : } // namespace ug
477 :
478 : #endif // UG__LIB_GRID__REFINEMENT__PROJECTORS__NEURITE_PROJECTOR_H
|
