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 : #include "neurite_projector.h"
36 : #include "common/error.h"
37 : #include "lib_grid/global_attachments.h" // GlobalAttachments
38 : #include "lib_grid/algorithms/debug_util.h" // ElementDebugInfo
39 : #include <cmath>
40 : #include <iterator> // std::distance
41 : #include <boost/lexical_cast.hpp>
42 :
43 : namespace ug {
44 :
45 : static number VecProd(const vector3& a, const vector3& b)
46 : {
47 0 : return a[0]*b[0] + a[1]*b[1] + a[2]*b[2];
48 : }
49 :
50 : static number VecNormSquared(const vector3& a)
51 : {
52 0 : return a[0]*a[0] + a[1]*a[1] + a[2]*a[2];
53 : }
54 :
55 : static number VecNorm(const vector3& a)
56 : {
57 0 : return sqrt(a[0]*a[0] + a[1]*a[1] + a[2]*a[2]);
58 : }
59 :
60 :
61 0 : NeuriteProjector::NeuriteProjector() // : m_quadOrder(80)
62 : {
63 : typedef Attachment<NeuriteProjector::SurfaceParams> NPSurfParam;
64 0 : if (!GlobalAttachments::is_declared("npSurfParams"))
65 0 : GlobalAttachments::declare_attachment<NPSurfParam>("npSurfParams", true);
66 :
67 : // If branching points extend 5r in each direction and we integrate over twice that size
68 : // we add around quadOrder/2 Gaussians with sigma=r over a total length of 20r.
69 : // So we should use about quadOrder=80 to ensure a smooth surface
70 : // that does not look like a pearl necklace.
71 0 : prepare_quadrature();
72 0 : }
73 :
74 :
75 0 : NeuriteProjector::NeuriteProjector(SPIGeometry3d geometry)
76 0 : : RefinementProjector(geometry) // ,m_quadOrder(80)
77 : {
78 : typedef Attachment<NeuriteProjector::SurfaceParams> NPSurfParam;
79 0 : if (!GlobalAttachments::is_declared("npSurfParams"))
80 0 : GlobalAttachments::declare_attachment<NPSurfParam>("npSurfParams", true);
81 :
82 0 : attach_surf_params();
83 0 : prepare_quadrature();
84 0 : }
85 :
86 :
87 0 : NeuriteProjector::~NeuriteProjector() {}
88 :
89 :
90 0 : void NeuriteProjector::set_geometry(SPIGeometry3d geometry)
91 : {
92 : // call base class method
93 0 : RefinementProjector::set_geometry(geometry);
94 0 : attach_surf_params();
95 0 : }
96 :
97 :
98 0 : number NeuriteProjector::new_vertex(Vertex* vrt, Vertex* parent)
99 : {
100 0 : m_aaSurfParams[vrt].neuriteID = m_aaSurfParams[parent].neuriteID;
101 0 : m_aaSurfParams[vrt].axial = m_aaSurfParams[parent].axial;
102 0 : m_aaSurfParams[vrt].angular = m_aaSurfParams[parent].angular;
103 0 : m_aaSurfParams[vrt].radial = m_aaSurfParams[parent].radial;
104 :
105 0 : set_pos(vrt, pos(parent));
106 0 : return 1.0;
107 : }
108 :
109 :
110 0 : number NeuriteProjector::new_vertex(Vertex* vrt, Edge* parent)
111 : {
112 0 : return push_into_place(vrt, parent);
113 : }
114 :
115 :
116 0 : number NeuriteProjector::new_vertex(Vertex* vrt, Face* parent)
117 : {
118 0 : return push_into_place(vrt, parent);
119 : }
120 :
121 :
122 0 : number NeuriteProjector::new_vertex(Vertex* vrt, Volume* parent)
123 : {
124 0 : return push_into_place(vrt, parent);
125 : }
126 :
127 :
128 0 : void NeuriteProjector::project(Vertex* vrt)
129 : {
130 0 : push_into_place(vrt, (IVertexGroup*) NULL);
131 0 : }
132 :
133 :
134 0 : void NeuriteProjector::direction_at_grid_object(vector3& dirOut, GridObject* o) const
135 : {
136 : // treat vertex separately as it is no vertex group
137 0 : if (o->base_object_id() == VERTEX)
138 : {
139 0 : Vertex* v = dynamic_cast<Vertex*>(o);
140 0 : UG_COND_THROW(!v, "Non-vertex with VERTEX base object id.")
141 :
142 : const SurfaceParams& sp = m_aaSurfParams[v];
143 0 : uint32_t nid = sp.neuriteID & ((1 << 20) - 1);
144 0 : float t = sp.axial;
145 :
146 0 : const Section& sec = *get_section_iterator(nid, t);
147 :
148 0 : number te = sec.endParam;
149 : const number* s = &sec.splineParamsX[0];
150 : number& v0 = dirOut[0];
151 0 : v0 = -3.0*s[0]*(te-t) - 2.0*s[1];
152 0 : v0 = v0*(te-t) - s[2];
153 :
154 : s = &sec.splineParamsY[0];
155 : number& v1 = dirOut[1];
156 0 : v1 = -3.0*s[0]*(te-t) - 2.0*s[1];
157 0 : v1 = v1*(te-t) - s[2];
158 :
159 : s = &sec.splineParamsZ[0];
160 : number& v2 = dirOut[2];
161 0 : v2 = -3.0*s[0]*(te-t) - 2.0*s[1];
162 0 : v2 = v2*(te-t) - s[2];
163 :
164 0 : return;
165 : }
166 :
167 0 : IVertexGroup* vrtGrp = dynamic_cast<IVertexGroup*>(o);
168 0 : UG_COND_THROW(!vrtGrp, "Non-vertex element which is not a vertex group.");
169 :
170 : uint32_t nid;
171 : float t, dummy1, dummy2;
172 0 : average_params(nid, t, dummy1, dummy2, vrtGrp);
173 :
174 0 : const Section& sec = *get_section_iterator(nid, t);
175 :
176 0 : number te = sec.endParam;
177 : const number* s = &sec.splineParamsX[0];
178 : number& v0 = dirOut[0];
179 0 : v0 = -3.0*s[0]*(te-t) - 2.0*s[1];
180 0 : v0 = v0*(te-t) - s[2];
181 :
182 : s = &sec.splineParamsY[0];
183 : number& v1 = dirOut[1];
184 0 : v1 = -3.0*s[0]*(te-t) - 2.0*s[1];
185 0 : v1 = v1*(te-t) - s[2];
186 :
187 : s = &sec.splineParamsZ[0];
188 : number& v2 = dirOut[2];
189 0 : v2 = -3.0*s[0]*(te-t) - 2.0*s[1];
190 0 : v2 = v2*(te-t) - s[2];
191 : }
192 :
193 :
194 :
195 0 : void NeuriteProjector::attach_surf_params()
196 : {
197 0 : Grid& grid = this->geometry()->grid();
198 0 : UG_COND_THROW(!GlobalAttachments::is_declared("npSurfParams"),
199 : "GlobalAttachment 'npSurfParams' not declared.");
200 0 : m_aSurfParams = GlobalAttachments::attachment<Attachment<SurfaceParams> >("npSurfParams");
201 :
202 : // make sure surfaceParams attachment is attached
203 0 : UG_COND_THROW(!grid.has_vertex_attachment(m_aSurfParams),
204 : "No surface parameter attachment for neurite projector attached to grid.");
205 0 : m_aaSurfParams.access(grid, m_aSurfParams);
206 0 : }
207 :
208 :
209 0 : void NeuriteProjector::debug_neurites() const
210 : {
211 : UG_LOGN("BPs of neurites in projector:")
212 : // print out branching region data
213 0 : for (size_t n = 0; n < m_vNeurites.size(); ++n)
214 : {
215 : UG_LOGN("Neurite " << n);
216 : const NeuriteProjector::Neurite& neurite = m_vNeurites[n];
217 0 : const std::vector<NeuriteProjector::BranchingRegion> vBR = neurite.vBR;
218 :
219 0 : for (size_t b = 0; b < vBR.size(); ++b)
220 : {
221 0 : UG_LOGN(" BR " << b << ": " << vBR[b].t);
222 : SmartPtr<NeuriteProjector::BranchingPoint> bp = vBR[b].bp;
223 :
224 : UG_LOGN(" associated BP data:")
225 : size_t bpSz = bp->vNid.size();
226 0 : if (bpSz != bp->vRegions.size())
227 : {
228 : UG_LOGN("Size mismatch: vNid " << bpSz << ", vRegions " << bp->vRegions.size());
229 : }
230 : else
231 : {
232 0 : for (size_t i = 0; i < bpSz; ++i)
233 : {
234 0 : UG_LOGN(" " << bp->vNid[i] << " (" << bp->vRegions[i]->t << ")");
235 : }
236 : }
237 : }
238 0 : }
239 0 : }
240 :
241 :
242 0 : std::vector<NeuriteProjector::Neurite>& NeuriteProjector::neurites()
243 : {
244 0 : return m_vNeurites;
245 : }
246 :
247 :
248 0 : const NeuriteProjector::Neurite& NeuriteProjector::neurite(uint32_t nid) const
249 : {
250 0 : UG_COND_THROW(nid >= m_vNeurites.size(),
251 : "Requested neurite index " << nid << " that is not present in the neurite.");
252 0 : return m_vNeurites[nid];
253 : }
254 :
255 :
256 : Grid::VertexAttachmentAccessor<Attachment<NeuriteProjector::SurfaceParams> >&
257 0 : NeuriteProjector::surface_params_accessor()
258 : {
259 0 : return m_aaSurfParams;
260 : }
261 :
262 :
263 : const Grid::VertexAttachmentAccessor<Attachment<NeuriteProjector::SurfaceParams> >&
264 0 : NeuriteProjector::surface_params_accessor() const
265 : {
266 0 : return m_aaSurfParams;
267 : }
268 :
269 :
270 0 : const std::vector<std::pair<number, number> >& NeuriteProjector::quadrature_points() const
271 : {
272 0 : return m_qPoints;
273 : };
274 :
275 :
276 : std::vector<NeuriteProjector::Section>::const_iterator
277 0 : NeuriteProjector::get_section_iterator(uint32_t nid, float t) const
278 : {
279 0 : Section cmpSec(t);
280 0 : const std::vector<Section>& vSections = m_vNeurites[nid].vSec;
281 : std::vector<Section>::const_iterator itSec =
282 0 : std::lower_bound(vSections.begin(), vSections.end(), cmpSec, CompareSections());
283 :
284 0 : UG_COND_THROW(itSec == vSections.end(),
285 : "Could not find section for parameter t = " << t << " in neurite " << nid << ".");
286 :
287 0 : return itSec;
288 : }
289 :
290 : std::vector<NeuriteProjector::Section>::const_iterator
291 0 : NeuriteProjector::get_soma_section_iterator(uint32_t nid, float t) const
292 : {
293 0 : Section cmpSec(t);
294 0 : const std::vector<Section>& vSections = m_vNeurites[nid].vSomaSec;
295 : std::vector<Section>::const_iterator itSec =
296 0 : std::lower_bound(vSections.begin(), vSections.end(), cmpSec, CompareSections());
297 :
298 0 : UG_COND_THROW(itSec == vSections.end(),
299 : "Could not find section for parameter t = " << t << " in neurite " << nid << ".");
300 :
301 0 : return itSec;
302 : }
303 :
304 :
305 0 : static bool cmpQPairs(const std::pair<number, number>& a, const std::pair<number, number>& b)
306 0 : {return a.first < b.first;}
307 :
308 0 : void NeuriteProjector::prepare_quadrature()
309 : {
310 : /*
311 : GaussLegendre gl(m_quadOrder);
312 : size_t nPts = gl.size();
313 : m_qPoints.resize(nPts);
314 : for (size_t i = 0; i < nPts; ++i)
315 : {
316 : m_qPoints[i].first = gl.point(i)[0];
317 : m_qPoints[i].second = gl.weight(i);
318 : }
319 : */
320 0 : m_qPoints.resize(40);
321 : /// points
322 0 : m_qPoints[0].first = 4.8061379124697458903340327798768835266e-1;
323 0 : m_qPoints[1].first = 5.1938620875302541096659672201231164734e-1;
324 0 : m_qPoints[2].first = 4.4195796466237239575827435779598794312e-1;
325 0 : m_qPoints[3].first = 5.5804203533762760424172564220401205688e-1;
326 0 : m_qPoints[4].first = 4.0365120964931445014224157396742505259e-1;
327 0 : m_qPoints[5].first = 5.9634879035068554985775842603257494741e-1;
328 0 : m_qPoints[6].first = 3.6592390749637315942940782759570190829e-1;
329 0 : m_qPoints[7].first = 6.3407609250362684057059217240429809171e-1;
330 0 : m_qPoints[8].first = 3.2900295458712076349625375941040284497e-1;
331 0 : m_qPoints[9].first = 6.7099704541287923650374624058959715503e-1;
332 0 : m_qPoints[10].first = 2.9311039781419749923756012709814315851e-1;
333 0 : m_qPoints[11].first = 7.0688960218580250076243987290185684149e-1;
334 0 : m_qPoints[12].first = 2.584620991569106435457167128775884977e-1;
335 0 : m_qPoints[13].first = 7.415379008430893564542832871224115023e-1;
336 0 : m_qPoints[14].first = 2.2526643745243589896203434723524101488e-1;
337 0 : m_qPoints[15].first = 7.7473356254756410103796565276475898512e-1;
338 0 : m_qPoints[16].first = 1.9372305516600988102369377488465256131e-1;
339 0 : m_qPoints[17].first = 8.0627694483399011897630622511534743869e-1;
340 0 : m_qPoints[18].first = 1.6402165769291022581032274251925294501e-1;
341 0 : m_qPoints[19].first = 8.3597834230708977418967725748074705499e-1;
342 0 : m_qPoints[20].first = 1.3634087240503644835950177412253472572e-1;
343 0 : m_qPoints[21].first = 8.6365912759496355164049822587746527428e-1;
344 0 : m_qPoints[22].first = 1.1084717428674030615251422724675257599e-1;
345 0 : m_qPoints[23].first = 8.8915282571325969384748577275324742401e-1;
346 0 : m_qPoints[24].first = 8.7693884583344168401839884666950613046e-2;
347 0 : m_qPoints[25].first = 9.1230611541665583159816011533304938695e-1;
348 0 : m_qPoints[26].first = 6.7020248393870248089609095822690018215e-2;
349 0 : m_qPoints[27].first = 9.3297975160612975191039090417730998179e-1;
350 0 : m_qPoints[28].first = 4.8950596515562851635873334565753448208e-2;
351 0 : m_qPoints[29].first = 9.5104940348443714836412666543424655179e-1;
352 0 : m_qPoints[30].first = 3.3593595860661733319573916577397141783e-2;
353 0 : m_qPoints[31].first = 9.6640640413933826668042608342260285822e-1;
354 0 : m_qPoints[32].first = 2.1041590393104172097729500273620357453e-2;
355 0 : m_qPoints[33].first = 9.7895840960689582790227049972637964255e-1;
356 0 : m_qPoints[34].first = 1.1370025008112868668314858143548096511e-2;
357 0 : m_qPoints[35].first = 9.8862997499188713133168514185645190349e-1;
358 0 : m_qPoints[36].first = 4.6368806502714967734728238893139225189e-3;
359 0 : m_qPoints[37].first = 9.9536311934972850322652717611068607748e-1;
360 0 : m_qPoints[38].first = 8.8114514472039982518864878970675383211e-4;
361 0 : m_qPoints[39].first = 9.9911885485527960017481135121029324617e-1;
362 :
363 : /// weights
364 0 : m_qPoints[0].second = 3.8752973989212405631861981479163163482e-2;
365 0 : m_qPoints[1].second = 3.8752973989212405631861981479163163482e-2;
366 0 : m_qPoints[2].second = 3.8519909082123982794153767141905124262e-2;
367 0 : m_qPoints[3].second = 3.8519909082123982794153767141905124262e-2;
368 0 : m_qPoints[4].second = 3.8055180950313121185779037961247411506e-2;
369 0 : m_qPoints[5].second = 3.8055180950313121185779037961247411506e-2;
370 0 : m_qPoints[6].second = 3.7361584528984132100094668130662336596e-2;
371 0 : m_qPoints[7].second = 3.7361584528984132100094668130662336596e-2;
372 0 : m_qPoints[8].second = 3.6443291197902029530255341721258917929e-2;
373 0 : m_qPoints[9].second = 3.6443291197902029530255341721258917929e-2;
374 0 : m_qPoints[10].second = 3.530582369564338984774181542764341618e-2;
375 0 : m_qPoints[11].second = 3.530582369564338984774181542764341618e-2;
376 0 : m_qPoints[12].second = 3.3956022907616951912845054115961992992e-2;
377 0 : m_qPoints[13].second = 3.3956022907616951912845054115961992992e-2;
378 0 : m_qPoints[14].second = 3.2402006728300519037277264783376365016e-2;
379 0 : m_qPoints[15].second = 3.2402006728300519037277264783376365016e-2;
380 0 : m_qPoints[16].second = 3.0653121246464469583268998204199297951e-2;
381 0 : m_qPoints[17].second = 3.0653121246464469583268998204199297951e-2;
382 0 : m_qPoints[18].second = 2.87198845496957756833088654552129928e-2;
383 0 : m_qPoints[19].second = 2.87198845496957756833088654552129928e-2;
384 0 : m_qPoints[20].second = 2.6613923491968412177498239886130252278e-2;
385 0 : m_qPoints[21].second = 2.6613923491968412177498239886130252278e-2;
386 0 : m_qPoints[22].second = 2.4347903817536116030717080224073194034e-2;
387 0 : m_qPoints[23].second = 2.4347903817536116030717080224073194034e-2;
388 0 : m_qPoints[24].second = 2.1935454092836635995837343020857747906e-2;
389 0 : m_qPoints[25].second = 2.1935454092836635995837343020857747906e-2;
390 0 : m_qPoints[26].second = 1.9391083987236008819986015645223081127e-2;
391 0 : m_qPoints[27].second = 1.9391083987236008819986015645223081127e-2;
392 0 : m_qPoints[28].second = 1.6730097641273923696339091543205424489e-2;
393 0 : m_qPoints[29].second = 1.6730097641273923696339091543205424489e-2;
394 0 : m_qPoints[30].second = 1.3968503490011700549244578753860538651e-2;
395 0 : m_qPoints[31].second = 1.3968503490011700549244578753860538651e-2;
396 0 : m_qPoints[32].second = 1.1122924597083478630752162092104286604e-2;
397 0 : m_qPoints[33].second = 1.1122924597083478630752162092104286604e-2;
398 0 : m_qPoints[34].second = 8.2105291909539443564317424411819636462e-3;
399 0 : m_qPoints[35].second = 8.2105291909539443564317424411819636462e-3;
400 0 : m_qPoints[36].second = 5.2491422655764068073710855336398261884e-3;
401 0 : m_qPoints[37].second = 5.2491422655764068073710855336398261884e-3;
402 0 : m_qPoints[38].second = 2.2606385492665956292358664390926663639e-3;
403 0 : m_qPoints[39].second = 2.2606385492665956292358664390926663639e-3;
404 :
405 0 : std::sort(m_qPoints.begin(), m_qPoints.end(), cmpQPairs);
406 0 : }
407 :
408 0 : static void compute_first_section_end_points
409 : (
410 : vector3& pos0Out,
411 : vector3& pos1Out,
412 : std::vector<NeuriteProjector::Section>::const_iterator secIt
413 : )
414 : {
415 0 : number te = secIt->endParam;
416 : const number* s = &secIt->splineParamsX[0];
417 : number& p0 = pos0Out[0];
418 0 : p0 = s[0]*te + s[1];
419 0 : p0 = p0*te + s[2];
420 0 : p0 = p0*te + s[3];
421 0 : pos1Out[0] = s[3];
422 :
423 : s = &secIt->splineParamsY[0];
424 : number& p1 = pos0Out[1];
425 0 : p1 = s[0]*te + s[1];
426 0 : p1 = p1*te + s[2];
427 0 : p1 = p1*te + s[3];
428 0 : pos1Out[1] = s[3];
429 :
430 : s = &secIt->splineParamsZ[0];
431 : number& p2 = pos0Out[2];
432 0 : p2 = s[0]*te + s[1];
433 0 : p2 = p2*te + s[2];
434 0 : p2 = p2*te + s[3];
435 0 : pos1Out[2] = s[3];
436 0 : }
437 :
438 :
439 0 : static void compute_position_and_velocity_in_section
440 : (
441 : vector3& posAxOut,
442 : vector3& velOut,
443 : number& radiusOut,
444 : std::vector<NeuriteProjector::Section>::const_iterator secIt,
445 : float t
446 : )
447 : {
448 0 : number te = secIt->endParam;
449 : const number* s = &secIt->splineParamsX[0];
450 : number& p0 = posAxOut[0];
451 : number& v0 = velOut[0];
452 0 : p0 = s[0]*(te-t) + s[1]; v0 = -3.0*s[0]*(te-t) - 2.0*s[1];
453 0 : p0 = p0*(te-t) + s[2]; v0 = v0*(te-t) - s[2];
454 0 : p0 = p0*(te-t) + s[3];
455 :
456 : s = &secIt->splineParamsY[0];
457 : number& p1 = posAxOut[1];
458 : number& v1 = velOut[1];
459 0 : p1 = s[0]*(te-t) + s[1]; v1 = -3.0*s[0]*(te-t) - 2.0*s[1];
460 0 : p1 = p1*(te-t) + s[2]; v1 = v1*(te-t) - s[2];
461 0 : p1 = p1*(te-t) + s[3];
462 :
463 : s = &secIt->splineParamsZ[0];
464 : number& p2 = posAxOut[2];
465 : number& v2 = velOut[2];
466 0 : p2 = s[0]*(te-t) + s[1]; v2 = -3.0*s[0]*(te-t) - 2.0*s[1];
467 0 : p2 = p2*(te-t) + s[2]; v2 = v2*(te-t) - s[2];
468 0 : p2 = p2*(te-t) + s[3];
469 :
470 : s = &secIt->splineParamsR[0];
471 0 : radiusOut = s[0]*(te-t) + s[1];
472 0 : radiusOut = radiusOut*(te-t) + s[2];
473 0 : radiusOut = radiusOut*(te-t) + s[3];
474 0 : }
475 :
476 :
477 : static void compute_radius_in_section
478 : (
479 : number& radiusOut,
480 : std::vector<NeuriteProjector::Section>::const_iterator secIt,
481 : float t
482 : )
483 : {
484 0 : number te = secIt->endParam;
485 : const number* s = &secIt->splineParamsR[0];
486 0 : radiusOut = s[0]*(te-t) + s[1];
487 0 : radiusOut = radiusOut*(te-t) + s[2];
488 0 : radiusOut = radiusOut*(te-t) + s[3];
489 : }
490 :
491 :
492 0 : static inline void compute_ONB
493 : (
494 : vector3& firstOut,
495 : vector3& secondOut,
496 : vector3& thirdOut,
497 : const vector3& firstIn,
498 : const vector3& refDir
499 : )
500 : {
501 0 : VecNormalize(firstOut, firstIn);
502 : number fac = VecProd(refDir, firstOut);
503 0 : VecScaleAdd(secondOut, 1.0, refDir, -fac, firstOut);
504 0 : VecNormalize(secondOut, secondOut);
505 :
506 0 : VecCross(thirdOut, firstOut, secondOut);
507 0 : }
508 :
509 :
510 : // computes the polar angle phi of a vector
511 : // given a rotational axis direction and two additional
512 : // directions, together being a orthoNORMAL basis
513 0 : static inline number compute_angle
514 : (
515 : const vector3& axis,
516 : const vector3& secondAxis,
517 : const vector3& thirdAxis,
518 : const vector3& posFromAxis
519 : )
520 : {
521 : // eliminate component in axis direction
522 : vector3 posFromAxisCopy;
523 0 : VecScaleAdd(posFromAxisCopy, 1.0, posFromAxis, -VecProd(posFromAxis, axis), axis);
524 :
525 : // compute coord vector (components 2 and 3) w.r.t. ONB
526 : vector2 relCoord;
527 0 : relCoord[0] = VecProd(posFromAxis, secondAxis);
528 0 : VecScaleAdd(posFromAxisCopy, 1.0, posFromAxisCopy, -relCoord[0], secondAxis);
529 0 : relCoord[1] = VecProd(posFromAxisCopy, thirdAxis);
530 0 : VecNormalize(relCoord, relCoord);
531 :
532 : // get angle from coord vector
533 : number angle;
534 0 : if (fabs(relCoord[0]) < 1e-8)
535 0 : angle = relCoord[1] < 0 ? 1.5*PI : 0.5*PI;
536 : else
537 0 : angle = relCoord[0] < 0 ? PI - atan(-relCoord[1]/relCoord[0]) : atan(relCoord[1]/relCoord[0]);
538 :
539 : // angle should be in [0,2*PI)
540 0 : if (angle < 0) angle += 2.0*PI;
541 :
542 0 : return angle;
543 : }
544 :
545 :
546 0 : void NeuriteProjector::average_params
547 : (
548 : uint32_t& neuriteID,
549 : float& t,
550 : float& angle,
551 : float& radius,
552 : const IVertexGroup* parent
553 : ) const
554 : {
555 : // find neurite ID parent belongs to
556 0 : size_t nVrt = parent->num_vertices();
557 : //std::vector<std::vector<std::pair<uint32_t, uint32_t> > > candidates(nVrt);
558 : std::map<uint32_t, size_t> candidates;
559 : std::map<uint32_t, uint32_t> branchInfo;
560 0 : for (size_t i = 0; i < nVrt; ++i)
561 : {
562 0 : const SurfaceParams& surfParams = m_aaSurfParams[parent->vertex(i)];
563 0 : uint32_t thisNID = surfParams.neuriteID & ((1 << 20) - 1);
564 0 : ++candidates[thisNID];
565 :
566 0 : if (surfParams.neuriteID > (1 << 20)) // vrt also belongs to branch(es)
567 : {
568 0 : const uint32_t br = (surfParams.neuriteID & (255 << 20)) >> 20;
569 0 : const std::vector<uint32_t>& vNID = m_vNeurites[thisNID].vBR[br].bp->vNid;
570 0 : const uint32_t children = surfParams.neuriteID >> 28;
571 0 : for (int j = 0; j < 4; ++j)
572 : {
573 0 : if (children & (1 << j))
574 : {
575 0 : ++candidates[vNID[j+1]];
576 0 : branchInfo[vNID[j+1]] = (br << 20) + (1 << (28+j));
577 : }
578 : }
579 : }
580 : }
581 :
582 : // vertex gets all NIDs that have maximal counts
583 : std::vector<uint32_t> maxCands;
584 : size_t maxCnt = 0;
585 :
586 : std::map<uint32_t, size_t>::const_iterator it = candidates.begin();
587 : std::map<uint32_t, size_t>::const_iterator itEnd = candidates.end();
588 0 : for (; it != itEnd; ++it)
589 0 : if (it->second > maxCnt)
590 : maxCnt = it->second;
591 0 : for (it = candidates.begin(); it != itEnd; ++it)
592 0 : if (it->second == maxCnt)
593 0 : maxCands.push_back(it->first);
594 :
595 : const size_t nMaxCands = maxCands.size();
596 0 : if (nMaxCands == 1)
597 0 : neuriteID = maxCands[0];
598 : else
599 : {
600 0 : neuriteID = *std::min_element(maxCands.begin(), maxCands.end());
601 0 : for (size_t i = 0; i < nMaxCands; ++i)
602 : {
603 : std::map<uint32_t, uint32_t>::const_iterator it2 = branchInfo.find(maxCands[i]);
604 0 : if (it2 != branchInfo.end())
605 0 : neuriteID = neuriteID | it2->second;
606 : }
607 : }
608 :
609 :
610 : // special case:
611 : // when refining the initial BP volume the center vertex needs to belong
612 : // to all branches and the parent
613 0 : if (branchInfo.size() && nVrt == 8)
614 : {
615 : // all must have the same radius and parent
616 : // and all children need to have count 4
617 0 : const number firstRad = m_aaSurfParams[parent->vertex(0)].radial;
618 0 : const uint32_t firstParent = m_aaSurfParams[parent->vertex(0)].neuriteID & ((1 << 20) - 1);
619 : bool allSame = true;
620 0 : for (size_t i = 1; i < 8; ++i)
621 : {
622 0 : if (m_aaSurfParams[parent->vertex(i)].radial != firstRad
623 0 : || (m_aaSurfParams[parent->vertex(i)].neuriteID & ((1 << 20) - 1)) != firstParent)
624 : {
625 : allSame = false;
626 : break;
627 : }
628 : }
629 :
630 0 : if (allSame)
631 : {
632 : bool allChildCounts4 = true;
633 :
634 : std::map<uint32_t, uint32_t>::const_iterator it2 = branchInfo.begin();
635 : std::map<uint32_t, uint32_t>::const_iterator it2End = branchInfo.end();
636 0 : for (; it2 != it2End; ++it2)
637 : {
638 0 : if (candidates[it2->first] != 4)
639 : {
640 : allChildCounts4 = false;
641 : break;
642 : }
643 : }
644 :
645 0 : if (allChildCounts4)
646 : {
647 0 : neuriteID = candidates.begin()->first;
648 0 : for (it2 = branchInfo.begin(); it2 != it2End; ++it2)
649 0 : neuriteID = neuriteID | it2->second;
650 : }
651 : }
652 : }
653 :
654 : // special case:
655 : // when refining the initial connecting quadrilateral face of a branch
656 : // the center vertex needs to belong exclusively to the branch, not to the parent
657 0 : if (nMaxCands == 2 && nVrt == 4 && maxCnt == 4)
658 0 : neuriteID = std::max(maxCands[0], maxCands[1]);
659 :
660 0 : const uint32_t plainNID = neuriteID & ((1 << 20) - 1);
661 :
662 :
663 : // average t, phi, r
664 0 : t = 0.0;
665 0 : float minRad = 1.0;
666 0 : float maxRad = 0.0;
667 : vector2 v(0.0);
668 : size_t nForeignParents = 0;
669 0 : std::vector<number> foreignParentRad(8, 0.0);
670 : vector3 foreignParentPosCenter = 0.0;
671 : number foreignParentsCenterRadial = 0.0;
672 :
673 : size_t nForeignChildren = 0;
674 : vector3 foreignChildrenPosCenter = 0.0;
675 :
676 0 : for (size_t i = 0; i < nVrt; ++i)
677 : {
678 0 : const SurfaceParams& surfParams = m_aaSurfParams[parent->vertex(i)];
679 0 : uint32_t nid = surfParams.neuriteID & ((1 << 20) - 1);
680 0 : minRad = std::min(minRad, surfParams.radial);
681 0 : maxRad = std::max(maxRad, surfParams.radial);
682 :
683 : // parent vertex belongs to our branch
684 0 : if (nid == plainNID)
685 : {
686 0 : t += surfParams.axial;
687 0 : if (surfParams.radial > 0 && surfParams.axial < 2.0)
688 : {
689 0 : const float phi = surfParams.angular;
690 0 : VecAdd(v, v, vector2(surfParams.radial*cos(phi), surfParams.radial*sin(phi)));
691 : }
692 : }
693 :
694 : // parent vertex belongs to parent branch
695 0 : else if (nid < plainNID)
696 : {
697 0 : foreignParentPosCenter += this->pos(parent->vertex(i));
698 0 : foreignParentRad[nForeignParents] = surfParams.radial;
699 0 : ++nForeignParents;
700 : }
701 :
702 : // parent vertex belongs to child branch
703 : else
704 : {
705 0 : foreignChildrenPosCenter += this->pos(parent->vertex(i));
706 0 : ++nForeignChildren;
707 : }
708 : }
709 :
710 : // for parent vertices, calculate the t and phi params of their center
711 : // w.r.t. branching neurite
712 0 : if (nForeignParents)
713 : {
714 0 : foreignParentPosCenter /= nForeignParents;
715 :
716 0 : std::vector<Section>::const_iterator secIt = m_vNeurites[plainNID].vSec.begin();
717 : vector3 pos0, pos1;
718 0 : compute_first_section_end_points(pos0, pos1, secIt);
719 : VecSubtract(pos1, pos1, pos0);
720 : VecSubtract(pos0, foreignParentPosCenter, pos0);
721 0 : number axPos = VecProd(pos0, pos1) / VecNormSquared(pos1) * secIt->endParam;
722 :
723 : vector3 vel;
724 : number surfRadius;
725 0 : compute_position_and_velocity_in_section(pos1, vel, surfRadius, secIt, axPos);
726 0 : const vector3& refDir = m_vNeurites[plainNID].refDir;
727 : vector3 posRefDir, thirdDir;
728 0 : compute_ONB(vel, posRefDir, thirdDir, vel, refDir);
729 : vector3 posFromAxis;
730 : VecSubtract(posFromAxis, foreignParentPosCenter, pos1);
731 0 : number phi = compute_angle(vel, posRefDir, thirdDir, posFromAxis);
732 :
733 0 : t += nForeignParents*axPos;
734 0 : foreignParentsCenterRadial = VecNorm(posFromAxis) / surfRadius;
735 0 : VecScaleAdd(v, 1.0, v, nForeignParents, vector2(foreignParentsCenterRadial*cos(phi), foreignParentsCenterRadial*sin(phi)));
736 : }
737 :
738 : // for child vertices, calculate the t and phi params of their center
739 : // w.r.t. parent neurite TODO: ATM this is not done as it is not overly important.
740 0 : if (nForeignChildren)
741 : {
742 0 : t *= ((number) nVrt) / (nVrt - nForeignChildren);
743 : VecScale(v, v, ((number) nVrt) / (nVrt - nForeignChildren));
744 : }
745 :
746 :
747 0 : t /= nVrt;
748 0 : radius = 0.5*(maxRad + minRad);
749 :
750 : // special case initial face of a branch:
751 : // move frontier between branch and parent into branch
752 0 : if (nForeignParents == 4 && nVrt == 4)
753 0 : t *= sqrt(2.0);
754 :
755 : // special case tip of neurite
756 0 : if (fabs(t-2.0) < 1e-8)
757 : return;
758 :
759 0 : if (VecTwoNorm(v) / nVrt < 0.4 * radius && (nVrt == 4 || nVrt == 8))
760 : {
761 0 : if (VecTwoNorm(v) / nVrt > 0.2 * radius)
762 : {
763 0 : UG_LOGN("HERE: radius " << radius << " -> 0 (real: " << VecTwoNorm(v) / nVrt << ")");
764 : //UG_LOGN(" plainNID " << plainNID);
765 : UG_LOGN(" nVrt " << nVrt);
766 : UG_LOGN(" parent vertex positions");
767 0 : for (size_t i = 0; i < nVrt; ++i)
768 0 : UG_LOGN(" " << this->position(parent->vertex(i)));
769 : }
770 : angle = 0.0;
771 0 : radius = 0.0;
772 : }
773 :
774 0 : if (fabs(v[0]) < 1e-4)
775 : {
776 0 : if (fabs(v[1]) < 1e-4) // do not alter radius in tip
777 : {
778 0 : angle = 0.0; // center angle is 0.0 by default
779 0 : radius = 0.0;
780 : }
781 : else
782 0 : angle = v[1] < 0 ? 1.5*PI : 0.5*PI;
783 : }
784 : else
785 0 : angle = v[0] < 0 ? PI - atan(-v[1]/v[0]) : atan(v[1]/v[0]);
786 :
787 0 : if (angle < 0)
788 0 : angle += 2.0*PI;
789 :
790 0 : }
791 :
792 :
793 0 : void NeuriteProjector::average_pos_from_parent(vector3& posOut, const IVertexGroup* const parent) const
794 : {
795 : posOut = 0.0;
796 0 : size_t nVrt = parent->num_vertices();
797 0 : for (size_t i = 0; i < nVrt; ++i)
798 0 : posOut += this->pos(parent->vertex(i));
799 0 : posOut /= (number) nVrt;
800 0 : }
801 :
802 :
803 : /*
804 : void NeuriteProjector::average_pos_from_parent_weighted(vector3& posOut, const IVertexGroup* parent) const
805 : {
806 : MultiGrid* mg = dynamic_cast<MultiGrid*>(&this->geometry()->grid());
807 : UG_COND_THROW(!mg, "Underlying grid is not a multigrid.");
808 : int lv = mg->get_level<Vertex>(parent->vertex(0));
809 : if (lv == 0)
810 : lv = 1;
811 :
812 : posOut = 0.0;
813 : size_t nVrt = parent->num_vertices();
814 : std::vector<int> vSheathID(nVrt);
815 : vSheathID[0] = 0;
816 : int sheathID[2] = {std::round(m_aaSurfParams[parent->vertex(0)].radial*(1 << (lv-1))), -1};
817 : int count[2] = {1, 0};
818 : for (size_t i = 1; i < nVrt; ++i)
819 : {
820 : int thisSheathID = std::round(m_aaSurfParams[parent->vertex(i)].radial*(1 << (lv-1)));
821 : if (thisSheathID != sheathID[0] && sheathID[1] == -1)
822 : sheathID[1] = thisSheathID;
823 : if (thisSheathID == sheathID[0])
824 : {
825 : ++count[0];
826 : vSheathID[i] = 0;
827 : }
828 : else if (thisSheathID == sheathID[1])
829 : {
830 : ++count[1];
831 : vSheathID[i] = 1;
832 : }
833 : else
834 : UG_THROW("An element must only be part of at most 2 sheaths.");
835 : }
836 :
837 : if (!count[1])
838 : count[1] = 1;
839 :
840 : for (size_t i = 0; i < nVrt; ++i)
841 : VecScaleAppend(posOut, count[1-vSheathID[i]], this->pos(parent->vertex(i)));
842 :
843 : posOut /= 2*count[0]*count[1];
844 : }
845 : */
846 :
847 :
848 0 : vector3 NeuriteProjector::position(Vertex* vrt) const
849 : {
850 0 : return this->pos(vrt);
851 : }
852 :
853 0 : number NeuriteProjector::axial_range_around_branching_region
854 : (
855 : uint32_t nid,
856 : size_t brInd,
857 : number numberOfRadii
858 : ) const
859 : {
860 0 : const Neurite& neurite = m_vNeurites[nid];
861 : const BranchingRegion& br = neurite.vBR[brInd];
862 :
863 : std::vector<Section>::const_iterator secIt;
864 0 : try {secIt = get_section_iterator(nid, br.t);}
865 0 : UG_CATCH_THROW("Could not get section iterator to branching region " << brInd << " at t = " << br.t
866 : << " of neurite " << nid << " (which has " << neurite.vBR.size() << " branching regions).");
867 :
868 : vector3 secStartPos;
869 : vector3 secEndPos;
870 0 : const number te = secIt->endParam;
871 : number ts = 0.0;
872 0 : if (secIt == m_vNeurites[nid].vSec.begin())
873 0 : compute_first_section_end_points(secStartPos, secEndPos, secIt);
874 : else
875 : {
876 0 : secEndPos[0] = secIt->splineParamsX[3];
877 0 : secEndPos[1] = secIt->splineParamsY[3];
878 0 : secEndPos[2] = secIt->splineParamsZ[3];
879 : --secIt;
880 0 : ts = secIt->endParam;
881 0 : secStartPos[0] = secIt->splineParamsX[3];
882 0 : secStartPos[1] = secIt->splineParamsY[3];
883 0 : secStartPos[2] = secIt->splineParamsZ[3];
884 : ++secIt;
885 : }
886 0 : const number neuriteLengthApprox = VecDistance(secEndPos, secStartPos) / (te - ts);
887 :
888 : const number* s = &secIt->splineParamsR[0];
889 0 : number radius = s[0]*(te-br.t) + s[1];
890 0 : radius = radius*(te-br.t) + s[2];
891 0 : radius = radius*(te-br.t) + s[3];
892 :
893 0 : return numberOfRadii*radius / neuriteLengthApprox;
894 : }
895 :
896 :
897 0 : void NeuriteProjector::print_surface_params(const Vertex* const v) const
898 : {
899 : const SurfaceParams& sp = m_aaSurfParams[v];
900 0 : UG_LOGN("neuriteID: " << sp.neuriteID);
901 0 : UG_LOGN("axial: " << sp.axial);
902 0 : UG_LOGN("angular: " << sp.angular);
903 0 : UG_LOGN("radial: " << sp.radial);
904 0 : }
905 :
906 :
907 0 : static void bp_defect_and_gradient
908 : (
909 : number& defectOut,
910 : vector3& gradientOut,
911 : const std::vector<NeuriteProjector::BPProjectionHelper>& bpList,
912 : const vector3& x,
913 : number rad,
914 : const vector3& constAngleSurfNormal,
915 : const NeuriteProjector* np
916 : )
917 : {
918 : // get integration rule points and weights
919 0 : const std::vector<std::pair<number, number> >& qPoints = np->quadrature_points();
920 : size_t nPts = qPoints.size();
921 :
922 : // evaluate integrand at points
923 0 : defectOut = 0.0;
924 : gradientOut = 0.0;
925 : size_t nParts = bpList.size();
926 0 : for (size_t i = 0; i < nParts; ++i)
927 : {
928 : const NeuriteProjector::BPProjectionHelper& helper = bpList[i];
929 : const number& start = helper.start;
930 : const number& end = helper.end;
931 0 : std::vector<NeuriteProjector::Section>::const_iterator secIt = helper.sec_start;
932 :
933 : // iterate IPs
934 0 : for (size_t j = 0; j < nPts; ++j)
935 : {
936 : // transform point coords to current context
937 0 : float t = (float) (start + qPoints[j].first*(end-start));
938 0 : while (secIt->endParam < t) ++secIt; // this should be safe as last section must end with 1.0
939 :
940 : vector3 posAx, vel, blub;
941 : number radius;
942 0 : compute_position_and_velocity_in_section(posAx, vel, radius, secIt, t);
943 :
944 0 : number radInv = 1.0 / (rad*radius);
945 : VecSubtract(posAx, posAx, x);
946 : number posNormSq = VecNormSquared(posAx);
947 0 : number velNorm = sqrt(VecNormSquared(vel));
948 0 : number integrandVal = radInv*exp(-posNormSq*radInv*radInv) * velNorm;
949 0 : number gradIntegrandVal = 2.0*radInv*radInv * integrandVal;
950 : VecScale(blub, posAx, gradIntegrandVal);
951 :
952 0 : number w = qPoints[j].second;
953 0 : defectOut += integrandVal * w * (end-start);
954 0 : VecScaleAdd(gradientOut, 1.0, gradientOut, w * (end-start), blub);
955 : #if 0
956 : UG_LOGN(" Ival: " << integrandVal << ", iExp: " << -posNormSq*radInv*radInv
957 : << ", velNorm: " << velNorm << ", weight: " << w
958 : << ", t: " << t << ", intvl: " << end-start
959 : << ", posAx: " << posAx << "x: " << x);
960 :
961 : UG_LOGN(" gradientOut: " << gradientOut);
962 : #endif
963 : }
964 : }
965 :
966 0 : defectOut -= sqrt(PI)*exp(-1.0);
967 :
968 : // project gradient to surface of constant angle
969 0 : VecScaleAppend(gradientOut, -VecDot(constAngleSurfNormal, gradientOut), constAngleSurfNormal);
970 0 : }
971 :
972 : #if 0
973 : static void pos_on_surface_soma
974 : (
975 : vector3& posOut,
976 : const NeuriteProjector::Neurite& neurite,
977 : const NeuriteProjector* np,
978 : const IVertexGroup* parent
979 : )
980 : {
981 : if (parent) np->average_pos_from_parent(posOut, parent);
982 : }
983 :
984 : static void bp_defect_and_gradient_soma
985 : (
986 : number& defectOut,
987 : vector3& gradientOut,
988 : const std::vector<NeuriteProjector::BPProjectionHelper>& bpList,
989 : const vector3& x,
990 : const NeuriteProjector* np,
991 : const vector3& constAngleSurfNormal
992 : )
993 : {
994 : // get integration rule points and weights
995 : const std::vector<std::pair<number, number> >& qPoints = np->quadrature_points();
996 : size_t nPts = qPoints.size();
997 :
998 : // evaluate integrand at points
999 : defectOut = 0.0;
1000 : gradientOut = 0.0;
1001 : const NeuriteProjector::BPProjectionHelper& helper = bpList[0];
1002 : const number& start = helper.start;
1003 : const number& end = helper.end;
1004 : const number& somaRadius = helper.sr->somaPt.get()->radius;
1005 : const vector3& posSoma = helper.sr->somaPt.get()->soma;
1006 : std::vector<NeuriteProjector::Section>::const_iterator secIt = helper.sec_start;
1007 :
1008 : // iterate IPs
1009 : for (size_t j = 0; j < nPts; ++j)
1010 : {
1011 : /// NEURITE part
1012 : // transform point coords to current context
1013 : float t = (float) (start + qPoints[j].first*(end-start));
1014 :
1015 : vector3 posAx, vel, blub;
1016 : number radius;
1017 : /// First section is always BP to soma PM or ER surface of inner sphere
1018 : // /t/=10000;
1019 : /// if (t < 0) { t = -t; }
1020 : UG_LOGN("Section information... t: " << t << ", secIt->endParam(): " << secIt->endParam);
1021 : compute_position_and_velocity_in_section(posAx, vel, radius, secIt, t);
1022 : UG_LOGN("position in neurite: " << posAx);
1023 :
1024 : /// FIXME: add back rad in denominator (rad*radius) from current vertex's aaSurfParams
1025 : UG_LOGN("x: " << x);
1026 : number radInv = 1.0 / (radius*radius);
1027 : VecSubtract(posAx, posAx, x);
1028 : number posNormSq = VecNormSquared(posAx);
1029 : number velNorm = sqrt(VecNormSquared(vel));
1030 : number integrandVal = radInv*exp(-posNormSq*radInv*radInv) * velNorm;
1031 : number gradIntegrandVal = 2.0*radInv*radInv * integrandVal;
1032 : VecScale(blub, posAx, gradIntegrandVal);
1033 : number w = qPoints[j].second;
1034 :
1035 : defectOut += integrandVal * w * (end-start);
1036 : VecScaleAdd(gradientOut, 1.0, gradientOut, w * (end-start), blub);
1037 :
1038 : UG_LOGN("posNeurite: " << posAx);
1039 : UG_LOGN("Neurite part information... radius: " << radius << ", "
1040 : << "x: " << x << ", vel: " << vel << ", posAx: " << posAx
1041 : << ",radInv: "<< radInv << ", velNorm: " << velNorm
1042 : << ",posNormSq: " << posNormSq << ", w:" << w << ", j: "
1043 : << j << ", start: " << start << ", end: " << end
1044 : << ",defectOut neurite: " << integrandVal * w * (end-start)
1045 : << ",integrand neurite:" << integrandVal
1046 : << ",gradintegrand neurite: " << gradIntegrandVal
1047 : << ",blub: " << blub);
1048 : }
1049 : UG_LOGN("gradientOut (neurite): " << gradientOut);
1050 :
1051 : /// SOMA surface PM part or ER surface of inner sphere part
1052 : vector3 posAxSoma;
1053 : UG_LOGN("posSoma (initial): " << posSoma);
1054 : /// VecSubtract(posAxSoma, posSoma, x);
1055 : /// TODO: Verify: Soma position means: in direction from soma center towards
1056 : /// current point position (x) scaled with t (axial) equals position!
1057 : vector3 dir;
1058 : VecSubtract(dir, x, posSoma);
1059 : // FIXME verify this: soma outer sphere axial = 0 inner ER sphere -1:
1060 : /// set position of soma between end and start region of integration
1061 : const number t = std::fabs(helper.sr->t) + helper.sr->t+(start-end)/2;
1062 : VecScale(dir, dir, 1+t/helper.sr->radius);
1063 : VecAdd(posAxSoma, posSoma, dir);
1064 : VecSubtract(posAxSoma, posAxSoma, x);
1065 : vector3 velSoma = posAxSoma;
1066 : number velSomaNorm = sqrt(VecNormSquared(velSoma));
1067 : number posNormSomaSq = VecNormSquared(posAxSoma);
1068 : UG_LOGN("t (Soma): " << t);
1069 : UG_LOGN("somaPos: " << posAxSoma);
1070 : number radInvSoma = 1.0 / ((1+t)*somaRadius);
1071 : number gradIntegrandValSoma = radInvSoma * exp(-posNormSomaSq*radInvSoma*radInvSoma) * velSomaNorm;
1072 : vector3 blubSoma;
1073 : VecScale(blubSoma, posAxSoma, gradIntegrandValSoma);
1074 : VecScaleAdd(gradientOut, 1.0, gradientOut, 1.0 * (end-start), blubSoma);
1075 :
1076 : defectOut += gradIntegrandValSoma * 1.0 * (end-start);
1077 : UG_LOGN("defectOut (soma):" << gradIntegrandValSoma * 1.0 * (end-start));
1078 :
1079 : UG_LOGN("Soma part information... radius: " << (1+t)*somaRadius << ", "
1080 : << "x: " << x << ", vel: " << velSoma << ", posAx: " << posAxSoma
1081 : << ",radInv: "<< radInvSoma << ", velNorm: " << velSomaNorm
1082 : << ",posNormSq: " << posNormSomaSq << ", w:" << 1.0 << ", j: "
1083 : << 1 << ", start: " << start << ", end: " << end
1084 : << ",defectOut neurite: " << gradIntegrandValSoma * 1.0 * (end-start)
1085 : << ",gradIntegrand neurite:" << gradIntegrandValSoma
1086 : << ",blub: " << blubSoma);
1087 :
1088 : /// TODO: Is the defect calculation correct? It seems we do not account for the added soma point to the integrand?
1089 : defectOut -= sqrt(PI)*exp(-1.0);
1090 : UG_LOGN("defectOut (total): " << defectOut);
1091 : UG_LOGN("gradientOut (total): " << gradientOut);
1092 :
1093 : // TODO: project gradient to surface of constant angle: Is this still correct to use after adding the soma point to the integrand?
1094 : VecScaleAppend(gradientOut, -VecDot(constAngleSurfNormal, gradientOut), constAngleSurfNormal);
1095 : }
1096 :
1097 :
1098 : static void newton_for_soma_bp_projection
1099 : (
1100 : vector3& posOut,
1101 : const std::vector<NeuriteProjector::BPProjectionHelper>& vProjHelp,
1102 : const NeuriteProjector* np,
1103 : const vector3& constAngleSurfNormal
1104 : )
1105 : {
1106 : vector3 posStart = posOut;
1107 : UG_LOGN("posOut (iter):" << posOut);
1108 : // calculate start defect and gradient
1109 : number def;
1110 : number def_init;
1111 : vector3 grad;
1112 : bp_defect_and_gradient_soma(def, grad, vProjHelp, posOut, np, constAngleSurfNormal);
1113 : def_init = fabs(def);
1114 : // if initially gradient is zero, better starting position has to be chosen
1115 : UG_LOGN("def initial: " << def);
1116 : UG_LOGN("gradient initial: " << grad);
1117 :
1118 : // perform some kind of Newton search for the correct position
1119 : size_t maxIt = 100;
1120 : number minDef = 1e-8*sqrt(PI)*exp(-1.0);
1121 : size_t iter = 0;
1122 : number corr = 1.0;
1123 :
1124 : while (fabs(def) > minDef && fabs(def) > 1e-8 * def_init && ++iter <= maxIt)
1125 : {
1126 : UG_LOGN("First while loop")
1127 : vector3 posTest;
1128 : vector3 gradTest;
1129 : number defTest;
1130 : // start with reduced step size to prevent overshooting in case the start position is too bad
1131 : corr /= 2;
1132 : UG_LOGN("-def: " << -def);
1133 : UG_LOGN("VecNormSquared(grad): " << VecNormSquared(grad));
1134 : UG_LOGN("grad: " << grad);
1135 : VecScaleAdd(posTest, 1.0, posOut, -(1.0-corr)*def / VecNormSquared(grad), grad);
1136 : bp_defect_and_gradient_soma(defTest, gradTest, vProjHelp, posTest, np, constAngleSurfNormal);
1137 : UG_COND_THROW(std::fabs(VecNorm(gradTest)) < SMALL, "Quasi-zero gradient: Iteration will always fail.");
1138 :
1139 : // line search
1140 : number lambda = 0.5;
1141 : number bestDef = defTest;
1142 : vector3 bestGrad = gradTest;
1143 : vector3 bestPos = posTest;
1144 : while (fabs(defTest) >= fabs(def) && lambda > 0.001)
1145 : {
1146 : UG_LOGN("Second while loop")
1147 : VecScaleAdd(posTest, 1.0, posOut, -lambda*def / VecNormSquared(grad), grad);
1148 : bp_defect_and_gradient_soma(defTest, gradTest, vProjHelp, posTest, np, constAngleSurfNormal);
1149 : if (fabs(defTest) < fabs(bestDef))
1150 : {
1151 : bestDef = defTest;
1152 : bestGrad = gradTest;
1153 : bestPos = posTest;
1154 : }
1155 : lambda *= 0.5;
1156 : }
1157 : def = bestDef;
1158 : grad = gradTest;
1159 : posOut = bestPos;
1160 : }
1161 : UG_COND_THROW(def != def,
1162 : "Newton iteration did not converge for soma branching point projection at " << posStart << ". Defect is NaN.")
1163 : UG_COND_THROW(fabs(def) > minDef && fabs(def) > 1e-8 * def_init,
1164 : "Newton iteration did not converge for soma branching point projection at " << posStart << ".")
1165 : }
1166 : #endif
1167 :
1168 :
1169 : #if 0
1170 : static void pos_on_neurite_spine
1171 : (
1172 : vector3& posOut,
1173 : const NeuriteProjector::Neurite& neurite,
1174 : size_t neuriteID,
1175 : float& t
1176 : )
1177 : {
1178 : const std::vector<NeuriteProjector::Section>& vSections = neurite.vSec;
1179 :
1180 : // find correct section
1181 : NeuriteProjector::Section cmpSec(t);
1182 : std::vector<NeuriteProjector::Section>::const_iterator it =
1183 : std::lower_bound(vSections.begin(), vSections.end(), cmpSec, NeuriteProjector::CompareSections());
1184 :
1185 : UG_COND_THROW(it == vSections.end(),
1186 : "Could not find section for parameter t = " << t << " in neurite " << neuriteID << ".");
1187 :
1188 : // calculate correct axial position
1189 : // and velocity dir
1190 : vector3 vel_dummy;
1191 : number radius_dummy;
1192 : compute_position_and_velocity_in_section(posOut, vel_dummy, radius_dummy, it, t);
1193 : }
1194 : #endif
1195 :
1196 :
1197 0 : static void newton_for_bp_projection
1198 : (
1199 : vector3& posOut,
1200 : const std::vector<NeuriteProjector::BPProjectionHelper>& vProjHelp,
1201 : float rad,
1202 : const vector3& constAngleSurfNormal,
1203 : const NeuriteProjector* np
1204 : )
1205 : {
1206 : vector3 posStart = posOut;
1207 :
1208 : // calculate start defect and gradient
1209 : number def;
1210 : number def_init;
1211 : vector3 grad;
1212 :
1213 0 : bp_defect_and_gradient(def, grad, vProjHelp, posOut, rad, constAngleSurfNormal, np);
1214 0 : def_init = fabs(def);
1215 :
1216 : // perform some kind of Newton search for the correct position
1217 : size_t maxIt = 10;
1218 : number minDef = 1e-8*sqrt(PI)*exp(-1.0);
1219 : size_t iter = 0;
1220 : number corr = 1.0;
1221 :
1222 : //UG_LOGN(iter << " " << def << " " << grad << " " << posOut);
1223 0 : while (fabs(def) > minDef && fabs(def) > 1e-8 * def_init && ++iter <= maxIt)
1224 : {
1225 0 : corr /= 2; // start with reduced step size to prevent overshooting in case the start position is too bad
1226 : vector3 posTest;
1227 : vector3 gradTest;
1228 : number defTest;
1229 0 : VecScaleAdd(posTest, 1.0, posOut, -(1.0-corr)*def / VecNormSquared(grad), grad);
1230 0 : bp_defect_and_gradient(defTest, gradTest, vProjHelp, posTest, rad, constAngleSurfNormal, np);
1231 :
1232 : // line search
1233 : number lambda = 0.5;
1234 0 : number bestDef = defTest;
1235 : vector3 bestGrad = gradTest;
1236 : vector3 bestPos = posTest;
1237 0 : while (fabs(defTest) >= fabs(def) && lambda > 0.001)
1238 : {
1239 : //UG_LOGN(" line search with lambda = " << lambda);
1240 0 : VecScaleAdd(posTest, 1.0, posOut, -lambda*def / VecNormSquared(grad), grad);
1241 0 : bp_defect_and_gradient(defTest, gradTest, vProjHelp, posTest, rad, constAngleSurfNormal, np);
1242 0 : if (fabs(defTest) < fabs(bestDef))
1243 : {
1244 : bestDef = defTest;
1245 : bestGrad = gradTest;
1246 : bestPos = posTest;
1247 : }
1248 0 : lambda *= 0.5;
1249 : }
1250 0 : def = bestDef;
1251 : grad = gradTest;
1252 : posOut = bestPos;
1253 :
1254 : //UG_LOGN(iter << " " << def << " " << grad << " " << posOut);
1255 : }
1256 0 : UG_COND_THROW(def != def,
1257 : "Newton iteration did not converge for branching point projection at " << posStart << ". Defect is NaN.")
1258 0 : UG_COND_THROW(fabs(def) > minDef && fabs(def) > 1e-8 * def_init,
1259 : "Newton iteration did not converge for branching point projection at " << posStart << ".")
1260 0 : }
1261 :
1262 :
1263 0 : static void pos_in_neurite
1264 : (
1265 : vector3& posOut,
1266 : const NeuriteProjector::Neurite& neurite,
1267 : size_t neuriteID,
1268 : float t,
1269 : float angle,
1270 : float rad
1271 : )
1272 : {
1273 : /// TODO: A similiar procedure as below has to be used for the soma/neurite BP iteration
1274 : const std::vector<NeuriteProjector::Section>& vSections = neurite.vSec;
1275 :
1276 : // find correct section
1277 0 : NeuriteProjector::Section cmpSec(t);
1278 : std::vector<NeuriteProjector::Section>::const_iterator it =
1279 0 : std::lower_bound(vSections.begin(), vSections.end(), cmpSec, NeuriteProjector::CompareSections());
1280 :
1281 0 : UG_COND_THROW(it == vSections.end(),
1282 : "Could not find section for parameter t = " << t << " in neurite " << neuriteID << ".");
1283 :
1284 : // calculate correct axial position
1285 : // and velocity dir
1286 : vector3 posAx, vel;
1287 : number radius;
1288 0 : compute_position_and_velocity_in_section(posAx, vel, radius, it, t);
1289 :
1290 : // calculate orthonormal basis vectors
1291 : vector3 projRefDir, thirdDir;
1292 0 : compute_ONB(vel, projRefDir, thirdDir, vel, neurite.refDir);
1293 :
1294 : // calculate new position
1295 0 : VecScaleAdd(posOut, 1.0, posAx, rad*radius*cos(angle), projRefDir, rad*radius*sin(angle), thirdDir);
1296 0 : }
1297 :
1298 :
1299 0 : static void bp_newton_start_pos
1300 : (
1301 : vector3& initOut,
1302 : uint32_t neuriteID,
1303 : float t,
1304 : float angle,
1305 : float rad,
1306 : const vector3& constAngleSurfaceNormal,
1307 : const IVertexGroup* parent,
1308 : const NeuriteProjector* np
1309 : )
1310 : {
1311 : // Usually, we will want to just use the average position v_avg of all
1312 : // parent element vertices as start position for the Newton iteration.
1313 : // However, if a dendrite is strongly bent in the vicinity of a BP,
1314 : // this position can be quite far away from the final position.
1315 : // In these cases, we rather want to start in the position v_spl defined
1316 : // in the dendrite by the axial, radial and angular coordinates.
1317 :
1318 : // In order to detect these cases, we calculate
1319 : // q = ||v_avg - v_spl|| / r(t),
1320 : // where r(t) is the dendrite radius at the axial position of interest,
1321 : // as well as
1322 : // p = diam(parent) / r(t)
1323 : // where diam(parent) is the diameter of the parent element.
1324 : // This is used to distinguish a situation as described above from one
1325 : // where we are in the middle of a BP and q is large, but we still want
1326 : // to use simple averaging.
1327 :
1328 : // To avoid hard switching between cases, we calculate a continuous weighting
1329 : // function w and use
1330 : // w*v_spl + (1-w)*v_avg
1331 : // as the starting position.
1332 :
1333 : vector3 pos_avg_onNeurite;
1334 0 : pos_in_neurite(pos_avg_onNeurite, np->neurite(neuriteID), neuriteID, t, angle, rad);
1335 :
1336 : // special case: no parents (happens when a vertex is being projected without refinement)
1337 0 : if (!parent)
1338 : {
1339 : // start at position given by the parameterization
1340 : initOut = pos_avg_onNeurite;
1341 0 : return;
1342 : }
1343 :
1344 : // find correct section
1345 0 : const std::vector<NeuriteProjector::Section>& vSections = np->neurite(neuriteID).vSec;
1346 0 : NeuriteProjector::Section cmpSec(t);
1347 : std::vector<NeuriteProjector::Section>::const_iterator it =
1348 0 : std::lower_bound(vSections.begin(), vSections.end(), cmpSec, NeuriteProjector::CompareSections());
1349 0 : UG_COND_THROW(it == vSections.end(),
1350 : "Could not find section for parameter t = " << t << " in neurite " << neuriteID << ".");
1351 :
1352 : // calculate radius
1353 : number radius = 0.0;
1354 : compute_radius_in_section(radius, it, t);
1355 :
1356 : // calculate element diameter
1357 0 : number diamSq = 0.0;
1358 : const size_t nVrt = parent->size();
1359 0 : for (size_t i = 0; i < nVrt; ++i)
1360 : {
1361 0 : const vector3& posi = np->geometry()->pos(parent->vertex(i));
1362 0 : for (size_t j = i+1; j < nVrt; ++j)
1363 0 : diamSq = std::max(diamSq, VecDistanceSq(posi, np->geometry()->pos(parent->vertex(j))));
1364 : }
1365 :
1366 : vector3 pos_avg;
1367 0 : np->average_pos_from_parent(pos_avg, parent);
1368 :
1369 0 : number distSq = VecDistanceSq(pos_avg_onNeurite, pos_avg);
1370 :
1371 : number qSq = 0.0;
1372 0 : if (radius*radius <= 1e-8*distSq)
1373 : qSq = 1e8;
1374 : else
1375 0 : qSq = distSq / (radius*radius);
1376 0 : const number w1 = qSq*qSq / (qSq*qSq + 1e-4); // in [0,1] with w1 = 0.5 for q = 0.1;
1377 :
1378 : number pSq = 0.0;
1379 0 : if (radius*radius <= 1e-8*diamSq)
1380 : pSq = 1e8;
1381 : else
1382 0 : pSq = diamSq / (radius*radius);
1383 0 : const number w2 = pSq*pSq / (pSq*pSq + 256); // in [0,1] with w2 = 0.5 for p = 4
1384 :
1385 0 : VecScaleAdd(initOut, w1*w2, pos_avg_onNeurite, 1.0-w1*w2, pos_avg);
1386 :
1387 : // project onto plane of constant angle
1388 : vector3 ptToSurf;
1389 : VecSubtract(ptToSurf, pos_avg_onNeurite, initOut);
1390 : VecScaleAdd(initOut, 1.0, initOut, VecDot(ptToSurf, constAngleSurfaceNormal), constAngleSurfaceNormal);
1391 : }
1392 :
1393 :
1394 : #if 0
1395 : static void pos_on_surface_soma_bp
1396 : (
1397 : vector3& posOut,
1398 : const NeuriteProjector::Neurite& neurite,
1399 : size_t neuriteID,
1400 : float& t,
1401 : float& angle,
1402 : const IVertexGroup* parent,
1403 : const NeuriteProjector* np,
1404 : float rad,
1405 : Vertex* vrt,
1406 : const NeuriteProjector::SomaBranchingRegion& sr,
1407 : const Grid::VertexAttachmentAccessor<Attachment<NeuriteProjector::SurfaceParams> >& aaSurfParams
1408 : ) {
1409 : // if vertex is on the neurite backbone, simply project onto spline
1410 : // (Newton iteration will not succeed for Jacobian is zero)
1411 : if (rad < 1e-5)
1412 : {
1413 : pos_in_neurite(posOut, np->neurite(neuriteID), neuriteID, t, angle, rad);
1414 : return;
1415 : }
1416 :
1417 : // 1. preparations for Newton method:
1418 : // save integration start and end positions of soma connection
1419 : // also save a section iterator to start from
1420 : std::vector<NeuriteProjector::BPProjectionHelper> vProjHelp(1);
1421 : number& start = vProjHelp[0].start;
1422 : number& end = vProjHelp[0].end;
1423 : std::vector<NeuriteProjector::Section>::const_iterator& sec_start = vProjHelp[0].sec_start;
1424 : number range = 0;
1425 : /// FIXME: This method generates a range (which is negative) far too large, c
1426 : /// reate only range up to first sections secEnd param at most and not negative
1427 : try {range = np->axial_range_around_soma_region(sr, rad, neuriteID, vrt); }
1428 : UG_CATCH_THROW("Range around soma region could not be calculated.");
1429 : UG_LOGN("range: " << range);
1430 : /// Test: Too small range will never succeed in integration. Why?!
1431 : /*
1432 : start = aaSurfParams[vrt].axial;
1433 : end = aaSurfParams[vrt].axial+0.00001;
1434 : */
1435 : /// FIXME: Integrates only outside of the soma for now now, but need to
1436 : // integrate around soma or ER start, e.g. [start-range, start+range]
1437 : /// soma starts at 0 ER starts at -1
1438 : start = 0;
1439 : end = 0.001;
1440 : UG_LOGN("Range: (start: " << start << ", end: " << end << ")");
1441 : vProjHelp[0].sr = &sr;
1442 : sec_start = (&np->neurite(neuriteID).vSec)->begin();
1443 :
1444 : // 2. determine suitable start position for Newton iteration
1445 : pos_on_surface_soma(posOut, np->neurite(neuriteID), np, parent);
1446 : /*if (parent) {
1447 : bp_newton_start_pos(posOut, neuriteID, t, angle, rad, parent, np);
1448 : }
1449 : */
1450 :
1451 : // also calculate normal of constant-angle surface
1452 : vector3 s, v;
1453 : number radius;
1454 : compute_position_and_velocity_in_section(s, v, radius, sec_start, t);
1455 : vector3 projRefDir;
1456 : vector3 thirdDir;
1457 : vector3 constAngleSurfNormal;
1458 : compute_ONB(constAngleSurfNormal, projRefDir, thirdDir, v, neurite.refDir);
1459 : VecScaleAdd(constAngleSurfNormal, sin(angle), projRefDir, -cos(angle), thirdDir);
1460 : UG_LOGN("posOut (initial): " << posOut);
1461 :
1462 : // 3. perform Newton iteration
1463 : try {newton_for_soma_bp_projection(posOut, vProjHelp, np, constAngleSurfNormal);}
1464 : UG_CATCH_THROW("Newton iteration for neurite projection at branching point failed.");
1465 : }
1466 : #endif
1467 :
1468 :
1469 0 : static void pos_in_bp
1470 : (
1471 : vector3& posOut,
1472 : const NeuriteProjector::Neurite& neurite,
1473 : uint32_t neuriteID,
1474 : float& t,
1475 : float& angle,
1476 : float rad,
1477 : std::vector<NeuriteProjector::BranchingRegion>::const_iterator& it,
1478 : const IVertexGroup* parent,
1479 : const NeuriteProjector* np
1480 : )
1481 : {
1482 : // if vertex is on the neurite backbone, simply project onto spline
1483 : // (Newton iteration will not succeed for Jacobian is zero)
1484 0 : if (rad < 1e-5)
1485 : {
1486 0 : pos_in_neurite(posOut, np->neurite(neuriteID), neuriteID, t, angle, rad);
1487 0 : return;
1488 : }
1489 :
1490 : // preparations for Newton method:
1491 : // save integration start and end positions of all BRs of this BP,
1492 : // also save a section iterator to start from
1493 : SmartPtr<NeuriteProjector::BranchingPoint> bp = it->bp;
1494 : size_t nParts = bp->vRegions.size();
1495 0 : std::vector<NeuriteProjector::BPProjectionHelper> vProjHelp(nParts);
1496 : const std::vector<NeuriteProjector::Section>* secs;
1497 0 : for (size_t i = 0; i < nParts; ++i)
1498 : {
1499 : number& start = vProjHelp[i].start;
1500 : number& end = vProjHelp[i].end;
1501 : std::vector<NeuriteProjector::Section>::const_iterator& sec_start = vProjHelp[i].sec_start;
1502 :
1503 0 : uint32_t nid = bp->vNid[i];
1504 0 : const NeuriteProjector::BranchingRegion* br = bp->vRegions[i];
1505 :
1506 0 : UG_COND_THROW(!br, "Requested branching region not accessible.");
1507 :
1508 : // integrate over twice the size of the original BR (i.e., 10 times rad*radius in each direction)
1509 : number range = 0.0;
1510 0 : try {range = np->axial_range_around_branching_region(nid,
1511 0 : (size_t) (br - &np->neurite(nid).vBR[0]), 10.0*rad);}
1512 0 : UG_CATCH_THROW("Range around branching region could not be calculated.");
1513 0 : start = std::max(br->t - range, 0.0);
1514 0 : end = std::min(br->t + range, 1.0);
1515 :
1516 0 : secs = &np->neurite(nid).vSec;
1517 0 : if (start == 0.0)
1518 0 : sec_start = secs->begin();
1519 : else
1520 : {
1521 : NeuriteProjector::Section cmpSec(start);
1522 0 : sec_start = std::lower_bound(secs->begin(), secs->end(), cmpSec, NeuriteProjector::CompareSections());
1523 0 : UG_COND_THROW(sec_start == secs->end(),
1524 : "Could not find section for start parameter t = " << start << ".");
1525 : }
1526 : }
1527 :
1528 : // also calculate normal of constant-angle surface
1529 : vector3 constAngleSurfNormal;
1530 0 : NeuriteProjector::Section cmpSec(t);
1531 : const std::vector<NeuriteProjector::Section>& vSections = neurite.vSec;
1532 : std::vector<NeuriteProjector::Section>::const_iterator secIt =
1533 0 : std::lower_bound(vSections.begin(), vSections.end(), cmpSec, NeuriteProjector::CompareSections());
1534 :
1535 : vector3 s, v;
1536 : number radius;
1537 0 : compute_position_and_velocity_in_section(s, v, radius, secIt, t);
1538 : vector3 projRefDir;
1539 : vector3 thirdDir;
1540 0 : compute_ONB(constAngleSurfNormal, projRefDir, thirdDir, v, neurite.refDir);
1541 0 : VecScaleAdd(constAngleSurfNormal, sin(angle), projRefDir, -cos(angle), thirdDir);
1542 :
1543 : // determine suitable start position for Newton iteration
1544 0 : bp_newton_start_pos(posOut, neuriteID, t, angle, rad, constAngleSurfNormal, parent, np);
1545 : //if (parent)
1546 : // np->average_pos_from_parent(posOut, parent);
1547 : //pos_in_neurite(posOut, np->neurite(neuriteID), neuriteID, t, angle, rad);
1548 :
1549 : // perform Newton iteration
1550 0 : try {newton_for_bp_projection(posOut, vProjHelp, rad, constAngleSurfNormal, np);}
1551 0 : UG_CATCH_THROW("Newton iteration for neurite projection at branching point failed"
1552 : " in neurite " << neuriteID << ".");
1553 :
1554 :
1555 : // update the surface param info to the new position
1556 : // In order not to have to compute zeros of a 5th order polynomial,
1557 : // we make a linearity approximation here:
1558 : // s(t) = s(t0) + v(t0)*(t-t0);
1559 : // v(t) = v(t0)
1560 : // Then from v(t) * (s(t)-pos) = 0, we get:
1561 : // t = t0 + v(t0) * (pos - s(t0)) / (v(t0)*v(t0))
1562 : VecSubtract(s, posOut, s);
1563 0 : t += VecProd(v,s) / VecProd(v,v);
1564 :
1565 : // and now angle
1566 0 : if (rad > 1e-5)
1567 0 : angle = compute_angle(v, projRefDir, thirdDir, s);
1568 :
1569 : // and radius
1570 : //rad = VecNorm(s) / radius;
1571 :
1572 : /*
1573 : // if we have ended up outside of the BP, then use the usual positioning
1574 : if (t > it->tend)
1575 : {
1576 : vector3 posAx;
1577 : number radius;
1578 : compute_position_and_velocity_in_section(posAx, v, radius, secIt, t);
1579 : VecScaleAdd(posOut, 1.0, posAx, rad*radius*cos(angle), projRefDir, rad*radius*sin(angle), thirdDir);
1580 : }
1581 : */
1582 0 : }
1583 :
1584 :
1585 0 : static void pos_on_surface_tip
1586 : (
1587 : vector3& posOut,
1588 : const NeuriteProjector::Neurite& neurite,
1589 : const IVertexGroup* parent,
1590 : const NeuriteProjector* np,
1591 : number rad
1592 : )
1593 : {
1594 : const std::vector<NeuriteProjector::Section>& vSections = neurite.vSec;
1595 :
1596 : // initial position: regular refinement
1597 : // in case there is no parent (projection of an existing vertex), keep the position
1598 0 : if (parent)
1599 0 : np->average_pos_from_parent(posOut, parent);
1600 :
1601 : // project to half-sphere with given radius over tip center
1602 : // TODO: One might think about something more sophisticated,
1603 : // as the current approach ensures only continuity of the radius
1604 : // at tips, but not differentiability.
1605 : const NeuriteProjector::Section& sec = vSections[vSections.size()-1];
1606 0 : vector3 tipCenter(sec.splineParamsX[3], sec.splineParamsY[3], sec.splineParamsZ[3]);
1607 0 : number radius = sec.splineParamsR[3];
1608 :
1609 : vector3 radialVec;
1610 : VecSubtract(radialVec, posOut, tipCenter);
1611 0 : number radNorm = sqrt(VecProd(radialVec, radialVec));
1612 :
1613 : // only project if we are not in the center of the half sphere
1614 0 : if (radNorm >= 1e-5*rad*radius)
1615 : {
1616 0 : VecScale(radialVec, radialVec, rad*radius/radNorm);
1617 : VecAdd(posOut, tipCenter, radialVec);
1618 : }
1619 0 : }
1620 :
1621 0 : number NeuriteProjector::axial_range_around_soma_region
1622 : (
1623 : const SomaBranchingRegion& sr,
1624 : const size_t numRadii,
1625 : const size_t nid,
1626 : Vertex* vrt
1627 : ) const
1628 : {
1629 : // soma parameters
1630 : ///const float somaRad = sr.radius;
1631 : vector3 bp = sr.bp;
1632 : std::vector<NeuriteProjector::Section>::const_iterator secIt;
1633 0 : try {secIt = get_soma_section_iterator(nid, sr.t);}
1634 0 : UG_CATCH_THROW("Could not get section iterator to soma region: " << sr.t);
1635 : vector3 secStartPos;
1636 : vector3 secEndPos;
1637 0 : const number te = secIt->endParam;
1638 : number ts = 0.0;
1639 0 : compute_first_section_end_points(secStartPos, secEndPos, secIt);
1640 0 : const number neuriteLengthApprox = VecDistance(secEndPos, secStartPos) / (te - ts);
1641 : const number* s = &secIt->splineParamsR[0];
1642 0 : number radius = s[0]*(te-sr.t) + s[1];
1643 0 : radius = radius*(te-sr.t) + s[2];
1644 0 : radius = radius*(te-sr.t) + s[3];
1645 0 : const number bpToVrt = VecDistance(bp, position(vrt));
1646 0 : const number neurite_length_to_soma_bp= neuriteLengthApprox + bpToVrt;
1647 0 : const number rad = (numRadii * radius)/neurite_length_to_soma_bp;
1648 0 : return rad;
1649 : }
1650 :
1651 0 : bool NeuriteProjector::is_in_axial_range_around_soma_region
1652 : (
1653 : const SomaBranchingRegion& sr,
1654 : size_t nid,
1655 : Vertex* vrt
1656 : ) const {
1657 : ///const number radius = axial_range_around_soma_region(sr, 0.1, nid, vrt);
1658 : const number radius = 0.1;
1659 0 : UG_LOGN("Soma Branching Region center: " << sr.center);
1660 : UG_LOGN("Soma Branching Region raidus: " << radius);
1661 0 : return IsElementInsideSphere<Vertex>(vrt, sr.center, radius);
1662 : }
1663 :
1664 0 : number NeuriteProjector::push_into_place(Vertex* vrt, const IVertexGroup* parent)
1665 : {
1666 : // average axial and angular params from parent;
1667 : // also decide on neuriteID for aaSurfParams.
1668 : // the mapping parameters 1d<->2d are not needed on higher levels of refinment
1669 : uint32_t neuriteID;
1670 : float t;
1671 : float angle;
1672 : float rad;
1673 0 : if (parent) {
1674 0 : average_params(neuriteID, t, angle, rad, parent);
1675 : } else {
1676 0 : neuriteID = m_aaSurfParams[vrt].neuriteID;
1677 0 : t = m_aaSurfParams[vrt].axial;
1678 0 : angle = m_aaSurfParams[vrt].angular;
1679 0 : rad = m_aaSurfParams[vrt].radial;
1680 : }
1681 :
1682 : //UG_LOGN("averaged params: nid " << neuriteID << " t " << t << " angle "
1683 : // << angle << " rad " << rad);
1684 :
1685 0 : const uint32_t plainNID = (neuriteID << 12) >> 12;
1686 :
1687 0 : UG_COND_THROW(plainNID >= m_vNeurites.size(), "Requested neurite ID which is not present.");
1688 :
1689 : const Neurite& neurite = m_vNeurites[plainNID];
1690 : const std::vector<BranchingRegion>& vBR = neurite.vBR;
1691 : const std::vector<SomaBranchingRegion>& vSBR = neurite.vSBR;
1692 :
1693 : ///UG_LOGN("Number of soma branching regions for neurite with id: " << plainNID << ": " << vSBR.size());
1694 :
1695 :
1696 :
1697 : // vector for new position
1698 0 : vector3 pos(position(vrt));
1699 :
1700 : // FIVE CASES can occur:
1701 : // 1. We are at a branching point.
1702 : // 2. We are at a soma branching point
1703 : // 3. We are well inside a regular piece of neurite.
1704 : // 4. We are well inside a regular piece of soma.
1705 : // 5. We are at the tip of a neurite.
1706 :
1707 : /// This works due to the fact that vBR is sorted ascending
1708 : bool isBP = false;
1709 0 : BranchingRegion cmpBR(t);
1710 : std::vector<BranchingRegion>::const_iterator it =
1711 0 : std::upper_bound(vBR.begin(), vBR.end(), cmpBR, CompareBranchingRegionEnds());
1712 :
1713 : /// This works due to the fact that vSBR is sorted ascending
1714 : bool isSP = false;
1715 0 : SomaBranchingRegion cmpSBR(t);
1716 : std::vector<SomaBranchingRegion>::const_iterator it2 =
1717 0 : std::lower_bound(vSBR.begin(), vSBR.end(), cmpSBR, CompareSomaBranchingRegionsEnd());
1718 :
1719 0 : if (it2 != vSBR.end()) {
1720 0 : isSP = is_in_axial_range_around_soma_region(*it2, plainNID, vrt);
1721 0 : if (isSP) {
1722 0 : UG_LOGN("Soma branching point: " << this->pos(vrt));
1723 : }
1724 : }
1725 :
1726 0 : if (it != vBR.end())
1727 : {
1728 : number range = 0.0;
1729 0 : try {range = axial_range_around_branching_region(plainNID, std::distance(vBR.begin(), it), 5.0*rad);}
1730 0 : UG_CATCH_THROW("Range around branching region could not be determined.");
1731 0 : if (it->t - t < range)
1732 : {
1733 : isBP = true;
1734 : }
1735 : }
1736 :
1737 0 : if (!isBP && it != vBR.begin())
1738 : {
1739 : --it;
1740 : number range = 0.0;
1741 0 : try {range = axial_range_around_branching_region(plainNID, std::distance(vBR.begin(), it), 5.0*rad);}
1742 0 : UG_CATCH_THROW("Range around branching region could not be determined.");
1743 0 : if (t - it->t < range)
1744 : {
1745 : isBP = true;
1746 : }
1747 : }
1748 :
1749 : // case 1: branching point
1750 0 : if (isBP)
1751 : {
1752 0 : pos_in_bp(pos, neurite, plainNID, t, angle, rad, it, parent, this);
1753 : }
1754 :
1755 : // case 2: soma/neurite or er/neurite branching point
1756 0 : else if (isSP)
1757 : {
1758 : /// This should integrate along the neurite [0, END_SOMA_BRANCHING_REGION]
1759 : /// but not [-range, range] to avoid dints in the soma interior close to
1760 : /// the surface. A new position for the vertex based on the average of the
1761 : /// (first) neurite section and the soma sphere position and radius will
1762 : /// be calculated based on the SomaRegion information stored in a struct.
1763 : /// FIXME: Optimize iteration: Find better starting position for iteration start
1764 : ///pos_on_surface_soma_bp(pos, neurite, neuriteID, t, angle, parent, this, rad, vrt, *it2, m_aaSurfParams);
1765 : }
1766 :
1767 : // case 3: normal neurite position
1768 0 : else if (t >= 0.0 && t <= 1.0)
1769 : {
1770 0 : pos_in_neurite(pos, neurite, plainNID, t, angle, rad);
1771 : }
1772 :
1773 : // case 4: normal soma position
1774 0 : else if (t < 0 && t >= -1.0)
1775 : {
1776 : ///pos_on_surface_soma(pos, neurite, this, parent);
1777 : }
1778 :
1779 : // case 5: tip of neurite
1780 : else
1781 : {
1782 0 : pos_on_surface_tip(pos, neurite, parent, this, rad);
1783 : }
1784 :
1785 : // save new surface params for new vertex
1786 0 : m_aaSurfParams[vrt].neuriteID = neuriteID;
1787 0 : m_aaSurfParams[vrt].axial = t;
1788 0 : m_aaSurfParams[vrt].angular = angle;
1789 0 : m_aaSurfParams[vrt].radial = rad;
1790 :
1791 : // set position
1792 : set_pos(vrt, pos);
1793 :
1794 0 : return 1.0;
1795 : }
1796 :
1797 :
1798 :
1799 : // -------------------------------------------------------- //
1800 : // DO NOT CHANGE below this line! Needed for serialization. //
1801 :
1802 0 : std::ostream& operator<<(std::ostream& os, const NeuriteProjector::SurfaceParams& surfParams)
1803 : {
1804 : using std::ostringstream;
1805 0 : ostringstream strs;
1806 0 : strs << surfParams.neuriteID << " ";
1807 0 : strs << surfParams.axial << " ";
1808 0 : strs << surfParams.angular << " ";
1809 0 : strs << surfParams.radial;
1810 0 : os << strs.str();
1811 0 : return os;
1812 0 : }
1813 :
1814 0 : std::istream& operator>>(std::istream& in, NeuriteProjector::SurfaceParams& surfParams)
1815 : {
1816 : std::string temp;
1817 : using boost::lexical_cast;
1818 :
1819 0 : in >> temp;
1820 0 : surfParams.neuriteID = lexical_cast<uint32>(temp);
1821 : temp.clear();
1822 :
1823 0 : in >> temp;
1824 0 : surfParams.axial = (lexical_cast<float>(temp));
1825 : temp.clear();
1826 :
1827 0 : in >> temp;
1828 0 : surfParams.angular = lexical_cast<float>(temp);
1829 : temp.clear();
1830 :
1831 0 : in >> temp;
1832 0 : surfParams.radial = lexical_cast<float>(temp);
1833 : temp.clear();
1834 :
1835 0 : return in;
1836 : }
1837 :
1838 0 : std::ostream& operator<<(std::ostream& os, const NeuriteProjector::Mapping& mapping)
1839 : {
1840 : /// Standard precision for UGX coord export is 18. Is this even correct to use?
1841 : /// Shouldn't one use std::numeric_limits<number>::digits10+1?
1842 : using std::ostringstream;
1843 0 : ostringstream strs;
1844 0 : for (size_t i = 0; i < mapping.v1.size(); i++) {
1845 0 : strs << std::setprecision(18) << mapping.v1[i] << " ";
1846 : }
1847 0 : for (size_t i = 0; i < mapping.v1.size(); i++) {
1848 0 : strs << std::setprecision(18) << mapping.v2[i] << " ";
1849 : }
1850 0 : strs << mapping.lambda;
1851 0 : os << strs.str();
1852 0 : return os;
1853 0 : }
1854 :
1855 0 : std::istream& operator>>(std::istream& in, NeuriteProjector::Mapping& mapping)
1856 : {
1857 : std::string temp;
1858 : using boost::lexical_cast;
1859 :
1860 0 : for (size_t i = 0; i < mapping.v1.size(); i++) {
1861 0 : in >> temp;
1862 0 : mapping.v1[i] = (lexical_cast<number>(temp));
1863 : }
1864 :
1865 0 : for (size_t i = 0; i < mapping.v2.size(); i++) {
1866 0 : in >> temp;
1867 0 : mapping.v2[i] = (lexical_cast<number>(temp));
1868 : }
1869 :
1870 0 : in >> temp;
1871 0 : mapping.lambda = (lexical_cast<number>(temp));
1872 : temp.clear();
1873 :
1874 0 : return in;
1875 : }
1876 :
1877 :
1878 : } // namespace ug
|
