Line data Source code
1 : /*
2 : * Copyright (c) 2018: G-CSC, Goethe University Frankfurt
3 : * Author: Stephan Grein
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):
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 :
33 : #ifndef __H__UG_cylinder_soma_projector
34 : #define __H__UG_cylinder_soma_projector
35 :
36 : #include "common/math/misc/math_util.h"
37 : #include "refinement_projector.h"
38 : #include "lib_grid/tools/copy_attachment_handler.h"
39 : #include "lib_grid/global_attachments.h"
40 :
41 : namespace ug {
42 : /// Projects new vertices onto a sphere during refinement.
43 : /** For projection during refinement the radius property is ignored. Instead
44 : * the distance to the center of a newly inserted vertex is calculated
45 : * as the average distance of the vertices of the parent element to the center.
46 : * The radius property thus defaults to -1.
47 : *
48 : * You may still specify a radius. This radius can be used for auto-fitting of
49 : * the center and for reprojecting a set of vertices onto the sphere.
50 : *
51 : * Only vertices which are at the soma (axial parameter: -1) are projected
52 : * on the surface. This works well for good natured connecting neurites.
53 : */
54 : class SomaProjector : public RefinementProjector {
55 : private:
56 : Attachment<NeuriteProjector::SurfaceParams> m_aSurfParams;
57 : Grid::VertexAttachmentAccessor<Attachment<NeuriteProjector::SurfaceParams> > m_aaSurfParams;
58 : CopyAttachmentHandler<Vertex, Attachment<NeuriteProjector::SurfaceParams> > m_cah;
59 : protected:
60 0 : virtual void set_geometry(SPIGeometry3d geometry)
61 : {
62 : // call base class method
63 0 : RefinementProjector::set_geometry(geometry);
64 0 : attach_surf_params();
65 0 : }
66 : public:
67 0 : SomaProjector () :
68 : m_center (0, 0, 0),
69 : m_axis (0, 0, 1),
70 : m_soma (0, 0, 0),
71 0 : m_somaRad(0)
72 : {
73 : typedef Attachment<NeuriteProjector::SurfaceParams> NPSurfParam;
74 0 : if (!GlobalAttachments::is_declared("npSurfParams"))
75 0 : GlobalAttachments::declare_attachment<NPSurfParam>("npSurfParams", true);
76 :
77 0 : }
78 :
79 : SomaProjector (const vector3& center,
80 : const vector3& axis,
81 : const vector3& soma,
82 : const number& rad) :
83 :
84 : m_center (center),
85 : m_axis (axis),
86 : m_soma (soma),
87 : m_somaRad(rad)
88 : {
89 : typedef Attachment<NeuriteProjector::SurfaceParams> NPSurfParam;
90 : if (!GlobalAttachments::is_declared("npSurfParams"))
91 : GlobalAttachments::declare_attachment<NPSurfParam>("npSurfParams", true);
92 : }
93 :
94 : public:
95 : /** \sa ug::RefinementProjector::RefinementProjector*/
96 : SomaProjector (SPIGeometry3d geometry,
97 : const vector3& center,
98 : const vector3& axis,
99 : const vector3& soma,
100 : const number& rad) :
101 : RefinementProjector (geometry),
102 : m_center (center),
103 : m_axis (axis),
104 : m_soma (soma),
105 : m_somaRad(rad)
106 : {
107 : attach_surf_params();
108 : }
109 :
110 :
111 : public:
112 0 : virtual ~SomaProjector () {}
113 :
114 : /// called when a new vertex was created from an old edge.
115 0 : virtual number new_vertex(Vertex* vrt, Edge* parent)
116 : {
117 0 : return perform_projection(vrt, parent);
118 : }
119 :
120 : /// called when a new vertex was created from an old face.
121 0 : virtual number new_vertex(Vertex* vrt, Face* parent)
122 : {
123 0 : return perform_projection(vrt, parent);
124 : }
125 :
126 : /// called when a new vertex was created from an old volume.
127 0 : virtual number new_vertex(Vertex* vrt, Volume* parent)
128 : {
129 0 : return perform_projection(vrt, parent);
130 : }
131 :
132 : private:
133 : /// check if global attachment was declared
134 0 : void check_attachment() {
135 : typedef Attachment<NeuriteProjector::SurfaceParams> NPSurfParam;
136 0 : if (!GlobalAttachments::is_declared("npSurfParams"))
137 0 : GlobalAttachments::declare_attachment<NPSurfParam>("npSurfParams", true);
138 0 : }
139 :
140 :
141 : /// helper method to attach surface parameters and do a consistency check
142 0 : void attach_surf_params() {
143 0 : check_attachment();
144 :
145 0 : Grid& grid = this->geometry()->grid();
146 :
147 : // make sure attachment was declared
148 0 : UG_COND_THROW(!GlobalAttachments::is_declared("npSurfParams"), "GlobalAttachment 'npSurfParams' not declared.");
149 0 : m_aSurfParams = GlobalAttachments::attachment<Attachment<NeuriteProjector::SurfaceParams> >("npSurfParams");
150 :
151 : // make sure surfaceParams attachment is attached
152 0 : UG_COND_THROW(!grid.has_vertex_attachment(m_aSurfParams),
153 : "No surface parameter attachment for neurite projector attached to grid.");
154 0 : m_aaSurfParams.access(grid, m_aSurfParams);
155 :
156 : // handle attachment values also on higher grid levels (if required)
157 0 : MultiGrid* mg = dynamic_cast<MultiGrid*>(&grid);
158 0 : if (mg) {
159 : SmartPtr<MultiGrid> spMG(mg);
160 :
161 : // never destroy the grid from here - we did not create it
162 0 : ++(*spMG.refcount_ptr());
163 :
164 : // attach copy attachment handler to propagate attachment to higher levels
165 : m_cah.set_attachment(m_aSurfParams);
166 0 : m_cah.set_grid(spMG);
167 : }
168 0 : }
169 :
170 :
171 : template <class TElem>
172 0 : number perform_projection(Vertex* vrt, TElem* parent)
173 : {
174 : // calculate the new position by linear interpolation and project that point
175 : // onto the cylinder.
176 0 : typename TElem::ConstVertexArray vrts = parent->vertices();
177 0 : size_t numVrts = parent->num_vertices();
178 :
179 0 : if(numVrts == 0){
180 0 : set_pos(vrt, vector3(0, 0, 0));
181 0 : return 1;
182 : }
183 :
184 : number avDist = 0;
185 : vector3 parentCenter (0, 0, 0);
186 :
187 0 : for(size_t i = 0; i < numVrts; ++i){
188 0 : vector3 p = pos(vrts[i]);
189 0 : avDist += DistancePointToRay(p, m_center, m_axis);
190 : parentCenter += p;
191 : }
192 :
193 0 : avDist /= (number)numVrts;
194 0 : VecScale(parentCenter, parentCenter, 1. / (number)numVrts);
195 :
196 : vector3 proj, v;
197 0 : ProjectPointToRay(proj, parentCenter, m_center, m_axis);
198 : VecSubtract(v, parentCenter, proj);
199 : number len = VecLength(v);
200 0 : if(len > SMALL * avDist){ // if avDist is very small, len may be small, too
201 0 : VecScale(v, v, avDist/len);
202 : proj += v;
203 : set_pos(vrt, proj);
204 : }
205 : else
206 : set_pos(vrt, parentCenter);
207 :
208 0 : if (vertex_at_soma_surf(vrt, parent))
209 0 : project_to_soma_surface(vrt, parent);
210 :
211 : return 1;
212 : }
213 :
214 : template <class TElem>
215 0 : bool vertex_at_soma_surf(Vertex* vrt, TElem* parent) {
216 : size_t numSomaVerts = 0;
217 0 : size_t nVrt = parent->num_vertices();
218 0 : for (size_t i = 0; i < nVrt; ++i) {
219 0 : if(m_aaSurfParams[parent->vertex(i)].axial == -1.0) {
220 0 : numSomaVerts++;
221 : }
222 : }
223 0 : return numSomaVerts >= 2;
224 : }
225 :
226 : template <class TElem>
227 0 : void project_to_soma_surface(Vertex* vrt, TElem* parent) {
228 : /// old vertex position
229 : vector3 v0 = this->pos(vrt);
230 : /// 1. P = (x',y',z') = (x - x0, y - y0, z - z0)
231 : vector3 P;
232 : VecSubtract(P, v0, m_soma);
233 : /// 2. |P| = sqrt(x'^2 + y'^2 + z'^2)
234 : number Plength = VecLength(P);
235 : /// 3. Q = (radius/|P|)*P
236 : vector3 Q = P;
237 0 : VecScale(Q, P, m_somaRad / Plength);
238 : /// 4. R = Q + (x0,y0,z0)
239 : vector3 R;
240 : VecAdd(R, Q, m_soma);
241 : /// 5. Set new position
242 : set_pos(vrt, R);
243 : /// 6. set the surface parameters for the new vertex
244 0 : m_aaSurfParams[vrt].axial = -1;
245 0 : }
246 :
247 :
248 : friend class boost::serialization::access;
249 :
250 : template <class Archive>
251 0 : void serialize( Archive& ar, const unsigned int version)
252 : {
253 0 : ar & make_nvp("center", m_center);
254 0 : ar & make_nvp("axis", m_axis);
255 0 : ar & make_nvp("soma", m_soma);;
256 0 : ar & make_nvp("somaRad", m_somaRad);;
257 : UG_EMPTY_BASE_CLASS_SERIALIZATION(SomaProjector, RefinementProjector);
258 0 : }
259 :
260 : vector3 m_center;
261 : vector3 m_axis;
262 : vector3 m_soma;
263 : number m_somaRad;
264 : };
265 :
266 : }// end of namespace
267 :
268 : #endif //__H__UG_cylinder_soma_projector
|
