| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 // This example is provided by the Geant4-DNA
27 // Any report or published results obtained us
28 // shall cite the following Geant4-DNA collabo
29 // Med. Phys. 37 (2010) 4692-4708
30 // and papers
31 // M. Batmunkh et al. J Radiat Res Appl Sci 8
32 // O. Belov et al. Physica Medica 32 (2016) 15
33 // The Geant4-DNA web site is available at htt
34 //
35 // -------------------------------------------
36 // November 2016
37 // -------------------------------------------
38 //
39 //
40 /// \file NeuronLoadDataFile.cc
41 /// \brief Implementation of the NeuronLoadDat
42
43 #include "NeuronLoadDataFile.hh"
44
45 #include "CommandLineParser.hh"
46
47 #include "G4Colour.hh"
48 #include "G4LogicalVolume.hh"
49 #include "G4PhysicalConstants.hh"
50 #include "G4RotationMatrix.hh"
51 #include "G4SystemOfUnits.hh"
52 #include "G4UImanager.hh"
53 #include "G4UnitsTable.hh"
54 #include "G4VPhysicalVolume.hh"
55 #include "G4VisAttributes.hh"
56 #include "G4ios.hh"
57 #include "globals.hh"
58
59 #include <algorithm>
60 #include <fstream>
61 #include <iostream>
62 #include <limits>
63 #include <sstream>
64
65 using namespace std;
66 using namespace G4DNAPARSER;
67
68 //....oooOO0OOooo........oooOO0OOooo........oo
69
70 NeuronLoadDataFile::NeuronLoadDataFile()
71 {
72 Command* commandLine(nullptr);
73
74 // 1. Load single neuron morphology and obta
75 // Default SWC file name of neuron
76 fNeuronFileNameSWC = "GranuleCell-Nr2.CNG.sw
77
78 // 2. Load neural network and obtain paramet
79 // Default prepared data filename of neural
80 // Small network of 10 pyramidal neurons wit
81 fNeuronFileNameDATA = "NeuralNETWORK.dat";
82
83 // Load/change SWC or DAT as "CommandLinePar
84 if ((commandLine = CommandLineParser::GetPar
85 fNeuronFileNameSWC = G4String(commandLine-
86 }
87 if ((commandLine = CommandLineParser::GetPar
88 auto ss = G4String(commandLine->GetOption(
89 if ("" != ss) {
90 fNeuronFileNameDATA = ss;
91 }
92 NeuralNetworkDATAfile(fNeuronFileNameDATA)
93 }
94 else {
95 SingleNeuronSWCfile(fNeuronFileNameSWC);
96 }
97 }
98
99 //....oooOO0OOooo........oooOO0OOooo........oo
100
101 void NeuronLoadDataFile::SingleNeuronSWCfile(c
102 {
103 // -----------
104 // 12 November 2012 - code created
105 // -----------------------------------------
106 // November 2012: First model of neuron[*] a
107 // from Claiborne`s database[
108 // February 2013: Loading SWC file from Neur
109 // suggested by L. Bayarchime
110 // [*] http://lt-jds.jinr.ru/record/62124/fi
111 // [**] http://www.utsa.edu/claibornelab/
112 // [***] http://neuromorpho.org
113 // -----------------------------------------
114
115 G4String sLine = "";
116 std::ifstream infile;
117 infile.open(filename.c_str());
118 if (!infile.is_open()) {
119 G4ExceptionDescription ed;
120 ed << "Datafile " << filename << " is not
121 G4Exception("NeuronLoadDataFile::SingleNeu
122 "Check file path");
123 return;
124 }
125 G4int nrows, nlines;
126 nrows = 0;
127 nlines = 0;
128 for (;;) {
129 getline(infile, sLine);
130 if (infile.eof()) {
131 break;
132 }
133 ++nrows;
134 }
135 infile.close();
136
137 G4cout << "NeuronLoadDataFile::SingleNeuronS
138 << filename << G4endl;
139
140 infile.open(filename.c_str());
141
142 G4double TotVolSoma, TotVolDend, TotVolAxon,
143 TotVolSoma = TotVolDend = TotVolAxon = TotVo
144 G4double TotSurfSoma, TotSurfDend, TotSurfAx
145 TotSurfSoma = TotSurfDend = TotSurfAxon = To
146 G4int nNcomp; // current index of neuronal
147 G4int typeNcomp; // type of neuron structur
148 G4double x, y, z; // cartesian coordinates
149 G4double radius; // radius of each compartm
150 G4int pNcomp; // linked compartment, indica
151 G4double minX, minY, minZ;
152 G4double maxX, maxY, maxZ;
153 G4double maxRad = -1e+09;
154 minX = minY = minZ = 1e+09;
155 maxX = maxY = maxZ = -1e+09;
156 G4double density = 1.0 * (g / cm3); // wate
157 G4double Piconst = (4.0 / 3.0) * pi;
158
159 fMassSomacomp.resize(nrows, 0);
160 fPosSomacomp.resize(nrows);
161 fRadSomacomp.resize(nrows, 0);
162 std::vector<G4ThreeVector> PosDendcomp(nrows
163 fRadDendcomp.resize(nrows, 0);
164 fHeightDendcomp.resize(nrows, 0);
165 fMassDendcomp.resize(nrows, 0);
166 fDistADendSoma.resize(nrows, 0);
167 fDistBDendSoma.resize(nrows, 0);
168 fPosDendcomp.resize(nrows);
169 fRotDendcomp.resize(nrows);
170 std::vector<G4ThreeVector> PosAxoncomp(nrows
171 fRadAxoncomp.resize(nrows, 0);
172 fHeightAxoncomp.resize(nrows, 0);
173 fMassAxoncomp.resize(nrows, 0);
174 fDistAxonsoma.resize(nrows, 0);
175 fPosAxoncomp.resize(nrows);
176 fRotAxoncomp.resize(nrows);
177 fMassSpinecomp.resize(nrows, 0);
178 fPosSpinecomp.resize(nrows);
179 fRadSpinecomp.resize(nrows, 0);
180 fRadNeuroncomp.resize(nrows, 0);
181 fHeightNeuroncomp.resize(nrows, 0);
182 fDistNeuronsoma.resize(nrows, 0);
183 fPosNeuroncomp.resize(nrows);
184 fRotNeuroncomp.resize(nrows);
185 fPosNeuroncomp.resize(nrows);
186 fRadNeuroncomp.resize(nrows, 0);
187 fTypeN.resize(nrows, 0);
188 G4ThreeVector base;
189
190 // to read datafile containing numbers, alph
191 for (;;) {
192 getline(infile, sLine);
193 if (infile.eof()) {
194 break;
195 }
196 if ("#" == sLine.substr(0, 1)) {
197 continue;
198 };
199
200 std::istringstream form(sLine);
201 form >> nNcomp >> typeNcomp >> x >> y >> z
202 /*
203 G4cout << "NeuronLoadDataFile::SingleNeuronS
204 << typeNcomp << " nNcomp=" << nNcomp
205 << fnbSomacomp << " N2=" << fnbAxonco
206 */
207 // =======================================
208 // to find the largest and the smallest va
209 // for parameters of bounding slice, spher
210 if (minX > x) minX = x;
211 if (minY > y) minY = y;
212 if (minZ > z) minZ = z;
213 if (maxX < x) maxX = x;
214 if (maxY < y) maxY = y;
215 if (maxZ < z) maxZ = z;
216 // max diameter of compartments
217 if (maxRad < radius) maxRad = radius;
218
219 // =======================================
220 // Soma compartments represented as Sphere
221 if (typeNcomp == 1) {
222 // Sphere volume and surface area
223 G4double VolSomacomp = Piconst * pow(rad
224 TotVolSoma = TotVolSoma + VolSomacomp;
225 G4double SurSomacomp = 3. * Piconst * po
226 TotSurfSoma = TotSurfSoma + SurSomacomp;
227 fMassSomacomp[fnbSomacomp] = density * V
228 fMassSomaTot = fMassSomaTot + fMassSomac
229 G4ThreeVector vSoma(x, y, z);
230 fPosSomacomp[fnbSomacomp] = vSoma;
231 fRadSomacomp[fnbSomacomp] = radius;
232 if (0 == fnbSomacomp) {
233 base = G4ThreeVector(fRadSomacomp[0],
234 }
235 ++fnbSomacomp;
236 }
237
238 // =======================================
239 // Apical and basal dendritic compartments
240 if (typeNcomp == 3 || typeNcomp == 4) {
241 G4ThreeVector vDend(x, y, z);
242 // Position and Radius of compartments
243 PosDendcomp[fnbDendritecomp] = vDend;
244 fRadDendcomp[fnbDendritecomp] = radius;
245 // To join two tracing points along the
246 // To calculate length, center and rotat
247 fPosDendcomp[fnbDendritecomp] = vDend;
248 // delta of position A and position B of
249 G4ThreeVector dend;
250 // primary dendritic branch should be co
251 if (0 == fnbDendritecomp) {
252 dend = PosDendcomp[fnbDendritecomp] -
253 }
254 else {
255 dend = PosDendcomp[fnbDendritecomp] -
256 }
257 // Height of compartment
258 G4double lengthDendcomp = dend.mag();
259 fHeightDendcomp[fnbDendritecomp] = lengt
260
261 // Distance from Soma
262 G4ThreeVector dendDis = fPosSomacomp[0]
263 if (typeNcomp == 3) fDistADendSoma[fnbDe
264 if (typeNcomp == 4) fDistBDendSoma[fnbDe
265
266 // Cylinder volume and surface area
267 G4double VolDendcomp = pi * pow(radius *
268 TotVolDend = TotVolDend + VolDendcomp;
269 G4double SurDendcomp = 2. * pi * radius
270 TotSurfDend = TotSurfDend + SurDendcomp;
271 fMassDendcomp[fnbDendritecomp] = density
272 fMassDendTot = fMassDendTot + fMassDendc
273
274 dend = dend.unit();
275
276 // Euler angles of each compartment
277 G4double theta_eulerDend = dend.theta();
278 G4double phi_eulerDend = dend.phi();
279 G4double psi_eulerDend = 0;
280
281 // Rotation Matrix, Euler constructor bu
282 G4RotationMatrix rotmDendInv =
283 G4RotationMatrix(phi_eulerDend + pi /
284 fRotDendcomp[fnbDendritecomp] = rotmDend
285 ++fnbDendritecomp;
286 }
287
288 // =======================================
289 // Axon compartments represented as cylind
290 if (typeNcomp == 2 || typeNcomp == 7) {
291 G4ThreeVector vAxon(x, y, z);
292 // Position and Radius of compartments
293 PosAxoncomp[fnbAxoncomp] = vAxon;
294 fRadAxoncomp[fnbAxoncomp] = radius;
295 // To join two tracing points in loaded
296 // To calculate length, center and rotat
297
298 // delta of position A and position B of
299 G4ThreeVector Axon;
300 // primary axon point should be connect
301 if (0 == fnbAxoncomp) {
302 Axon = PosAxoncomp[fnbAxoncomp] - fPos
303 }
304 else {
305 Axon = PosAxoncomp[fnbAxoncomp] - PosA
306 }
307 G4double lengthAxoncomp = Axon.mag();
308 // Height of compartment
309 fHeightAxoncomp[fnbAxoncomp] = lengthAxo
310
311 // Distance from Soma
312 G4ThreeVector AxonDis = fPosSomacomp[0]
313 fDistAxonsoma[fnbAxoncomp] = AxonDis.mag
314
315 // Cylinder volume and surface area
316 G4double VolAxoncomp = pi * pow(radius *
317 TotVolAxon = TotVolAxon + VolAxoncomp;
318 G4double SurAxoncomp = 2. * pi * radius
319 TotSurfAxon = TotSurfAxon + SurAxoncomp;
320 fMassAxoncomp[fnbAxoncomp] = density * V
321 fMassAxonTot += fMassAxoncomp[fnbAxoncom
322 Axon = Axon.unit();
323
324 // Euler angles of each compartment
325 G4double theta_eulerAxon = Axon.theta();
326 G4double phi_eulerAxon = Axon.phi();
327 G4double psi_eulerAxon = 0;
328
329 // Rotation Matrix, Euler constructor bu
330 G4RotationMatrix rotmAxonInv =
331 G4RotationMatrix(phi_eulerAxon + pi /
332 G4RotationMatrix rotmAxon = rotmAxonInv.
333 fRotAxoncomp[fnbAxoncomp] = rotmAxon;
334 ++fnbAxoncomp;
335 }
336 // =======================================
337 // checking additional types
338 if (typeNcomp != 1 && typeNcomp != 2 && ty
339 G4cout << " Additional types:--> " << t
340 }
341
342 // If tracing points including spines, use
343 // including stubby, mushroom, thin, long
344 // branched with heads and necks!
345
346 if (typeNcomp == 5) {
347 // Sphere volume and surface area
348 G4double VolSpinecomp = Piconst * pow(ra
349 TotVolSpine = TotVolSpine + VolSpinecomp
350 G4double SurSpinecomp = 3. * Piconst * p
351 TotSurfSpine = TotSurfSpine + SurSpineco
352 fMassSpinecomp[fnbSpinecomp] = density *
353 fMassSpineTot = fMassSpineTot + fMassSpi
354 // OR
355 // Ellipsoid volume and Approximate for
356 // ...
357 G4ThreeVector vSpine(x, y, z);
358 fPosSpinecomp[fnbSpinecomp] = vSpine;
359 fRadSpinecomp[fnbSpinecomp] = radius;
360 ++fnbSpinecomp;
361 }
362 ++nlines;
363 }
364 infile.close();
365 // =========================================
366
367 fnbNeuroncomp = nlines;
368 G4cout << " Total number of compartments int
369 G4cout << G4endl;
370
371 // to calculate SHIFT value for neuron trans
372 fshiftX = (minX + maxX) / 2.;
373 fshiftY = (minY + maxY) / 2.;
374 fshiftZ = (minZ + maxZ) / 2.;
375
376 // width, height, depth of bounding slice vo
377 fwidthB = std::fabs(minX - maxX) + maxRad;
378 fheightB = std::fabs(minY - maxY) + maxRad;
379 fdepthB = std::fabs(minZ - maxZ) + maxRad;
380
381 // diagonal length of bounding slice, that g
382 // for particle direction and fluence!
383 fdiagnlLength = std::sqrt(fwidthB * fwidthB
384
385 fTotVolNeuron = TotVolSoma + TotVolDend + To
386 fTotSurfNeuron = TotSurfSoma + TotSurfDend +
387 fTotMassNeuron = fMassSomaTot + fMassDendTot
388
389 fTotVolSlice = fwidthB * um * fheightB * um
390 fTotSurfSlice =
391 2 * (fwidthB * um * fheightB * um + fheigh
392 fTotMassSlice = 1.0 * (g / cm3) * fTotVolSli
393
394 fTotVolMedium = Piconst * pow(fdiagnlLength
395 fTotSurfMedium = 3. * Piconst * pow(fdiagnlL
396 fTotMassMedium = 1.0 * (g / cm3) * fTotVolMe
397 }
398
399 //....oooOO0OOooo........oooOO0OOooo........oo
400
401 // Load prepared data file of neural network w
402 void NeuronLoadDataFile::NeuralNetworkDATAfile
403 {
404 G4String sLine = "";
405 std::ifstream infile;
406 infile.open(filename.c_str());
407 if (!infile.is_open()) {
408 G4ExceptionDescription ed;
409 ed << "Datafile " << filename << " is not
410 G4Exception("NeuronLoadDataFile::NeuralNet
411 "Check file path");
412 return;
413 }
414 G4cout << "NeuronLoadDataFile::NeuralNetwork
415
416 G4int nlines, nbSoma, nbDendrite;
417 nlines = 0;
418 fnbSomacomp = 0; // total number of compart
419 fnbDendritecomp = 0; // total number of com
420 fnbAxoncomp = 0; // total number of compart
421 fnbSpinecomp = 0; // total number of compar
422 G4double TotVolSoma, TotVolDend, TotVolAxon;
423 TotVolSoma = TotVolDend = TotVolAxon = 0.;
424 G4double TotSurfSoma, TotSurfDend, TotSurfAx
425 TotSurfSoma = TotSurfDend = TotSurfAxon = 0.
426 G4int typeNcomp; // types of structure: som
427 G4double x1, y1, z1, x2, y2, z2; // cartesi
428 G4double radius; // radius of each compartm
429 G4double height; // height of each compartm
430 G4double maxRad = -1e+09;
431 G4double density = 1.0 * (g / cm3); // wate
432 G4double Piconst = (4.0 / 3.0) * pi;
433
434 for (;;) {
435 getline(infile, sLine);
436 if (infile.eof()) {
437 break;
438 }
439 std::istringstream form(sLine);
440 if (nlines == 0) {
441 // to read total number of compartments
442 form >> fnbNeuroncomp >> nbSoma >> nbDen
443 fMassSomacomp.resize(nbSoma, 0);
444 fPosSomacomp.resize(nbSoma);
445 fRadSomacomp.resize(nbSoma, 0);
446 fRadDendcomp.resize(nbDendrite, 0);
447 fHeightDendcomp.resize(nbDendrite, 0);
448 fMassDendcomp.resize(nbDendrite, 0);
449 fDistADendSoma.resize(nbDendrite, 0);
450 fDistBDendSoma.resize(nbDendrite, 0);
451 fPosDendcomp.resize(nbDendrite);
452 fRotDendcomp.resize(nbDendrite);
453 }
454 // =======================================
455 // Soma compartments represented as Sphere
456 if (nlines > 0 && nlines <= nbSoma) {
457 form >> typeNcomp >> x1 >> y1 >> z1 >> r
458 if (typeNcomp != 1) break;
459 // max diameter of compartments
460 if (maxRad < radius) maxRad = radius;
461 // Sphere volume and surface area
462 G4double VolSomacomp = Piconst * pow(rad
463 TotVolSoma = TotVolSoma + VolSomacomp;
464 G4double SurSomacomp = 3. * Piconst * po
465 TotSurfSoma = TotSurfSoma + SurSomacomp;
466 fMassSomacomp[fnbSomacomp] = density * V
467 fMassSomaTot = fMassSomaTot + fMassSomac
468
469 G4ThreeVector vSoma(x1, y1, z1);
470 fPosSomacomp[fnbSomacomp] = vSoma;
471 fRadSomacomp[fnbSomacomp] = radius;
472 ++fnbSomacomp;
473 }
474 // =======================================
475 // Apical and basal dendritic compartments
476 if (nlines > nbSoma && nlines <= fnbNeuron
477 form >> typeNcomp >> x1 >> y1 >> z1 >> x
478 if (typeNcomp != 3) break; // || typeNc
479
480 // To calculate length, center and rotat
481 // Center-position of each cylinder
482 G4double Dendxx = x1 + x2;
483 G4double Dendyy = y1 + y2;
484 G4double Dendzz = z1 + z2;
485 G4ThreeVector translmDend = G4ThreeVecto
486 fPosDendcomp[fnbDendritecomp] = translmD
487 fRadDendcomp[fnbDendritecomp] = radius;
488 G4double lengthDendcomp = height;
489 // Height of compartment
490 fHeightDendcomp[fnbDendritecomp] = lengt
491 // Distance from Soma
492
493 // Cylinder volume and surface area
494 G4double VolDendcomp = pi * pow(radius *
495 TotVolDend = TotVolDend + VolDendcomp;
496 G4double SurDendcomp = 2. * pi * radius
497 TotSurfDend = TotSurfDend + SurDendcomp;
498 fMassDendcomp[fnbDendritecomp] = density
499 fMassDendTot = fMassDendTot + fMassDendc
500
501 G4double Dendx = x1 - x2;
502 G4double Dendy = y1 - y2;
503 G4double Dendz = z1 - z2;
504 Dendx = Dendx / lengthDendcomp;
505 Dendy = Dendy / lengthDendcomp;
506 Dendz = Dendz / lengthDendcomp;
507
508 // Euler angles of each compartment
509 G4ThreeVector directionDend = G4ThreeVec
510 G4double theta_eulerDend = directionDend
511 G4double phi_eulerDend = directionDend.p
512 G4double psi_eulerDend = 0;
513
514 // Rotation Matrix, Euler constructor bu
515 G4RotationMatrix rotmDendInv =
516 G4RotationMatrix(phi_eulerDend + pi /
517 G4RotationMatrix rotmDend = rotmDendInv.
518
519 fRotDendcomp[fnbDendritecomp] = rotmDend
520 ++fnbDendritecomp;
521 }
522 ++nlines;
523 }
524
525 // =========================================
526
527 G4cout << " Total number of compartments int
528
529 // to calculate SHIFT value for neuron trans
530 fshiftX = 0.; //(minX + maxX)/2. ;
531 fshiftY = 0.; //(minY + maxY)/2. ;
532 fshiftZ = 0.; //(minZ + maxZ)/2. ;
533
534 // width, height, depth of bounding slice vo
535 fwidthB = 640.;
536 fheightB = 280.;
537 fdepthB = 25.;
538 // diagonal length of bounding slice, that g
539 // for particle direction and fluence!
540 fdiagnlLength = std::sqrt(fwidthB * fwidthB
541
542 fTotVolNeuron = TotVolSoma + TotVolDend + To
543 fTotSurfNeuron = TotSurfSoma + TotSurfDend +
544 fTotMassNeuron = fMassSomaTot + fMassDendTot
545
546 fTotVolSlice = fwidthB * um * fheightB * um
547 fTotSurfSlice =
548 2 * (fwidthB * um * fheightB * um + fheigh
549 fTotMassSlice = 1.0 * (g / cm3) * fTotVolSli
550
551 fTotVolMedium = Piconst * pow(fdiagnlLength
552 fTotSurfMedium = 3. * Piconst * pow(fdiagnlL
553 fTotMassMedium = 1.0 * (g / cm3) * fTotVolMe
554
555 infile.close();
556 }
557
558 //....oooOO0OOooo........oooOO0OOooo........oo
559
560 void NeuronLoadDataFile::ComputeTransformation
561 {
562 // to calculate Euler angles from Rotation M
563 //
564 G4RotationMatrix rotmNeuron = G4RotationMatr
565 G4double cosX = std::sqrt(rotmNeuron.xx() *
566 G4double euX, euY, euZ;
567 if (cosX > 16 * FLT_EPSILON) {
568 euX = std::atan2(rotmNeuron.zy(), rotmNeur
569 euY = std::atan2(-rotmNeuron.zx(), cosX);
570 euZ = std::atan2(rotmNeuron.yx(), rotmNeur
571 }
572 else {
573 euX = std::atan2(-rotmNeuron.yz(), rotmNeu
574 euY = std::atan2(-rotmNeuron.zx(), cosX);
575 euZ = 0.;
576 }
577 G4RotationMatrix* rot = new G4RotationMatrix
578 rot->rotateX(euX);
579 rot->rotateY(euY);
580 rot->rotateZ(euZ);
581
582 physVol->SetRotation(rot);
583
584 // shift of cylinder compartments
585 G4ThreeVector originNeuron((fPosNeuroncomp[c
586 (fPosNeuroncomp[c
587 (fPosNeuroncomp[c
588 physVol->SetTranslation(originNeuron);
589 }
590
591 //....oooOO0OOooo........oooOO0OOooo........oo
592
593 void NeuronLoadDataFile::ComputeDimensions(G4T
594 con
595 {
596 fcylinderComp.SetInnerRadius(0 * um);
597 fcylinderComp.SetOuterRadius(fRadNeuroncomp[
598 fcylinderComp.SetZHalfLength(fHeightNeuronco
599 fcylinderComp.SetStartPhiAngle(0. * deg);
600 fcylinderComp.SetDeltaPhiAngle(360. * deg);
601 }
602
603 //....oooOO0OOooo........oooOO0OOooo........oo
604