| 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 //
27 //
28 //
29 // File I/O manager class for writing or readi
30 // distribution and some event information
31 //
32 // Created: Mar. 31, 2009 Akinori Kimura : r
33 //
34 // Akinori Kimur
35 // gMocren home
36 // http://geant4
37 //
38 //
39 #include "G4GMocrenIO.hh"
40 #include <iostream>
41 #include <ctime>
42 #include <sstream>
43 #include <iomanip>
44 #include <cstdlib>
45 #include <cstring>
46
47 #include "globals.hh"
48 #include "G4VisManager.hh"
49
50 #if defined(_WIN32)
51 #define LITTLE_ENDIAN 1234
52 #define BYTE_ORDER LITTLE_ENDIAN
53 #endif
54
55 const int DOSERANGE = 25000;
56
57 //----- GMocrenDataPrimitive class in the GMoc
58 template <typename T>
59 GMocrenDataPrimitive<T>::GMocrenDataPrimitive
60 clear();
61 }
62 template <typename T>
63 GMocrenDataPrimitive<T>::~GMocrenDataPrimitive
64 /*
65 std::vector<short *>::iterator itr = image
66 for(; itr != image.end(); itr++) {
67 delete [] *itr;
68 }
69 */
70 }
71
72 template <typename T> GMocrenDataPrimitive<T>
73 GMocrenDataPrimitive<T>::operator = (const GMo
74 if (this == &_right) return *this;
75 for(int i = 0; i < 3; i++) {
76 kSize[i] = _right.kSize[i];
77 kCenter[i] = _right.kCenter[i];
78 }
79 kScale = _right.kScale;
80 for(int i = 0; i < 2; i++) kMinmax[i] = _rig
81 int num = kSize[0]*kSize[1];
82 kImage.clear();
83 for(int z = 0; z < kSize[2]; z++) {
84 T * img = new T[num];
85 for(int i = 0; i < num; i++) img[i] =_righ
86 kImage.push_back(img);
87 }
88 return *this;
89 }
90
91 template <typename T> GMocrenDataPrimitive<T>
92 GMocrenDataPrimitive<T>::operator + (const GMo
93
94 GMocrenDataPrimitive<T> rprim;
95 bool stat = true;
96 for(int i = 0; i < 3; i++) {
97 if(kSize[i] != _right.kSize[i]) stat = fal
98 if(kCenter[i] != _right.kCenter[i]) stat =
99 }
100 if(!stat) {
101 if (G4VisManager::GetVerbosity() >= G4VisM
102 G4cout << "Warning: operator + "
103 << " Cannot do the operator +"
104 << G4endl;
105 return *this;
106 }
107
108 rprim.setSize(kSize);
109 rprim.setCenterPosition(kCenter);
110
111 T mms[2] = {9e100,-9e100};
112 //if(mms[0] > _right.minmax[0]) mms[0] = _ri
113 //if(mms[1] < _right.minmax[1]) mms[1] = _ri
114
115 int num = kSize[0]*kSize[1];
116 for(int z = 0; z < kSize[2]; z++) {
117 T * img = new T[num];
118 for(int xy = 0; xy < num; xy++) {
119 img[xy] = kImage[z][xy] + _right.kImage[
120 if(mms[0] > img[xy]) mms[0] = img[xy];
121 if(mms[1] < img[xy]) mms[1] = img[xy];
122 }
123 rprim.addImage(img);
124 }
125 rprim.setMinMax(mms);
126
127 T scl = mms[1]/DOSERANGE;
128 rprim.setScale(scl);
129
130 return rprim;
131 }
132
133 template <typename T> GMocrenDataPrimitive<T>
134 GMocrenDataPrimitive<T>::operator += (const GM
135
136 bool stat = true;
137 for(int i = 0; i < 3; i++) {
138 if(kSize[i] != _right.kSize[i]) stat = fal
139 if(kCenter[i] != _right.kCenter[i]) stat =
140 }
141 if(!stat) {
142 if (G4VisManager::GetVerbosity() >= G4VisM
143 G4cout << "Warning: operator += " << G4e
144 << " Cannot do the operator +="
145 << G4endl;
146 return *this;
147 }
148
149 if(kMinmax[0] > _right.kMinmax[0]) kMinmax[0
150 if(kMinmax[1] < _right.kMinmax[1]) kMinmax[1
151
152 int num = kSize[0]*kSize[1];
153 for(int z = 0; z < kSize[2]; z++) {
154 for(int xy = 0; xy < num; xy++) {
155 kImage[z][xy] += _right.kImage[z][xy];
156 if(kMinmax[0] > kImage[z][xy]) kMinmax[0
157 if(kMinmax[1] < kImage[z][xy]) kMinmax[1
158 }
159 }
160
161 kScale = kMinmax[1]/DOSERANGE;
162
163 return *this;
164 }
165
166 template <typename T>
167 void GMocrenDataPrimitive<T>::clear() {
168 for(int i = 0; i < 3; i++) {
169 kSize[i] = 0;
170 kCenter[i] = 0.;
171 }
172 kScale = 1.;
173 kMinmax[0] = (T)32109;
174 kMinmax[1] = (T)-32109;
175
176 clearImage();
177 }
178 template <typename T>
179 void GMocrenDataPrimitive<T>::clearImage() {
180 typename std::vector<T *>::iterator itr;
181 for(itr = kImage.begin(); itr != kImage.end(
182 delete [] *itr;
183 }
184 kImage.clear();
185 }
186 template <typename T>
187 void GMocrenDataPrimitive<T>::setSize(int _siz
188 for(int i = 0; i < 3; i++) kSize[i] = _size[
189 }
190 template <typename T>
191 void GMocrenDataPrimitive<T>::getSize(int _siz
192 for(int i = 0; i < 3; i++) _size[i] = kSize[
193 }
194 template <typename T>
195 void GMocrenDataPrimitive<T>::setScale(double
196 kScale = _scale;
197 }
198 template <typename T>
199 double GMocrenDataPrimitive<T>::getScale() {
200 return kScale;
201 }
202 template <typename T>
203 void GMocrenDataPrimitive<T>::setMinMax(T _min
204 for(int i = 0; i < 2; i++) kMinmax[i] = _min
205 }
206 template <typename T>
207 void GMocrenDataPrimitive<T>::getMinMax(T _min
208 for(int i = 0; i < 2; i++) _minmax[i] = kMin
209
210 }
211 template <typename T>
212 void GMocrenDataPrimitive<T>::setImage(std::ve
213 kImage = _image;
214 }
215 template <typename T>
216 void GMocrenDataPrimitive<T>::addImage(T * _im
217 kImage.push_back(_image);
218 }
219 template <typename T>
220 std::vector<T *> & GMocrenDataPrimitive<T>::ge
221 return kImage;
222 }
223 template <typename T>
224 T * GMocrenDataPrimitive<T>::getImage(int _z)
225 if(_z >= (int)kImage.size()) return 0;
226 return kImage[_z];
227 }
228 template <typename T>
229 void GMocrenDataPrimitive<T>::setCenterPositio
230 for(int i = 0; i < 3; i++) kCenter[i] = _cen
231 }
232 template <typename T>
233 void GMocrenDataPrimitive<T>::getCenterPositio
234 for(int i = 0; i < 3; i++) _center[i] = kCen
235 }
236 template <typename T>
237 void GMocrenDataPrimitive<T>::setName(std::str
238 kDataName = _name;
239 }
240 template <typename T>
241 std::string GMocrenDataPrimitive<T>::getName()
242 return kDataName;
243 }
244
245
246
247
248
249 GMocrenTrack::GMocrenTrack() {
250 kTrack.clear();
251 for(int i = 0; i < 3; i++) kColor[i] = 0;
252 }
253
254 void GMocrenTrack::addStep(float _startx, floa
255 float _endx, float _endy, float _endz
256 struct Step step;
257 step.startPoint[0] = _startx;
258 step.startPoint[1] = _starty;
259 step.startPoint[2] = _startz;
260 step.endPoint[0] = _endx;
261 step.endPoint[1] = _endy;
262 step.endPoint[2] = _endz;
263 kTrack.push_back(step);
264 }
265 void GMocrenTrack::getStep(float & _startx, fl
266 float & _endx, float & _endy, float &
267 int _num) {
268 if(_num >= (int)kTrack.size()) {
269 if (G4VisManager::GetVerbosity() >= G4VisM
270 G4cout << "GMocrenTrack::getStep(...) Er
271 << "invalid step # : " << _num << G4end
272 return;
273 }
274
275 _startx = kTrack[_num].startPoint[0];
276 _starty = kTrack[_num].startPoint[1];
277 _startz = kTrack[_num].startPoint[2];
278 _endx = kTrack[_num].endPoint[0];
279 _endy = kTrack[_num].endPoint[1];
280 _endz = kTrack[_num].endPoint[2];
281 }
282 void GMocrenTrack::translate(std::vector<float
283 std::vector<struct Step>::iterator itr = kTr
284 for(; itr != kTrack.end(); itr++) {
285 for(int i = 0; i < 3; i++ ) {
286 itr->startPoint[i] += _translate[i];
287 itr->endPoint[i] += _translate[i];
288 }
289 }
290 }
291
292
293
294
295
296
297
298
299
300 GMocrenDetector::GMocrenDetector() {
301 kDetector.clear();
302 for(int i = 0; i < 3; i++) kColor[i] = 0;
303 }
304
305 void GMocrenDetector::addEdge(float _startx, f
306 float _endx, float _endy, float _e
307 struct Edge edge;
308 edge.startPoint[0] = _startx;
309 edge.startPoint[1] = _starty;
310 edge.startPoint[2] = _startz;
311 edge.endPoint[0] = _endx;
312 edge.endPoint[1] = _endy;
313 edge.endPoint[2] = _endz;
314 kDetector.push_back(edge);
315 }
316 void GMocrenDetector::getEdge(float & _startx,
317 float & _endx, float & _endy, float &
318 int _num) {
319 if(_num >= (int)kDetector.size()) {
320 if (G4VisManager::GetVerbosity() >= G4VisM
321 G4cout << "GMocrenDetector::getEdge(...)
322 << "invalid edge # : " << _num << G4end
323 return;
324 }
325
326 _startx = kDetector[_num].startPoint[0];
327 _starty = kDetector[_num].startPoint[1];
328 _startz = kDetector[_num].startPoint[2];
329 _endx = kDetector[_num].endPoint[0];
330 _endy = kDetector[_num].endPoint[1];
331 _endz = kDetector[_num].endPoint[2];
332 }
333 void GMocrenDetector::translate(std::vector<fl
334 std::vector<struct Edge>::iterator itr = kDe
335 for(; itr != kDetector.end(); itr++) {
336 for(int i = 0; i < 3; i++) {
337 itr->startPoint[i] += _translate[i];
338 itr->endPoint[i] += _translate[i];
339 }
340 }
341 }
342
343
344
345
346
347
348
349
350
351 // file information
352 std::string G4GMocrenIO::kId;
353 std::string G4GMocrenIO::kVersion = "2.0.0";
354 int G4GMocrenIO::kNumberOfEvents = 0;
355 char G4GMocrenIO::kLittleEndianInput = true;
356
357 #if BYTE_ORDER == LITTLE_ENDIAN
358 char G4GMocrenIO::kLittleEndianOutput = true;
359 #else
360 char G4GMocrenIO::kLittleEndianOutput = false;
361 #endif
362 std::string G4GMocrenIO::kComment;
363 //
364 std::string G4GMocrenIO::kFileName = "dose.gdd
365
366 //
367 unsigned int G4GMocrenIO::kPointerToModalityDa
368 std::vector<unsigned int> G4GMocrenIO::kPointe
369 unsigned int G4GMocrenIO::kPointerToROIData =
370 unsigned int G4GMocrenIO::kPointerToTrackData
371 unsigned int G4GMocrenIO::kPointerToDetectorDa
372
373 // modality
374 float G4GMocrenIO::kVoxelSpacing[3] = {0., 0.,
375 class GMocrenDataPrimitive<short> G4GMocrenIO
376 std::vector<float> G4GMocrenIO::kModalityImage
377 std::string G4GMocrenIO::kModalityUnit = "g/cm
378
379 // dose
380 std::vector<class GMocrenDataPrimitive<double>
381 std::string G4GMocrenIO::kDoseUnit = "keV
382
383 // ROI
384 std::vector<class GMocrenDataPrimitive<short>
385
386 // track
387 std::vector<float *> G4GMocrenIO::kSteps;
388 std::vector<unsigned char *> G4GMocrenIO::kSte
389 std::vector<class GMocrenTrack> G4GMocrenIO::k
390
391 // detector
392 std::vector<class GMocrenDetector> G4GMocrenIO
393
394 // verbose
395 int G4GMocrenIO::kVerbose = 0;
396
397 const int IDLENGTH = 21;
398 const int VERLENGTH = 6;
399
400 // constructor
401 G4GMocrenIO::G4GMocrenIO()
402 : kTracksWillBeStored(true) {
403 ;
404 }
405
406 // destructor
407 G4GMocrenIO::~G4GMocrenIO() {
408 ;
409 }
410
411 // initialize
412 void G4GMocrenIO::initialize() {
413
414 kId.clear();
415 kVersion = "2.0.0";
416 kNumberOfEvents = 0;
417 kLittleEndianInput = true;
418 #if BYTE_ORDER == LITTLE_ENDIAN
419 kLittleEndianOutput = true;
420 #else // Big endian
421 kLittleEndianOutput = false;
422 #endif
423 kComment.clear();
424 kFileName = "dose.gdd";
425 kPointerToModalityData = 0;
426 kPointerToDoseDistData.clear();
427 kPointerToROIData = 0;
428 kPointerToTrackData = 0;
429 // modality
430 for(int i = 0; i < 3; i++) kVoxelSpacing[i]
431 kModality.clear();
432 kModalityImageDensityMap.clear();
433 kModalityUnit = "g/cm3 "; // 12 Bytes
434 // dose
435 kDose.clear();
436 kDoseUnit = "keV "; // 12 Bytes
437 // ROI
438 kRoi.clear();
439 // track
440 std::vector<float *>::iterator itr;
441 for(itr = kSteps.begin(); itr != kSteps.end(
442 kSteps.clear();
443 std::vector<unsigned char *>::iterator citr;
444 for(citr = kStepColors.begin(); citr != kSte
445 delete [] *citr;
446 kStepColors.clear();
447 kTracksWillBeStored = true;
448
449 // verbose
450 kVerbose = 0;
451 }
452
453 bool G4GMocrenIO::storeData() {
454 return storeData4();
455 }
456 //
457 bool G4GMocrenIO::storeData(char * _filename)
458 return storeData4(_filename);
459 }
460
461 bool G4GMocrenIO::storeData4() {
462
463 bool DEBUG = false;//
464
465 if(DEBUG || kVerbose > 0)
466 G4cout << ">>>>>>> store data (ver.4) <<<
467 if(DEBUG || kVerbose > 0)
468 G4cout << " " << kFileName << G4en
469
470 // output file open
471 std::ofstream ofile(kFileName.c_str(),
472 std::ios_base::out|std::ios_base::bi
473 if(DEBUG || kVerbose > 0)
474 G4cout << " file open status: " <<
475
476 // file identifier
477 ofile.write("gMocren ", 8);
478
479 // file version
480 unsigned char ver = 0x04;
481 ofile.write((char *)&ver, 1);
482
483 // endian
484 //ofile.write((char *)&kLittleEndianOutput,
485 char littleEndian = 0x01;
486 ofile.write((char *)&littleEndian, sizeof(ch
487 if(DEBUG || kVerbose > 0) {
488 //G4cout << "Endian: " << (int)kLittleEndi
489 G4cout << "Endian: " << (int)littleEndian
490 }
491
492 // for inverting the byte order
493 float ftmp[6];
494 int itmp[6];
495 short stmp[6];
496
497 // comment length (fixed size)
498 int commentLength = 1024;
499 if(kLittleEndianOutput) {
500 ofile.write((char *)&commentLength, 4);
501 } else {
502 invertByteOrder((char *)&commentLength, it
503 ofile.write((char *)itmp, 4);
504 }
505
506 // comment
507 char cmt[1025];
508 std::strncpy(cmt, kComment.c_str(), 1024);
509 cmt[1024] = '\0';
510 ofile.write(cmt, 1024);
511 if(DEBUG || kVerbose > 0) {
512 G4cout << "Data comment : "
513 << kComment << G4endl;
514 }
515
516 // voxel spacings for all images
517 if(kLittleEndianOutput) {
518 ofile.write((char *)kVoxelSpacing, 12);
519 } else {
520 for(int j = 0; j < 3; j++)
521 invertByteOrder((char *)&kVoxelSpacing[j
522 ofile.write((char *)ftmp, 12);
523 }
524 if(DEBUG || kVerbose > 0) {
525 G4cout << "Voxel spacing : ("
526 << kVoxelSpacing[0] << ", "
527 << kVoxelSpacing[1] << ", "
528 << kVoxelSpacing[2]
529 << ") mm " << G4endl;
530 }
531
532 calcPointers4();
533 if(!kTracksWillBeStored) kPointerToTrackData
534
535 // offset from file starting point to the mo
536 if(kLittleEndianOutput) {
537 ofile.write((char *)&kPointerToModalityDat
538 } else {
539 invertByteOrder((char *)&kPointerToModalit
540 ofile.write((char *)itmp, 4);
541 }
542
543 // # of dose distributions
544 //int nDoseDist = (int)pointerToDoseDistData
545 int nDoseDist = getNumDoseDist();
546 if(kLittleEndianOutput) {
547 ofile.write((char *)&nDoseDist, 4);
548 } else {
549 invertByteOrder((char *)&nDoseDist, itmp[0
550 ofile.write((char *)itmp, 4);
551 }
552
553 // offset from file starting point to the do
554 if(kLittleEndianOutput) {
555 for(int i = 0; i < nDoseDist; i++) {
556 ofile.write((char *)&kPointerToDoseDistD
557 }
558 } else {
559 for(int i = 0; i < nDoseDist; i++) {
560 invertByteOrder((char *)&kPointerToDoseD
561 ofile.write((char *)itmp, 4);
562 }
563 }
564
565 // offset from file starting point to the RO
566 if(kLittleEndianOutput) {
567 ofile.write((char *)&kPointerToROIData, 4)
568 } else {
569 invertByteOrder((char *)&kPointerToROIData
570 ofile.write((char *)itmp, 4);
571 }
572
573 // offset from file starting point to the tr
574 if(kLittleEndianOutput) {
575 ofile.write((char *)&kPointerToTrackData,
576 } else {
577 invertByteOrder((char *)&kPointerToTrackDa
578 ofile.write((char *)itmp, 4);
579 }
580
581 // offset from file starting point to the de
582 if(kLittleEndianOutput) {
583 ofile.write((char *)&kPointerToDetectorDat
584 } else {
585 invertByteOrder((char *)&kPointerToDetecto
586 ofile.write((char *)itmp, 4);
587 }
588
589 if(DEBUG || kVerbose > 0) {
590 G4cout << "Each pointer to data : "
591 << kPointerToModalityData << ", ";
592 for(int i = 0; i < nDoseDist; i++) {
593 G4cout << kPointerToDoseDistData[i] << "
594 }
595 G4cout << kPointerToROIData << ", "
596 << kPointerToTrackData << ", "
597 << kPointerToDetectorData
598 << G4endl;
599 }
600
601 //----- modality image -----//
602
603 int size[3];
604 float scale;
605 short minmax[2];
606 float fCenter[3];
607 int iCenter[3];
608 // modality image size
609 kModality.getSize(size);
610
611 if(kLittleEndianOutput) {
612 ofile.write((char *)size, 3*sizeof(int));
613 } else {
614 for(int j = 0; j < 3; j++)
615 invertByteOrder((char *)&size[j], itmp[j
616 ofile.write((char *)itmp, 12);
617 }
618
619 if(DEBUG || kVerbose > 0) {
620 G4cout << "Modality image size : ("
621 << size[0] << ", "
622 << size[1] << ", "
623 << size[2] << ")"
624 << G4endl;
625 }
626
627 // modality image max. & min.
628 kModality.getMinMax(minmax);
629 if(kLittleEndianOutput) {
630 ofile.write((char *)minmax, 4);
631 } else {
632 for(int j = 0; j < 2; j++)
633 invertByteOrder((char *)&minmax[j], stmp
634 ofile.write((char *)stmp, 4);
635 }
636
637 // modality image unit
638 char munit[13] = "g/cm3\0";
639 ofile.write((char *)munit, 12);
640
641 // modality image scale
642 scale = (float)kModality.getScale();
643 if(kLittleEndianOutput) {
644 ofile.write((char *)&scale, 4);
645 } else {
646 invertByteOrder((char *)&scale, ftmp[0]);
647 ofile.write((char *)ftmp, 4);
648 }
649 if(DEBUG || kVerbose > 0) {
650 G4cout << "Modality image min., max., scal
651 << minmax[0] << ", "
652 << minmax[1] << ", "
653 << scale << G4endl;
654 }
655
656 // modality image
657 int psize = size[0]*size[1];
658 if(DEBUG || kVerbose > 0) G4cout << "Modalit
659 for(int i = 0; i < size[2]; i++) {
660 short * image = kModality.getImage(i);
661 if(kLittleEndianOutput) {
662 ofile.write((char *)image, psize*sizeof(
663 } else {
664 for(int j = 0; j < psize; j++) {
665 invertByteOrder((char *)&image[j], stmp[0]);
666 ofile.write((char *)stmp, 2);
667 }
668 }
669
670 if(DEBUG || kVerbose > 0) G4cout << "[" <<
671 }
672 if(DEBUG || kVerbose > 0) G4cout << G4endl;
673
674 // modality desity map for CT value
675 size_t msize = minmax[1] - minmax[0]+1;
676 if(DEBUG || kVerbose > 0)
677 G4cout << "modality image : " << minmax[0]
678 float * pdmap = new float[msize];
679 for(int i = 0; i < (int)msize; i++) pdmap[i]
680
681 if(kLittleEndianOutput) {
682 ofile.write((char *)pdmap, msize*sizeof(fl
683 } else {
684 for(int j = 0; j < (int)msize; j++) {
685 invertByteOrder((char *)&pdmap[j], ftmp[
686 ofile.write((char *)ftmp, 4);
687 }
688 }
689
690 if(DEBUG || kVerbose > 0) {
691 G4cout << "density map : " << std::ends;
692 for(int i = 0; i < (int)msize; i+=50)
693 G4cout <<kModalityImageDensityMap[i] <<
694 G4cout << G4endl;
695 }
696 delete [] pdmap;
697
698
699 //----- dose distribution image -----//
700
701 if(!isDoseEmpty()) {
702
703 calcDoseDistScale();
704
705 for(int ndose = 0; ndose < nDoseDist; ndos
706 // dose distrbution image size
707 kDose[ndose].getSize(size);
708 if(kLittleEndianOutput) {
709 ofile.write((char *)size, 3*sizeof(int));
710 } else {
711 for(int j = 0; j < 3; j++)
712 invertByteOrder((char *)&size[j], itmp[j])
713 ofile.write((char *)itmp, 12);
714 }
715 if(DEBUG || kVerbose > 0) {
716 G4cout << "Dose dist. [" << ndose << "] imag
717 << size[0] << ", "
718 << size[1] << ", "
719 << size[2] << ")"
720 << G4endl;
721 }
722
723 // dose distribution max. & min.
724 getShortDoseDistMinMax(minmax, ndose);
725 if(kLittleEndianOutput) {
726 ofile.write((char *)minmax, 2*2); // sizeof(
727 } else {
728 for(int j = 0; j < 2; j++)
729 invertByteOrder((char *)&minmax[j], stmp[j
730 ofile.write((char *)stmp, 4);
731 }
732
733 // dose distribution unit
734 char cdunit[13];
735 std::strncpy(cdunit, kDoseUnit.c_str(),
736 cdunit[12] = '\0';
737 ofile.write(cdunit, 12);
738 if(DEBUG || kVerbose > 0) {
739 G4cout << "Dose dist. unit : " << kDoseUnit
740 }
741
742 // dose distribution scaling
743 double dscale;
744 dscale = getDoseDistScale(ndose);
745 scale = float(dscale);
746 if(kLittleEndianOutput) {
747 ofile.write((char *)&scale, 4);
748 } else {
749 invertByteOrder((char *)&scale, ftmp[0]);
750 ofile.write((char *)ftmp, 4);
751 }
752 if(DEBUG || kVerbose > 0) {
753 G4cout << "Dose dist. [" << ndose
754 << "] image min., max., scale : "
755 << minmax[0] << ", "
756 << minmax[1] << ", "
757 << scale << G4endl;
758 }
759
760 // dose distribution image
761 int dsize = size[0]*size[1];
762 short * dimage = new short[dsize];
763 for(int z = 0; z < size[2]; z++) {
764 getShortDoseDist(dimage, z, ndose);
765 if(kLittleEndianOutput) {
766 ofile.write((char *)dimage, dsize*2); //si
767 } else {
768 for(int j = 0; j < dsize; j++) {
769 invertByteOrder((char *)&dimage[j], stmp
770 ofile.write((char *)stmp, 2);
771 }
772 }
773
774 if(DEBUG || kVerbose > 0) {
775 for(int j = 0; j < dsize; j++) {
776 if(dimage[j] < 0)
777 G4cout << "[" << j << "," << z << "]"
778 << dimage[j] << ", ";
779 }
780 }
781 }
782 if(DEBUG || kVerbose > 0) G4cout << G4en
783 delete [] dimage;
784
785 // relative location of the dose distrib
786 // the modality image
787 getDoseDistCenterPosition(fCenter, ndose
788 for(int i = 0; i < 3; i++) iCenter[i] =
789 if(kLittleEndianOutput) {
790 ofile.write((char *)iCenter, 3*4); // 3*size
791 } else {
792 for(int j = 0; j < 3; j++)
793 invertByteOrder((char *)&iCenter[j], itmp[
794 ofile.write((char *)itmp, 12);
795 }
796 if(DEBUG || kVerbose > 0) {
797 G4cout << "Dose dist. [" << ndose
798 << "]image relative location : ("
799 << iCenter[0] << ", "
800 << iCenter[1] << ", "
801 << iCenter[2] << ")" << G4endl;
802 }
803
804 // dose distribution name
805 std::string name = getDoseDistName(ndose
806 if(name.size() == 0) name = "dose";
807 name.resize(80);
808 ofile.write((char *)name.c_str(), 80);
809 if(DEBUG || kVerbose > 0) {
810 G4cout << "Dose dist. name : " << name << G4
811 }
812
813 }
814 }
815
816 //----- ROI image -----//
817 if(!isROIEmpty()) {
818 // ROI image size
819 kRoi[0].getSize(size);
820 if(kLittleEndianOutput) {
821 ofile.write((char *)size, 3*sizeof(int))
822 } else {
823 for(int j = 0; j < 3; j++)
824 invertByteOrder((char *)&size[j], itmp[j]);
825 ofile.write((char *)itmp, 12);
826 }
827 if(DEBUG || kVerbose > 0) {
828 G4cout << "ROI image size : ("
829 << size[0] << ", "
830 << size[1] << ", "
831 << size[2] << ")"
832 << G4endl;
833 }
834
835 // ROI max. & min.
836 kRoi[0].getMinMax(minmax);
837 if(kLittleEndianOutput) {
838 ofile.write((char *)minmax, sizeof(short
839 } else {
840 for(int j = 0; j < 2; j++)
841 invertByteOrder((char *)&minmax[j], stmp[j])
842 ofile.write((char *)stmp, 4);
843 }
844
845 // ROI distribution scaling
846 scale = (float)kRoi[0].getScale();
847 if(kLittleEndianOutput) {
848 ofile.write((char *)&scale, sizeof(float
849 } else {
850 invertByteOrder((char *)&scale, ftmp[0])
851 ofile.write((char *)ftmp, 4);
852 }
853 if(DEBUG || kVerbose > 0) {
854 G4cout << "ROI image min., max., scale :
855 << minmax[0] << ", "
856 << minmax[1] << ", "
857 << scale << G4endl;
858 }
859
860 // ROI image
861 int rsize = size[0]*size[1];
862 for(int i = 0; i < size[2]; i++) {
863 short * rimage = kRoi[0].getImage(i);
864 if(kLittleEndianOutput) {
865 ofile.write((char *)rimage, rsize*sizeof(sho
866 } else {
867 for(int j = 0; j < rsize; j++) {
868 invertByteOrder((char *)&rimage[j], stmp[0
869 ofile.write((char *)stmp, 2);
870 }
871 }
872
873 }
874
875 // ROI relative location
876 kRoi[0].getCenterPosition(fCenter);
877 for(int i = 0; i < 3; i++) iCenter[i] = (i
878 if(kLittleEndianOutput) {
879 ofile.write((char *)iCenter, 3*sizeof(in
880 } else {
881 for(int j = 0; j < 3; j++)
882 invertByteOrder((char *)&iCenter[j], itmp[j]
883 ofile.write((char *)itmp, 12);
884 }
885 if(DEBUG || kVerbose > 0) {
886 G4cout << "ROI image relative location :
887 << iCenter[0] << ", "
888 << iCenter[1] << ", "
889 << iCenter[2] << ")" << G4endl;
890 }
891 }
892
893 //----- track information -----//
894 // number of track
895 if(kPointerToTrackData > 0) {
896
897 int ntrk = (int)kTracks.size();
898 if(kLittleEndianOutput) {
899 ofile.write((char *)&ntrk, sizeof(int));
900 } else {
901 invertByteOrder((char *)&ntrk, itmp[0]);
902 ofile.write((char *)itmp, 4);
903 }
904 if(DEBUG || kVerbose > 0) {
905 G4cout << "# of tracks : "
906 << ntrk << G4endl;
907 }
908
909 for(int nt = 0; nt < ntrk; nt++) {
910
911 // # of steps in a track
912 int nsteps = kTracks[nt].getNumberOfStep
913 if(kLittleEndianOutput) {
914 ofile.write((char *)&nsteps, sizeof(int));
915 } else {
916 invertByteOrder((char *)&nsteps, itmp[0]);
917 ofile.write((char *)itmp, 4);
918 }
919 if(DEBUG || kVerbose > 0) {
920 G4cout << "# of steps : " << nsteps << G4end
921 }
922
923 // track color
924 unsigned char tcolor[3];
925 kTracks[nt].getColor(tcolor);
926 ofile.write((char *)tcolor, 3);
927
928 // steps
929 float stepPoints[6];
930 for(int isteps = 0; isteps < nsteps; ist
931 kTracks[nt].getStep(stepPoints[0], stepPoint
932 stepPoints[3], stepPoints[4], stepPo
933 isteps);
934
935 if(kLittleEndianOutput) {
936 ofile.write((char *)stepPoints, sizeof(flo
937 } else {
938 for(int j = 0; j < 6; j++)
939 invertByteOrder((char *)&stepPoints[j],
940 ofile.write((char *)ftmp, 24);
941 }
942 }
943 }
944 }
945
946 //----- detector information -----//
947 // number of detectors
948 if(kPointerToDetectorData > 0) {
949 int ndet = (int)kDetectors.size();
950 if(kLittleEndianOutput) {
951 ofile.write((char *)&ndet, sizeof(int));
952 } else {
953 invertByteOrder((char *)&ndet, itmp[0]);
954 ofile.write((char *)itmp, 4);
955 }
956 if(DEBUG || kVerbose > 0) {
957 G4cout << "# of detectors : "
958 << ndet << G4endl;
959 }
960
961 for(int nd = 0; nd < ndet; nd++) {
962
963 // # of edges of a detector
964 int nedges = kDetectors[nd].getNumberOfE
965 if(kLittleEndianOutput) {
966 ofile.write((char *)&nedges, sizeof(int));
967 } else {
968 invertByteOrder((char *)&nedges, itmp[0]);
969 ofile.write((char *)itmp, 4);
970 }
971 if(DEBUG || kVerbose > 0) {
972 G4cout << "# of edges in a detector : " << n
973 }
974
975 // edges
976 float edgePoints[6];
977 for(int ne = 0; ne < nedges; ne++) {
978 kDetectors[nd].getEdge(edgePoints[0], edgePo
979 edgePoints[3], edgePoints[4], edg
980 ne);
981
982 if(kLittleEndianOutput) {
983 ofile.write((char *)edgePoints, sizeof(flo
984 } else {
985 for(int j = 0; j < 6; j++)
986 invertByteOrder((char *)&edgePoints[j],
987 ofile.write((char *)ftmp, 24);
988 }
989
990 if(DEBUG || kVerbose > 0) {
991 if(ne < 1) {
992 G4cout << " edge : (" << edgePoints[0] <
993 << edgePoints[1] << ", "
994 << edgePoints[2] << ") - ("
995 << edgePoints[3] << ", "
996 << edgePoints[4] << ", "
997 << edgePoints[5] << ")" << G4endl;
998 }
999 }
1000 }
1001
1002 // detector color
1003 unsigned char dcolor[3];
1004 kDetectors[nd].getColor(dcolor);
1005 ofile.write((char *)dcolor, 3);
1006 if(DEBUG || kVerbose > 0) {
1007 G4cout << " rgb : (" << (int)dcolor[0] << "
1008 << (int)dcolor[1] << ", "
1009 << (int)dcolor[2] << ")" << G4endl;
1010 }
1011
1012 // detector name
1013 std::string dname = kDetectors[nd].getN
1014 dname.resize(80);
1015 ofile.write((char *)dname.c_str(), 80);
1016 if(DEBUG || kVerbose > 0) {
1017 G4cout << " detector name : " << dname << G
1018
1019 }
1020 }
1021 }
1022
1023 // file end mark
1024 ofile.write("END", 3);
1025
1026 ofile.close();
1027 if(DEBUG || kVerbose > 0)
1028 G4cout << ">>>> closed gdd file: " << kFi
1029
1030 return true;
1031 }
1032 bool G4GMocrenIO::storeData3() {
1033
1034 if(kVerbose > 0) G4cout << ">>>>>>> store
1035 if(kVerbose > 0) G4cout << " " << k
1036
1037 bool DEBUG = false;//
1038
1039 // output file open
1040 std::ofstream ofile(kFileName.c_str(),
1041 std::ios_base::out|std::ios_base::b
1042
1043 // file identifier
1044 ofile.write("gMocren ", 8);
1045
1046 // file version
1047 unsigned char ver = 0x03;
1048 ofile.write((char *)&ver, 1);
1049
1050 // endian
1051 ofile.write((char *)&kLittleEndianOutput, s
1052
1053 // comment length (fixed size)
1054 int commentLength = 1024;
1055 ofile.write((char *)&commentLength, 4);
1056
1057 // comment
1058 char cmt[1025];
1059 std::strncpy(cmt, kComment.c_str(), 1024);
1060 ofile.write((char *)cmt, 1024);
1061 if(DEBUG || kVerbose > 0) {
1062 G4cout << "Data comment : "
1063 << kComment << G4endl;
1064 }
1065
1066 // voxel spacings for all images
1067 ofile.write((char *)kVoxelSpacing, 12);
1068 if(DEBUG || kVerbose > 0) {
1069 G4cout << "Voxel spacing : ("
1070 << kVoxelSpacing[0] << ", "
1071 << kVoxelSpacing[1] << ", "
1072 << kVoxelSpacing[2]
1073 << ") mm " << G4endl;
1074 }
1075
1076 calcPointers3();
1077
1078 // offset from file starting point to the m
1079 ofile.write((char *)&kPointerToModalityData
1080
1081 // # of dose distributions
1082 //int nDoseDist = (int)pointerToDoseDistDat
1083 int nDoseDist = getNumDoseDist();
1084 ofile.write((char *)&nDoseDist, 4);
1085
1086 // offset from file starting point to the d
1087 for(int i = 0; i < nDoseDist; i++) {
1088 ofile.write((char *)&kPointerToDoseDistDa
1089 }
1090
1091 // offset from file starting point to the R
1092 ofile.write((char *)&kPointerToROIData, 4);
1093
1094 // offset from file starting point to the t
1095 ofile.write((char *)&kPointerToTrackData, 4
1096 if(DEBUG || kVerbose > 0) {
1097 G4cout << "Each pointer to data : "
1098 << kPointerToModalityData << ", ";
1099 for(int i = 0; i < nDoseDist; i++) {
1100 G4cout << kPointerToDoseDistData[i] <<
1101 }
1102 G4cout << kPointerToROIData << ", "
1103 << kPointerToTrackData << G4endl;
1104 }
1105
1106 //----- modality image -----//
1107
1108 int size[3];
1109 float scale;
1110 short minmax[2];
1111 float fCenter[3];
1112 int iCenter[3];
1113 // modality image size
1114 kModality.getSize(size);
1115 ofile.write((char *)size, 3*sizeof(int));
1116 if(DEBUG || kVerbose > 0) {
1117 G4cout << "Modality image size : ("
1118 << size[0] << ", "
1119 << size[1] << ", "
1120 << size[2] << ")"
1121 << G4endl;
1122 }
1123
1124 // modality image max. & min.
1125 kModality.getMinMax(minmax);
1126 ofile.write((char *)minmax, 4);
1127
1128 // modality image unit
1129 char munit[13] = "g/cm3 ";
1130 ofile.write((char *)munit, 12);
1131
1132 // modality image scale
1133 scale = (float)kModality.getScale();
1134 ofile.write((char *)&scale, 4);
1135 if(DEBUG || kVerbose > 0) {
1136 G4cout << "Modality image min., max., sca
1137 << minmax[0] << ", "
1138 << minmax[1] << ", "
1139 << scale << G4endl;
1140 }
1141
1142 // modality image
1143 int psize = size[0]*size[1];
1144 if(DEBUG || kVerbose > 0) G4cout << "Modali
1145 for(int i = 0; i < size[2]; i++) {
1146 short * image = kModality.getImage(i);
1147 ofile.write((char *)image, psize*sizeof(s
1148
1149 if(DEBUG || kVerbose > 0) G4cout << "[" <
1150 }
1151 if(DEBUG || kVerbose > 0) G4cout << G4endl;
1152
1153 // modality desity map for CT value
1154 size_t msize = minmax[1] - minmax[0]+1;
1155 float * pdmap = new float[msize];
1156 for(int i = 0; i < (int)msize; i++) pdmap[i
1157 ofile.write((char *)pdmap, msize*sizeof(flo
1158 if(DEBUG || kVerbose > 0) {
1159 G4cout << "density map : " << std::ends;
1160 for(int i = 0; i < (int)msize; i+=50)
1161 G4cout <<kModalityImageDensityMap[i] <<
1162 G4cout << G4endl;
1163 }
1164 delete [] pdmap;
1165
1166
1167 //----- dose distribution image -----//
1168
1169 if(!isDoseEmpty()) {
1170
1171 calcDoseDistScale();
1172
1173 for(int ndose = 0; ndose < nDoseDist; ndo
1174 // dose distrbution image size
1175 kDose[ndose].getSize(size);
1176 ofile.write((char *)size, 3*sizeof(int)
1177 if(DEBUG || kVerbose > 0) {
1178 G4cout << "Dose dist. [" << ndose << "] ima
1179 << size[0] << ", "
1180 << size[1] << ", "
1181 << size[2] << ")"
1182 << G4endl;
1183 }
1184
1185 // dose distribution max. & min.
1186 getShortDoseDistMinMax(minmax, ndose);
1187 ofile.write((char *)minmax, 2*2); // si
1188
1189 // dose distribution unit
1190 ofile.write((char *)kDoseUnit.c_str(),
1191 if(DEBUG || kVerbose > 0) {
1192 G4cout << "Dose dist. unit : " << kDoseUnit
1193 }
1194
1195 // dose distribution scaling
1196 double dscale;
1197 dscale = getDoseDistScale(ndose);
1198 scale = float(dscale);
1199 ofile.write((char *)&scale, 4);
1200 if(DEBUG || kVerbose > 0) {
1201 G4cout << "Dose dist. [" << ndose
1202 << "] image min., max., scale : "
1203 << minmax[0] << ", "
1204 << minmax[1] << ", "
1205 << scale << G4endl;
1206 }
1207
1208 // dose distribution image
1209 int dsize = size[0]*size[1];
1210 short * dimage = new short[dsize];
1211 for(int z = 0; z < size[2]; z++) {
1212 getShortDoseDist(dimage, z, ndose);
1213 ofile.write((char *)dimage, dsize*2); //siz
1214
1215 if(DEBUG || kVerbose > 0) {
1216 for(int j = 0; j < dsize; j++) {
1217 if(dimage[j] < 0)
1218 G4cout << "[" << j << "," << z << "]"
1219 << dimage[j] << ", ";
1220 }
1221 }
1222 }
1223 if(DEBUG || kVerbose > 0) G4cout << G4e
1224 delete [] dimage;
1225
1226 // relative location of the dose distri
1227 // the modality image
1228 getDoseDistCenterPosition(fCenter, ndos
1229 for(int i = 0; i < 3; i++) iCenter[i] =
1230 ofile.write((char *)iCenter, 3*4); // 3
1231 if(DEBUG || kVerbose > 0) {
1232 G4cout << "Dose dist. [" << ndose
1233 << "]image relative location : ("
1234 << iCenter[0] << ", "
1235 << iCenter[1] << ", "
1236 << iCenter[2] << ")" << G4endl;
1237 }
1238 }
1239 }
1240
1241 //----- ROI image -----//
1242 if(!isROIEmpty()) {
1243 // ROI image size
1244 kRoi[0].getSize(size);
1245 ofile.write((char *)size, 3*sizeof(int));
1246 if(DEBUG || kVerbose > 0) {
1247 G4cout << "ROI image size : ("
1248 << size[0] << ", "
1249 << size[1] << ", "
1250 << size[2] << ")"
1251 << G4endl;
1252 }
1253
1254 // ROI max. & min.
1255 kRoi[0].getMinMax(minmax);
1256 ofile.write((char *)minmax, sizeof(short)
1257
1258 // ROI distribution scaling
1259 scale = (float)kRoi[0].getScale();
1260 ofile.write((char *)&scale, sizeof(float)
1261 if(DEBUG || kVerbose > 0) {
1262 G4cout << "ROI image min., max., scale
1263 << minmax[0] << ", "
1264 << minmax[1] << ", "
1265 << scale << G4endl;
1266 }
1267
1268 // ROI image
1269 int rsize = size[0]*size[1];
1270 for(int i = 0; i < size[2]; i++) {
1271 short * rimage = kRoi[0].getImage(i);
1272 ofile.write((char *)rimage, rsize*sizeo
1273
1274 }
1275
1276 // ROI relative location
1277 kRoi[0].getCenterPosition(fCenter);
1278 for(int i = 0; i < 3; i++) iCenter[i] = (
1279 ofile.write((char *)iCenter, 3*sizeof(int
1280 if(DEBUG || kVerbose > 0) {
1281 G4cout << "ROI image relative location
1282 << iCenter[0] << ", "
1283 << iCenter[1] << ", "
1284 << iCenter[2] << ")" << G4endl;
1285 }
1286 }
1287
1288 //----- track information -----//
1289 // number of track
1290 int ntrk = (int)kSteps.size();
1291 ofile.write((char *)&ntrk, sizeof(int));
1292 if(DEBUG || kVerbose > 0) {
1293 G4cout << "# of tracks : "
1294 << ntrk << G4endl;
1295 }
1296 // track position
1297 for(int i = 0; i < ntrk; i++) {
1298 float * tp = kSteps[i];
1299 ofile.write((char *)tp, sizeof(float)*6);
1300 }
1301 // track color
1302 int ntcolor = int(kStepColors.size());
1303 if(ntrk != ntcolor)
1304 if (G4VisManager::GetVerbosity() >= G4Vis
1305 G4cout << "# of track color information
1306 << G4endl;
1307 unsigned char white[3] = {255,255,255}; //
1308 for(int i = 0; i < ntrk; i++) {
1309 if(i < ntcolor) {
1310 unsigned char * tcolor = kStepColors[i]
1311 ofile.write((char *)tcolor, 3);
1312 } else {
1313 ofile.write((char *)white, 3);
1314 }
1315 }
1316
1317 // file end mark
1318 ofile.write("END", 3);
1319
1320 ofile.close();
1321
1322 return true;
1323 }
1324 //
1325 bool G4GMocrenIO::storeData4(char * _filename
1326 kFileName = _filename;
1327 return storeData4();
1328 }
1329
1330 // version 2
1331 bool G4GMocrenIO::storeData2() {
1332
1333 if(kVerbose > 0) G4cout << ">>>>>>> store
1334 if(kVerbose > 0) G4cout << " " << k
1335
1336 bool DEBUG = false;//
1337
1338 // output file open
1339 std::ofstream ofile(kFileName.c_str(),
1340 std::ios_base::out|std::ios_base::b
1341
1342 // file identifier
1343 ofile.write("GRAPE ", 8);
1344
1345 // file version
1346 unsigned char ver = 0x02;
1347 ofile.write((char *)&ver, 1);
1348 // file id for old file format support
1349 ofile.write(kId.c_str(), IDLENGTH);
1350 // file version for old file format support
1351 ofile.write(kVersion.c_str(), VERLENGTH);
1352 // endian
1353 ofile.write((char *)&kLittleEndianOutput, s
1354
1355 /*
1356 // event number
1357 ofile.write((char *)&numberOfEvents, sizeof
1358 float imageSpacing[3];
1359 imageSpacing[0] = modalityImageVoxelSpacing
1360 imageSpacing[1] = modalityImageVoxelSpacing
1361 imageSpacing[2] = modalityImageVoxelSpacing
1362 ofile.write((char *)imageSpacing, 12);
1363 */
1364
1365
1366 // voxel spacings for all images
1367 ofile.write((char *)kVoxelSpacing, 12);
1368 if(DEBUG || kVerbose > 0) {
1369 G4cout << "Voxel spacing : ("
1370 << kVoxelSpacing[0] << ", "
1371 << kVoxelSpacing[1] << ", "
1372 << kVoxelSpacing[2]
1373 << ") mm " << G4endl;
1374 }
1375
1376 calcPointers2();
1377 // offset from file starting point to the m
1378 ofile.write((char *)&kPointerToModalityData
1379
1380 // offset from file starting point to the d
1381 ofile.write((char *)&kPointerToDoseDistData
1382
1383 // offset from file starting point to the R
1384 ofile.write((char *)&kPointerToROIData, 4);
1385
1386 // offset from file starting point to the t
1387 ofile.write((char *)&kPointerToTrackData, 4
1388 if(DEBUG || kVerbose > 0) {
1389 G4cout << "Each pointer to data : "
1390 << kPointerToModalityData << ", "
1391 << kPointerToDoseDistData[0] << ", "
1392 << kPointerToROIData << ", "
1393 << kPointerToTrackData << G4endl;
1394 }
1395
1396 //----- modality image -----//
1397
1398 int size[3];
1399 float scale;
1400 short minmax[2];
1401 float fCenter[3];
1402 int iCenter[3];
1403 // modality image size
1404 kModality.getSize(size);
1405 ofile.write((char *)size, 3*sizeof(int));
1406 if(DEBUG || kVerbose > 0) {
1407 G4cout << "Modality image size : ("
1408 << size[0] << ", "
1409 << size[1] << ", "
1410 << size[2] << ")"
1411 << G4endl;
1412 }
1413
1414 // modality image max. & min.
1415 kModality.getMinMax(minmax);
1416 ofile.write((char *)minmax, 4);
1417
1418 // modality image unit
1419 //char munit[13] = "g/cm3 ";
1420 //ofile.write((char *)&munit, 12);
1421
1422 // modality image scale
1423 scale = (float)kModality.getScale();
1424 ofile.write((char *)&scale, 4);
1425 if(DEBUG || kVerbose > 0) {
1426 G4cout << "Modality image min., max., sca
1427 << minmax[0] << ", "
1428 << minmax[1] << ", "
1429 << scale << G4endl;
1430 }
1431
1432 // modality image
1433 int psize = size[0]*size[1];
1434 if(DEBUG || kVerbose > 0) G4cout << "Modali
1435 for(int i = 0; i < size[2]; i++) {
1436 short * image =kModality.getImage(i);
1437 ofile.write((char *)image, psize*sizeof(s
1438
1439 if(DEBUG || kVerbose > 0) G4cout << "[" <
1440 }
1441 if(DEBUG || kVerbose > 0) G4cout << G4endl;
1442
1443 // modality desity map for CT value
1444 size_t msize = minmax[1] - minmax[0]+1;
1445 float * pdmap = new float[msize];
1446 for(int i = 0; i < (int)msize; i++) pdmap[i
1447 ofile.write((char *)pdmap, msize*sizeof(flo
1448 if(DEBUG || kVerbose > 0) {
1449 G4cout << "density map : " << std::ends;
1450 for(int i = 0; i < (int)msize; i+=50)
1451 G4cout <<kModalityImageDensityMap[i] <<
1452 G4cout << G4endl;
1453 }
1454 delete [] pdmap;
1455
1456
1457 //----- dose distribution image -----//
1458
1459 if(!isDoseEmpty()) {
1460 calcDoseDistScale();
1461
1462 // dose distrbution image size
1463 kDose[0].getSize(size);
1464 ofile.write((char *)size, 3*sizeof(int));
1465 if(DEBUG || kVerbose > 0) {
1466 G4cout << "Dose dist. image size : ("
1467 << size[0] << ", "
1468 << size[1] << ", "
1469 << size[2] << ")"
1470 << G4endl;
1471 }
1472
1473 // dose distribution max. & min.
1474 getShortDoseDistMinMax(minmax);
1475 ofile.write((char *)minmax, sizeof(short)
1476
1477 // dose distribution scaling
1478 scale = (float)kDose[0].getScale();
1479 ofile.write((char *)&scale, sizeof(float)
1480 if(DEBUG || kVerbose > 0) {
1481 G4cout << "Dose dist. image min., max.,
1482 << minmax[0] << ", "
1483 << minmax[1] << ", "
1484 << scale << G4endl;
1485 }
1486
1487 // dose distribution image
1488 int dsize = size[0]*size[1];
1489 short * dimage = new short[dsize];
1490 for(int z = 0; z < size[2]; z++) {
1491 getShortDoseDist(dimage, z);
1492 ofile.write((char *)dimage, dsize*sizeo
1493
1494 if(DEBUG || kVerbose > 0) {
1495 for(int j = 0; j < dsize; j++) {
1496 if(dimage[j] < 0)
1497 G4cout << "[" << j << "," << z << "]"
1498 << dimage[j] << ", ";
1499 }
1500 }
1501 }
1502 if(DEBUG || kVerbose > 0) G4cout << G4end
1503 delete [] dimage;
1504
1505 // relative location of the dose distribu
1506 // the modality image
1507 kDose[0].getCenterPosition(fCenter);
1508 for(int i = 0; i < 3; i++) iCenter[i] = (
1509 ofile.write((char *)iCenter, 3*sizeof(int
1510 if(DEBUG || kVerbose > 0) {
1511 G4cout << "Dose dist. image relative lo
1512 << iCenter[0] << ", "
1513 << iCenter[1] << ", "
1514 << iCenter[2] << ")" << G4endl;
1515 }
1516
1517 }
1518
1519 //----- ROI image -----//
1520 if(!isROIEmpty()) {
1521 // ROI image size
1522 kRoi[0].getSize(size);
1523 ofile.write((char *)size, 3*sizeof(int));
1524 if(DEBUG || kVerbose > 0) {
1525 G4cout << "ROI image size : ("
1526 << size[0] << ", "
1527 << size[1] << ", "
1528 << size[2] << ")"
1529 << G4endl;
1530 }
1531
1532 // ROI max. & min.
1533 kRoi[0].getMinMax(minmax);
1534 ofile.write((char *)minmax, sizeof(short)
1535
1536 // ROI distribution scaling
1537 scale = (float)kRoi[0].getScale();
1538 ofile.write((char *)&scale, sizeof(float)
1539 if(DEBUG || kVerbose > 0) {
1540 G4cout << "ROI image min., max., scale
1541 << minmax[0] << ", "
1542 << minmax[1] << ", "
1543 << scale << G4endl;
1544 }
1545
1546 // ROI image
1547 int rsize = size[0]*size[1];
1548 for(int i = 0; i < size[2]; i++) {
1549 short * rimage = kRoi[0].getImage(i);
1550 ofile.write((char *)rimage, rsize*sizeo
1551
1552 }
1553
1554 // ROI relative location
1555 kRoi[0].getCenterPosition(fCenter);
1556 for(int i = 0; i < 3; i++) iCenter[i] = (
1557 ofile.write((char *)iCenter, 3*sizeof(int
1558 if(DEBUG || kVerbose > 0) {
1559 G4cout << "ROI image relative location
1560 << iCenter[0] << ", "
1561 << iCenter[1] << ", "
1562 << iCenter[2] << ")" << G4endl;
1563 }
1564 }
1565
1566
1567 //----- track information -----//
1568 // track
1569 int ntrk = (int)kSteps.size();
1570 ofile.write((char *)&ntrk, sizeof(int));
1571 if(DEBUG || kVerbose > 0) {
1572 G4cout << "# of tracks : "
1573 << ntrk << G4endl;
1574 }
1575 for(int i = 0; i < ntrk; i++) {
1576 float * tp = kSteps[i];
1577 ofile.write((char *)tp, sizeof(float)*6);
1578 }
1579
1580
1581 // file end mark
1582 ofile.write("END", 3);
1583
1584 ofile.close();
1585
1586 return true;
1587 }
1588 //
1589 bool G4GMocrenIO::storeData2(char * _filename
1590 kFileName = _filename;
1591 return storeData();
1592 }
1593
1594 bool G4GMocrenIO::retrieveData() {
1595
1596 // input file open
1597 std::ifstream ifile(kFileName.c_str(), std:
1598 if(!ifile) {
1599 if (G4VisManager::GetVerbosity() >= G4Vis
1600 G4cout << "Cannot open file: " << kFile
1601 << " in G4GMocrenIO::retrieveData()."
1602 return false;
1603 }
1604
1605 // file identifier
1606 char verid[9];
1607 ifile.read((char *)verid, 8);
1608 // file version
1609 unsigned char ver;
1610 ifile.read((char *)&ver, 1);
1611 ifile.close();
1612
1613 if(std::strncmp(verid, "gMocren", 7) == 0)
1614 if(ver == 0x03) {
1615 G4cout << ">>>>>>> retrieve data (ver.
1616 G4cout << " " << kFileName << G
1617 retrieveData3();
1618 } else if (ver == 0x04) {
1619 G4cout << ">>>>>>> retrieve data (ver.
1620 G4cout << " " << kFileName << G
1621 retrieveData4();
1622 } else {
1623 if (G4VisManager::GetVerbosity() >= G4V
1624 G4cout << "Error -- invalid file version :
1625 << G4endl;
1626 G4cout << " " << kFileName << G4end
1627 }
1628 G4Exception("G4GMocrenIO::retrieveDadta
1629 "gMocren2001", FatalExcepti
1630 "Error.");
1631
1632 }
1633 } else if(std::strncmp(verid, "GRAPE", 5) =
1634 G4cout << ">>>>>>> retrieve data (ver.2)
1635 G4cout << " " << kFileName << G4e
1636 retrieveData2();
1637 } else {
1638 if (G4VisManager::GetVerbosity() >= G4Vis
1639 G4cout << kFileName << " was not gdd fi
1640 return false;
1641 }
1642
1643 return true;
1644 }
1645
1646 bool G4GMocrenIO::retrieveData(char * _filena
1647 kFileName = _filename;
1648 return retrieveData();
1649 }
1650
1651 //
1652 bool G4GMocrenIO::retrieveData4() {
1653
1654 bool DEBUG = false;//
1655
1656 // input file open
1657 std::ifstream ifile(kFileName.c_str(), std:
1658 if(!ifile) {
1659 if (G4VisManager::GetVerbosity() >= G4Vis
1660 G4cout << "Cannot open file: " << kFile
1661 << " in G4GMocrenIO::retrieveData3()." <<
1662 return false;
1663 }
1664
1665 // data buffer
1666 char ctmp[24];
1667
1668 // file identifier
1669 char verid[9];
1670 ifile.read((char *)verid, 8);
1671
1672 // file version
1673 unsigned char ver;
1674 ifile.read((char *)&ver, 1);
1675 std::stringstream ss;
1676 ss << (int)ver;
1677 kVersion = ss.str();
1678 if(DEBUG || kVerbose > 0) G4cout << "File v
1679
1680 // endian
1681 ifile.read((char *)&kLittleEndianInput, siz
1682 if(DEBUG || kVerbose > 0) {
1683 G4cout << "Endian : ";
1684 if(kLittleEndianInput == 1)
1685 G4cout << " little" << G4endl;
1686 else {
1687 G4cout << " big" << G4endl;
1688 }
1689 }
1690
1691 // comment length (fixed size)
1692 int clength;
1693 ifile.read((char *)ctmp, 4);
1694 convertEndian(ctmp, clength);
1695 // comment
1696 char cmt[1025];
1697 ifile.read((char *)cmt, clength);
1698 std::string scmt = cmt;
1699 scmt += '\0';
1700 setComment(scmt);
1701 if(DEBUG || kVerbose > 0) {
1702 G4cout << "Data comment : "
1703 << kComment << G4endl;
1704 }
1705
1706 // voxel spacings for all images
1707 ifile.read((char *)ctmp, 12);
1708 convertEndian(ctmp, kVoxelSpacing[0]);
1709 convertEndian(ctmp+4, kVoxelSpacing[1]);
1710 convertEndian(ctmp+8, kVoxelSpacing[2]);
1711 if(DEBUG || kVerbose > 0) {
1712 G4cout << "Voxel spacing : ("
1713 << kVoxelSpacing[0] << ", "
1714 << kVoxelSpacing[1] << ", "
1715 << kVoxelSpacing[2]
1716 << ") mm " << G4endl;
1717 }
1718
1719
1720 // offset from file starting point to the m
1721 ifile.read((char *)ctmp, 4);
1722 convertEndian(ctmp, kPointerToModalityData)
1723
1724 // # of dose distributions
1725 ifile.read((char *)ctmp, 4);
1726 int nDoseDist;
1727 convertEndian(ctmp, nDoseDist);
1728
1729 // offset from file starting point to the d
1730 for(int i = 0; i < nDoseDist; i++) {
1731 ifile.read((char *)ctmp, 4);
1732 unsigned int dptr;
1733 convertEndian(ctmp, dptr);
1734 addPointerToDoseDistData(dptr);
1735 }
1736
1737 // offset from file starting point to the R
1738 ifile.read((char *)ctmp, 4);
1739 convertEndian(ctmp, kPointerToROIData);
1740
1741 // offset from file starting point to the t
1742 ifile.read((char *)ctmp, 4);
1743 convertEndian(ctmp, kPointerToTrackData);
1744
1745 // offset from file starting point to the d
1746 ifile.read((char *)ctmp, 4);
1747 convertEndian(ctmp, kPointerToDetectorData)
1748
1749 if(DEBUG || kVerbose > 0) {
1750 G4cout << "Each pointer to data : "
1751 << kPointerToModalityData << ", ";
1752 for(int i = 0; i < nDoseDist; i++)
1753 G4cout << kPointerToDoseDistData[i] <<
1754 G4cout << kPointerToROIData << ", "
1755 << kPointerToTrackData << ", "
1756 << kPointerToDetectorData
1757 << G4endl;
1758 }
1759
1760
1761
1762 if(kPointerToModalityData == 0 && kPointerT
1763 kPointerToROIData == 0 && kPointerToTrac
1764 if(DEBUG || kVerbose > 0) {
1765 G4cout << "No data." << G4endl;
1766 }
1767 return false;
1768 }
1769
1770 // event number
1771 /* ver 1
1772 ifile.read(ctmp, sizeof(int));
1773 convertEndian(ctmp, numberOfEvents);
1774 */
1775
1776 int size[3];
1777 float scale;
1778 double dscale;
1779 short minmax[2];
1780 float fCenter[3];
1781 int iCenter[3];
1782
1783 //----- Modality image -----//
1784 // modality image size
1785 ifile.read(ctmp, 3*sizeof(int));
1786 convertEndian(ctmp, size[0]);
1787 convertEndian(ctmp+sizeof(int), size[1]);
1788 convertEndian(ctmp+2*sizeof(int), size[2]);
1789 if(DEBUG || kVerbose > 0) {
1790 G4cout << "Modality image size : ("
1791 << size[0] << ", "
1792 << size[1] << ", "
1793 << size[2] << ")"
1794 << G4endl;
1795 }
1796 kModality.setSize(size);
1797
1798 // modality image voxel spacing
1799 /*
1800 ifile.read(ctmp, 3*sizeof(float));
1801 convertEndian(ctmp, modalityImageVoxelSpa
1802 convertEndian(ctmp+sizeof(float), modalit
1803 convertEndian(ctmp+2*sizeof(float), modal
1804 */
1805
1806 if(kPointerToModalityData != 0) {
1807
1808 // modality density max. & min.
1809 ifile.read((char *)ctmp, 4);
1810 convertEndian(ctmp, minmax[0]);
1811 convertEndian(ctmp+2, minmax[1]);
1812 kModality.setMinMax(minmax);
1813
1814 // modality image unit
1815 char munit[13];
1816 munit[12] = '\0';
1817 ifile.read((char *)munit, 12);
1818 std::string smunit = munit;
1819 setModalityImageUnit(smunit);
1820
1821 // modality density scale
1822 ifile.read((char *)ctmp, 4);
1823 convertEndian(ctmp, scale);
1824 kModality.setScale(dscale = scale);
1825 if(DEBUG || kVerbose > 0) {
1826 G4cout << "Modality image min., max., s
1827 << minmax[0] << ", "
1828 << minmax[1] << ", "
1829 << scale << G4endl;
1830 }
1831
1832 // modality density
1833 int psize = size[0]*size[1];
1834 if(DEBUG || kVerbose > 0) G4cout << "Moda
1835 char * cimage = new char[psize*sizeof(sho
1836 for(int i = 0; i < size[2]; i++) {
1837 ifile.read((char *)cimage, psize*sizeof
1838 short * mimage = new short[psize];
1839 for(int j = 0; j < psize; j++) {
1840 convertEndian(cimage+j*sizeof(short), mimag
1841 }
1842 kModality.addImage(mimage);
1843
1844 if(DEBUG || kVerbose > 0) G4cout << "["
1845 }
1846 if(DEBUG || kVerbose > 0) G4cout << G4end
1847 delete [] cimage;
1848
1849 // modality desity map for CT value
1850 size_t msize = minmax[1]-minmax[0]+1;
1851 if(DEBUG || kVerbose > 0) G4cout << "msiz
1852 char * pdmap = new char[msize*sizeof(floa
1853 ifile.read((char *)pdmap, msize*sizeof(fl
1854 float ftmp;
1855 for(int i = 0; i < (int)msize; i++) {
1856 convertEndian(pdmap+i*sizeof(float), ft
1857 kModalityImageDensityMap.push_back(ftmp
1858 }
1859 delete [] pdmap;
1860
1861 if(DEBUG || kVerbose > 0) {
1862 G4cout << "density map : " << std::ends
1863 for(int i = 0; i < 10; i++)
1864 G4cout <<kModalityImageDensityMap[i] << ",
1865 G4cout << G4endl;
1866 for(int i = 0; i < 10; i++) G4cout << "
1867 G4cout << G4endl;
1868 for(size_t i =kModalityImageDensityMap.
1869 G4cout <<kModalityImageDensityMap[i] << ",
1870 G4cout << G4endl;
1871 }
1872
1873 }
1874
1875
1876 //----- dose distribution image -----//
1877 for(int ndose = 0; ndose < nDoseDist; ndose
1878
1879 newDoseDist();
1880
1881 // dose distrbution image size
1882 ifile.read((char *)ctmp, 3*sizeof(int));
1883 convertEndian(ctmp, size[0]);
1884 convertEndian(ctmp+sizeof(int), size[1]);
1885 convertEndian(ctmp+2*sizeof(int), size[2]
1886 if(DEBUG || kVerbose > 0) {
1887 G4cout << "Dose dist. image size : ("
1888 << size[0] << ", "
1889 << size[1] << ", "
1890 << size[2] << ")"
1891 << G4endl;
1892 }
1893 kDose[ndose].setSize(size);
1894
1895 // dose distribution max. & min.
1896 ifile.read((char *)ctmp, sizeof(short)*2)
1897 convertEndian(ctmp, minmax[0]);
1898 convertEndian(ctmp+2, minmax[1]);
1899
1900 // dose distribution unit
1901 char dunit[13];
1902 dunit[12] = '\0';
1903 ifile.read((char *)dunit, 12);
1904 std::string sdunit = dunit;
1905 setDoseDistUnit(sdunit, ndose);
1906 if(DEBUG || kVerbose > 0) {
1907 G4cout << "Dose dist. unit : " << kDose
1908 }
1909
1910 // dose distribution scaling
1911 ifile.read((char *)ctmp, 4); // sizeof(fl
1912 convertEndian(ctmp, scale);
1913 kDose[ndose].setScale(dscale = scale);
1914
1915 double dminmax[2];
1916 for(int i = 0; i < 2; i++) dminmax[i] = m
1917 kDose[ndose].setMinMax(dminmax);
1918
1919 if(DEBUG || kVerbose > 0) {
1920 G4cout << "Dose dist. image min., max.,
1921 << dminmax[0] << ", "
1922 << dminmax[1] << ", "
1923 << scale << G4endl;
1924 }
1925
1926 // dose distribution image
1927 int dsize = size[0]*size[1];
1928 if(DEBUG || kVerbose > 0) G4cout << "Dose
1929 char * di = new char[dsize*sizeof(short)]
1930 short * shimage = new short[dsize];
1931 for(int z = 0; z < size[2]; z++) {
1932 ifile.read((char *)di, dsize*sizeof(sho
1933 double * dimage = new double[dsize];
1934 for(int xy = 0; xy < dsize; xy++) {
1935 convertEndian(di+xy*sizeof(short), shimage[
1936 dimage[xy] = shimage[xy]*dscale;
1937 }
1938 kDose[ndose].addImage(dimage);
1939
1940 if(DEBUG || kVerbose > 0) G4cout << "["
1941
1942 if(DEBUG || kVerbose > 0) {
1943 for(int j = 0; j < dsize; j++) {
1944 if(dimage[j] < 0)
1945 G4cout << "[" << j << "," << z << "]"
1946 << dimage[j] << ", ";
1947 }
1948 }
1949 }
1950 delete [] shimage;
1951 delete [] di;
1952 if(DEBUG || kVerbose > 0) G4cout << G4end
1953
1954 ifile.read((char *)ctmp, 3*4); // 3*sizeo
1955 convertEndian(ctmp, iCenter[0]);
1956 convertEndian(ctmp+4, iCenter[1]);
1957 convertEndian(ctmp+8, iCenter[2]);
1958 for(int i = 0; i < 3; i++) fCenter[i] = (
1959 kDose[ndose].setCenterPosition(fCenter);
1960
1961 if(DEBUG || kVerbose > 0) {
1962 G4cout << "Dose dist. image relative lo
1963 << fCenter[0] << ", "
1964 << fCenter[1] << ", "
1965 << fCenter[2] << ")" << G4endl;
1966 }
1967
1968
1969 // dose distribution name
1970 char cname[81];
1971 ifile.read((char *)cname, 80);
1972 std::string dosename = cname;
1973 setDoseDistName(dosename, ndose);
1974 if(DEBUG || kVerbose > 0) {
1975 G4cout << "Dose dist. name : " << dosen
1976 }
1977
1978 }
1979
1980 //----- ROI image -----//
1981 if(kPointerToROIData != 0) {
1982
1983 newROI();
1984
1985 // ROI image size
1986 ifile.read((char *)ctmp, 3*sizeof(int));
1987 convertEndian(ctmp, size[0]);
1988 convertEndian(ctmp+sizeof(int), size[1]);
1989 convertEndian(ctmp+2*sizeof(int), size[2]
1990 kRoi[0].setSize(size);
1991 if(DEBUG || kVerbose > 0) {
1992 G4cout << "ROI image size : ("
1993 << size[0] << ", "
1994 << size[1] << ", "
1995 << size[2] << ")"
1996 << G4endl;
1997 }
1998
1999 // ROI max. & min.
2000 ifile.read((char *)ctmp, sizeof(short)*2)
2001 convertEndian(ctmp, minmax[0]);
2002 convertEndian(ctmp+sizeof(short), minmax[
2003 kRoi[0].setMinMax(minmax);
2004
2005 // ROI distribution scaling
2006 ifile.read((char *)ctmp, sizeof(float));
2007 convertEndian(ctmp, scale);
2008 kRoi[0].setScale(dscale = scale);
2009 if(DEBUG || kVerbose > 0) {
2010 G4cout << "ROI image min., max., scale
2011 << minmax[0] << ", "
2012 << minmax[1] << ", "
2013 << scale << G4endl;
2014 }
2015
2016 // ROI image
2017 int rsize = size[0]*size[1];
2018 char * ri = new char[rsize*sizeof(short)]
2019 for(int i = 0; i < size[2]; i++) {
2020 ifile.read((char *)ri, rsize*sizeof(sho
2021 short * rimage = new short[rsize];
2022 for(int j = 0; j < rsize; j++) {
2023 convertEndian(ri+j*sizeof(short), rimage[j]
2024 }
2025 kRoi[0].addImage(rimage);
2026
2027 }
2028 delete [] ri;
2029
2030 // ROI relative location
2031 ifile.read((char *)ctmp, 3*sizeof(int));
2032 convertEndian(ctmp, iCenter[0]);
2033 convertEndian(ctmp+sizeof(int), iCenter[1
2034 convertEndian(ctmp+2*sizeof(int), iCenter
2035 for(int i = 0; i < 3; i++) fCenter[i] = i
2036 kRoi[0].setCenterPosition(fCenter);
2037 if(DEBUG || kVerbose > 0) {
2038 G4cout << "ROI image relative location
2039 << fCenter[0] << ", "
2040 << fCenter[1] << ", "
2041 << fCenter[2] << ")" << G4endl;
2042 }
2043
2044 }
2045
2046 //----- track information -----//
2047 if(kPointerToTrackData != 0) {
2048
2049 // track
2050 ifile.read((char *)ctmp, sizeof(int));
2051 int ntrk;
2052 convertEndian(ctmp, ntrk);
2053 if(DEBUG || kVerbose > 0) {
2054 G4cout << "# of tracks: " << ntrk << G4
2055 }
2056
2057 // track position
2058 unsigned char rgb[3];
2059 for(int i = 0; i < ntrk; i++) {
2060
2061
2062 // # of steps in a track
2063 ifile.read((char *)ctmp, sizeof(int));
2064 int nsteps;
2065 convertEndian(ctmp, nsteps);
2066
2067 // track color
2068 ifile.read((char *)rgb, 3);
2069
2070 std::vector<float *> steps;
2071 // steps
2072 for(int j = 0; j < nsteps; j++) {
2073
2074 float * steppoint = new float[6];
2075 ifile.read((char *)ctmp, sizeof(float)*6);
2076
2077 for(int k = 0; k < 6; k++) {
2078 convertEndian(ctmp+k*sizeof(float), stepp
2079 }
2080
2081 steps.push_back(steppoint);
2082 }
2083
2084 // add a track to the track container
2085 addTrack(steps, rgb);
2086
2087 if(DEBUG || kVerbose > 0) {
2088 if(i < 5) {
2089 G4cout << i << ": " ;
2090 for(int j = 0; j < 3; j++) G4cout << step
2091 int nstp = (int)steps.size();
2092 G4cout << "<-> ";
2093 for(int j = 3; j < 6; j++) G4cout << step
2094 G4cout << " rgb( ";
2095 for(int j = 0; j < 3; j++) G4cout << (int
2096 G4cout << ")" << G4endl;
2097 }
2098 }
2099 }
2100
2101
2102 }
2103
2104
2105 //----- detector information -----//
2106 if(kPointerToDetectorData != 0) {
2107
2108 // number of detectors
2109 ifile.read((char *)ctmp, sizeof(int));
2110 int ndet;
2111 convertEndian(ctmp, ndet);
2112
2113 if(DEBUG || kVerbose > 0) {
2114 G4cout << "# of detectors : "
2115 << ndet << G4endl;
2116 }
2117
2118 for(int nd = 0; nd < ndet; nd++) {
2119
2120 // # of edges of a detector
2121 ifile.read((char *)ctmp, sizeof(int));
2122 int nedges;
2123 convertEndian(ctmp, nedges);
2124 if(DEBUG || kVerbose > 0) {
2125 G4cout << "# of edges in a detector : " <<
2126 }
2127
2128 // edges
2129 std::vector<float *> detector;
2130 char cftmp[24];
2131 for(int ne = 0; ne < nedges; ne++) {
2132
2133 ifile.read((char *)cftmp, sizeof(float)*6);
2134 float * edgePoints = new float[6];
2135 for(int j = 0; j < 6; j++) convertEndian(&c
2136 detector.push_back(edgePoints);
2137
2138 }
2139
2140 if(DEBUG || kVerbose > 0) {
2141 G4cout << " first edge : (" << detector[0][
2142 << detector[0][1] << ", "
2143 << detector[0][2] << ") - ("
2144 << detector[0][3] << ", "
2145 << detector[0][4] << ", "
2146 << detector[0][5] << ")" << G4endl;
2147 }
2148
2149 // detector color
2150 unsigned char dcolor[3];
2151 ifile.read((char *)dcolor, 3);
2152 if(DEBUG || kVerbose > 0) {
2153 G4cout << " detector color : rgb("
2154 << (int)dcolor[0] << ", "
2155 << (int)dcolor[1] << ", "
2156 << (int)dcolor[2] << G4endl;
2157 }
2158
2159
2160 // detector name
2161 char cname[80];
2162 ifile.read((char *)cname, 80);
2163 std::string dname = cname;
2164 if(DEBUG || kVerbose > 0) {
2165 G4cout << " detector name : " << dname << G
2166 }
2167
2168
2169 addDetector(dname, detector, dcolor);
2170
2171 }
2172 }
2173
2174
2175 ifile.close();
2176
2177 return true;
2178 }
2179 bool G4GMocrenIO::retrieveData4(char * _filen
2180 kFileName = _filename;
2181 return retrieveData();
2182 }
2183
2184 //
2185 bool G4GMocrenIO::retrieveData3() {
2186
2187 bool DEBUG = false;//
2188
2189 // input file open
2190 std::ifstream ifile(kFileName.c_str(), std:
2191 if(!ifile) {
2192 if (G4VisManager::GetVerbosity() >= G4Vis
2193 G4cout << "Cannot open file: " << kFile
2194 << " in G4GMocrenIO::retrieveData3()." <<
2195 return false;
2196 }
2197
2198 // data buffer
2199 char ctmp[12];
2200
2201 // file identifier
2202 char verid[9];
2203 ifile.read((char *)verid, 8);
2204
2205 // file version
2206 unsigned char ver;
2207 ifile.read((char *)&ver, 1);
2208 std::stringstream ss;
2209 ss << (int)ver;
2210 kVersion = ss.str();
2211 if(DEBUG || kVerbose > 0) G4cout << "File v
2212
2213 // endian
2214 ifile.read((char *)&kLittleEndianInput, siz
2215 if(DEBUG || kVerbose > 0) {
2216 G4cout << "Endian : ";
2217 if(kLittleEndianInput == 1)
2218 G4cout << " little" << G4endl;
2219 else {
2220 G4cout << " big" << G4endl;
2221 }
2222 }
2223
2224 // comment length (fixed size)
2225 int clength;
2226 ifile.read((char *)ctmp, 4);
2227 convertEndian(ctmp, clength);
2228 // comment
2229 char cmt[1025];
2230 ifile.read((char *)cmt, clength);
2231 std::string scmt = cmt;
2232 setComment(scmt);
2233 if(DEBUG || kVerbose > 0) {
2234 G4cout << "Data comment : "
2235 << kComment << G4endl;
2236 }
2237
2238 // voxel spacings for all images
2239 ifile.read((char *)ctmp, 12);
2240 convertEndian(ctmp, kVoxelSpacing[0]);
2241 convertEndian(ctmp+4, kVoxelSpacing[1]);
2242 convertEndian(ctmp+8, kVoxelSpacing[2]);
2243 if(DEBUG || kVerbose > 0) {
2244 G4cout << "Voxel spacing : ("
2245 << kVoxelSpacing[0] << ", "
2246 << kVoxelSpacing[1] << ", "
2247 << kVoxelSpacing[2]
2248 << ") mm " << G4endl;
2249 }
2250
2251
2252 // offset from file starting point to the m
2253 ifile.read((char *)ctmp, 4);
2254 convertEndian(ctmp, kPointerToModalityData)
2255
2256 // # of dose distributions
2257 ifile.read((char *)ctmp, 4);
2258 int nDoseDist;
2259 convertEndian(ctmp, nDoseDist);
2260
2261 // offset from file starting point to the d
2262 for(int i = 0; i < nDoseDist; i++) {
2263 ifile.read((char *)ctmp, 4);
2264 unsigned int dptr;
2265 convertEndian(ctmp, dptr);
2266 addPointerToDoseDistData(dptr);
2267 }
2268
2269 // offset from file starting point to the R
2270 ifile.read((char *)ctmp, 4);
2271 convertEndian(ctmp, kPointerToROIData);
2272
2273 // offset from file starting point to the t
2274 ifile.read((char *)ctmp, 4);
2275 convertEndian(ctmp, kPointerToTrackData);
2276 if(DEBUG || kVerbose > 0) {
2277 G4cout << "Each pointer to data : "
2278 << kPointerToModalityData << ", ";
2279 for(int i = 0; i < nDoseDist; i++)
2280 G4cout << kPointerToDoseDistData[0] <<
2281 G4cout << kPointerToROIData << ", "
2282 << kPointerToTrackData << G4endl;
2283 }
2284
2285 if(kPointerToModalityData == 0 && kPointerT
2286 kPointerToROIData == 0 && kPointerToTrac
2287 if(DEBUG || kVerbose > 0) {
2288 G4cout << "No data." << G4endl;
2289 }
2290 return false;
2291 }
2292
2293 // event number
2294 /* ver 1
2295 ifile.read(ctmp, sizeof(int));
2296 convertEndian(ctmp, numberOfEvents);
2297 */
2298
2299 int size[3];
2300 float scale;
2301 double dscale;
2302 short minmax[2];
2303 float fCenter[3];
2304 int iCenter[3];
2305
2306 //----- Modality image -----//
2307 // modality image size
2308 ifile.read(ctmp, 3*sizeof(int));
2309 convertEndian(ctmp, size[0]);
2310 convertEndian(ctmp+sizeof(int), size[1]);
2311 convertEndian(ctmp+2*sizeof(int), size[2]);
2312 if(DEBUG || kVerbose > 0) {
2313 G4cout << "Modality image size : ("
2314 << size[0] << ", "
2315 << size[1] << ", "
2316 << size[2] << ")"
2317 << G4endl;
2318 }
2319 kModality.setSize(size);
2320
2321 // modality image voxel spacing
2322 /*
2323 ifile.read(ctmp, 3*sizeof(float));
2324 convertEndian(ctmp, modalityImageVoxelSpa
2325 convertEndian(ctmp+sizeof(float), modalit
2326 convertEndian(ctmp+2*sizeof(float), modal
2327 */
2328
2329 if(kPointerToModalityData != 0) {
2330
2331 // modality density max. & min.
2332 ifile.read((char *)ctmp, 4);
2333 convertEndian(ctmp, minmax[0]);
2334 convertEndian(ctmp+2, minmax[1]);
2335 kModality.setMinMax(minmax);
2336
2337 // modality image unit
2338 char munit[13];
2339 ifile.read((char *)munit, 12);
2340 std::string smunit = munit;
2341 setModalityImageUnit(smunit);
2342
2343 // modality density scale
2344 ifile.read((char *)ctmp, 4);
2345 convertEndian(ctmp, scale);
2346 kModality.setScale(dscale = scale);
2347 if(DEBUG || kVerbose > 0) {
2348 G4cout << "Modality image min., max., s
2349 << minmax[0] << ", "
2350 << minmax[1] << ", "
2351 << scale << G4endl;
2352 }
2353
2354 // modality density
2355 int psize = size[0]*size[1];
2356 if(DEBUG || kVerbose > 0) G4cout << "Moda
2357 char * cimage = new char[psize*sizeof(sho
2358 for(int i = 0; i < size[2]; i++) {
2359 ifile.read((char *)cimage, psize*sizeof
2360 short * mimage = new short[psize];
2361 for(int j = 0; j < psize; j++) {
2362 convertEndian(cimage+j*sizeof(short), mimag
2363 }
2364 kModality.addImage(mimage);
2365
2366 if(DEBUG || kVerbose > 0) G4cout << "["
2367 }
2368 if(DEBUG || kVerbose > 0) G4cout << G4end
2369 delete [] cimage;
2370
2371 // modality desity map for CT value
2372 size_t msize = minmax[1]-minmax[0]+1;
2373 if(DEBUG || kVerbose > 0) G4cout << "msiz
2374 char * pdmap = new char[msize*sizeof(floa
2375 ifile.read((char *)pdmap, msize*sizeof(fl
2376 float ftmp;
2377 for(int i = 0; i < (int)msize; i++) {
2378 convertEndian(pdmap+i*sizeof(float), ft
2379 kModalityImageDensityMap.push_back(ftmp
2380 }
2381 delete [] pdmap;
2382 if(DEBUG || kVerbose > 0) {
2383 G4cout << "density map : " << std::ends
2384 for(int i = 0; i < 10; i++)
2385 G4cout <<kModalityImageDensityMap[i] << ",
2386 G4cout << G4endl;
2387 for(int i = 0; i < 10; i++) G4cout << "
2388 G4cout << G4endl;
2389 for(size_t i =kModalityImageDensityMap.
2390 G4cout <<kModalityImageDensityMap[i] << ",
2391 G4cout << G4endl;
2392 }
2393
2394 }
2395
2396
2397 //----- dose distribution image -----//
2398 for(int ndose = 0; ndose < nDoseDist; ndose
2399
2400 newDoseDist();
2401
2402 // dose distrbution image size
2403 ifile.read((char *)ctmp, 3*sizeof(int));
2404 convertEndian(ctmp, size[0]);
2405 convertEndian(ctmp+sizeof(int), size[1]);
2406 convertEndian(ctmp+2*sizeof(int), size[2]
2407 if(DEBUG || kVerbose > 0) {
2408 G4cout << "Dose dist. image size : ("
2409 << size[0] << ", "
2410 << size[1] << ", "
2411 << size[2] << ")"
2412 << G4endl;
2413 }
2414 kDose[ndose].setSize(size);
2415
2416 // dose distribution max. & min.
2417 ifile.read((char *)ctmp, sizeof(short)*2)
2418 convertEndian(ctmp, minmax[0]);
2419 convertEndian(ctmp+2, minmax[1]);
2420
2421 // dose distribution unit
2422 char dunit[13];
2423 ifile.read((char *)dunit, 12);
2424 std::string sdunit = dunit;
2425 setDoseDistUnit(sdunit, ndose);
2426 if(DEBUG || kVerbose > 0) {
2427 G4cout << "Dose dist. unit : " << kDose
2428 }
2429
2430 // dose distribution scaling
2431 ifile.read((char *)ctmp, 4); // sizeof(fl
2432 convertEndian(ctmp, scale);
2433 kDose[ndose].setScale(dscale = scale);
2434
2435 double dminmax[2];
2436 for(int i = 0; i < 2; i++) dminmax[i] = m
2437 kDose[ndose].setMinMax(dminmax);
2438
2439 if(DEBUG || kVerbose > 0) {
2440 G4cout << "Dose dist. image min., max.,
2441 << dminmax[0] << ", "
2442 << dminmax[1] << ", "
2443 << scale << G4endl;
2444 }
2445
2446 // dose distribution image
2447 int dsize = size[0]*size[1];
2448 if(DEBUG || kVerbose > 0) G4cout << "Dose
2449 char * di = new char[dsize*sizeof(short)]
2450 short * shimage = new short[dsize];
2451 for(int z = 0; z < size[2]; z++) {
2452 ifile.read((char *)di, dsize*sizeof(sho
2453 double * dimage = new double[dsize];
2454 for(int xy = 0; xy < dsize; xy++) {
2455 convertEndian(di+xy*sizeof(short), shimage[
2456 dimage[xy] = shimage[xy]*dscale;
2457 }
2458 kDose[ndose].addImage(dimage);
2459
2460 if(DEBUG || kVerbose > 0) G4cout << "["
2461
2462 if(DEBUG || kVerbose > 0) {
2463 for(int j = 0; j < dsize; j++) {
2464 if(dimage[j] < 0)
2465 G4cout << "[" << j << "," << z << "]"
2466 << dimage[j] << ", ";
2467 }
2468 }
2469 }
2470 delete [] shimage;
2471 delete [] di;
2472 if(DEBUG || kVerbose > 0) G4cout << G4end
2473
2474 ifile.read((char *)ctmp, 3*4); // 3*sizeo
2475 convertEndian(ctmp, iCenter[0]);
2476 convertEndian(ctmp+4, iCenter[1]);
2477 convertEndian(ctmp+8, iCenter[2]);
2478 for(int i = 0; i < 3; i++) fCenter[i] = (
2479 kDose[ndose].setCenterPosition(fCenter);
2480
2481 if(DEBUG || kVerbose > 0) {
2482 G4cout << "Dose dist. image relative lo
2483 << fCenter[0] << ", "
2484 << fCenter[1] << ", "
2485 << fCenter[2] << ")" << G4endl;
2486 }
2487
2488
2489 }
2490
2491 //----- ROI image -----//
2492 if(kPointerToROIData != 0) {
2493
2494 newROI();
2495
2496 // ROI image size
2497 ifile.read((char *)ctmp, 3*sizeof(int));
2498 convertEndian(ctmp, size[0]);
2499 convertEndian(ctmp+sizeof(int), size[1]);
2500 convertEndian(ctmp+2*sizeof(int), size[2]
2501 kRoi[0].setSize(size);
2502 if(DEBUG || kVerbose > 0) {
2503 G4cout << "ROI image size : ("
2504 << size[0] << ", "
2505 << size[1] << ", "
2506 << size[2] << ")"
2507 << G4endl;
2508 }
2509
2510 // ROI max. & min.
2511 ifile.read((char *)ctmp, sizeof(short)*2)
2512 convertEndian(ctmp, minmax[0]);
2513 convertEndian(ctmp+sizeof(short), minmax[
2514 kRoi[0].setMinMax(minmax);
2515
2516 // ROI distribution scaling
2517 ifile.read((char *)ctmp, sizeof(float));
2518 convertEndian(ctmp, scale);
2519 kRoi[0].setScale(dscale = scale);
2520 if(DEBUG || kVerbose > 0) {
2521 G4cout << "ROI image min., max., scale
2522 << minmax[0] << ", "
2523 << minmax[1] << ", "
2524 << scale << G4endl;
2525 }
2526
2527 // ROI image
2528 int rsize = size[0]*size[1];
2529 char * ri = new char[rsize*sizeof(short)]
2530 for(int i = 0; i < size[2]; i++) {
2531 ifile.read((char *)ri, rsize*sizeof(sho
2532 short * rimage = new short[rsize];
2533 for(int j = 0; j < rsize; j++) {
2534 convertEndian(ri+j*sizeof(short), rimage[j]
2535 }
2536 kRoi[0].addImage(rimage);
2537
2538 }
2539 delete [] ri;
2540
2541 // ROI relative location
2542 ifile.read((char *)ctmp, 3*sizeof(int));
2543 convertEndian(ctmp, iCenter[0]);
2544 convertEndian(ctmp+sizeof(int), iCenter[1
2545 convertEndian(ctmp+2*sizeof(int), iCenter
2546 for(int i = 0; i < 3; i++) fCenter[i] = i
2547 kRoi[0].setCenterPosition(fCenter);
2548 if(DEBUG || kVerbose > 0) {
2549 G4cout << "ROI image relative location
2550 << fCenter[0] << ", "
2551 << fCenter[1] << ", "
2552 << fCenter[2] << ")" << G4endl;
2553 }
2554
2555 }
2556
2557 //----- track information -----//
2558 if(kPointerToTrackData != 0) {
2559
2560 // track
2561 ifile.read((char *)ctmp, sizeof(int));
2562 int ntrk;
2563 convertEndian(ctmp, ntrk);
2564 if(DEBUG || kVerbose > 0) {
2565 G4cout << "# of tracks: " << ntrk << G4
2566 }
2567
2568 // v4
2569 std::vector<float *> trkv4;
2570
2571 // track position
2572 for(int i = 0; i < ntrk; i++) {
2573 float * tp = new float[6];
2574
2575 ifile.read((char *)ctmp, sizeof(float)*
2576 if(DEBUG || kVerbose > 0) if(i < 10) G4
2577 for(int j = 0; j < 3; j++) {
2578 convertEndian(ctmp+j*sizeof(float), tp[j]);
2579 if(DEBUG || kVerbose > 0) if(i < 10) G4cout
2580 }
2581
2582 ifile.read((char *)ctmp, sizeof(float)*
2583 for(int j = 0; j < 3; j++) {
2584 convertEndian(ctmp+j*sizeof(float), tp[j+3]
2585 if(DEBUG || kVerbose > 0) if(i < 10) G4cout
2586 }
2587 addTrack(tp);
2588 if(DEBUG || kVerbose > 0) if(i < 10) G4
2589
2590 // v4
2591 trkv4.push_back(tp);
2592 }
2593
2594 //v4
2595 unsigned char trkcolorv4[3];
2596
2597 // track color
2598 for(int i = 0; i < ntrk; i++) {
2599 unsigned char * rgb = new unsigned char
2600 ifile.read((char *)rgb, 3);
2601 addTrackColor(rgb);
2602
2603 // v4
2604 for(int j = 0; j < 3; j++) trkcolorv4[j
2605 std::vector<float *> trk;
2606 trk.push_back(trkv4[i]);
2607 addTrack(trk, trkcolorv4);
2608
2609 }
2610
2611 }
2612
2613 ifile.close();
2614
2615 return true;
2616 }
2617 bool G4GMocrenIO::retrieveData3(char * _filen
2618 kFileName = _filename;
2619 return retrieveData();
2620 }
2621
2622 //
2623 bool G4GMocrenIO::retrieveData2() {
2624
2625 bool DEBUG = false;//
2626
2627 // input file open
2628 std::ifstream ifile(kFileName.c_str(), std:
2629 if(!ifile) {
2630 if (G4VisManager::GetVerbosity() >= G4Vis
2631 G4cout << "Cannot open file: " << kFile
2632 << " in G4GMocrenIO::retrieveData2()." <<
2633 return false;
2634 }
2635
2636 // data buffer
2637 char ctmp[12];
2638
2639 // file identifier
2640 char verid[9];
2641 ifile.read((char *)verid, 8);
2642
2643 // file version
2644 unsigned char ver;
2645 ifile.read((char *)&ver, 1);
2646 std::stringstream ss;
2647 ss << (int)ver;
2648 kVersion = ss.str();
2649 if(DEBUG || kVerbose > 0) G4cout << "File v
2650
2651 // id of version 1
2652 char idtmp[IDLENGTH];
2653 ifile.read((char *)idtmp, IDLENGTH);
2654 kId = idtmp;
2655 // version of version 1
2656 char vertmp[VERLENGTH];
2657 ifile.read((char *)vertmp, VERLENGTH);
2658
2659 // endian
2660 ifile.read((char *)&kLittleEndianInput, siz
2661 if(DEBUG || kVerbose > 0) {
2662 G4cout << "Endian : ";
2663 if(kLittleEndianInput == 1)
2664 G4cout << " little" << G4endl;
2665 else {
2666 G4cout << " big" << G4endl;
2667 }
2668 }
2669
2670 // voxel spacings for all images
2671 ifile.read((char *)ctmp, 12);
2672 convertEndian(ctmp, kVoxelSpacing[0]);
2673 convertEndian(ctmp+4, kVoxelSpacing[1]);
2674 convertEndian(ctmp+8, kVoxelSpacing[2]);
2675 if(DEBUG || kVerbose > 0) {
2676 G4cout << "Voxel spacing : ("
2677 << kVoxelSpacing[0] << ", "
2678 << kVoxelSpacing[1] << ", "
2679 << kVoxelSpacing[2]
2680 << ") mm " << G4endl;
2681 }
2682
2683
2684 // offset from file starting point to the m
2685 ifile.read((char *)ctmp, 4);
2686 convertEndian(ctmp, kPointerToModalityData)
2687
2688 // offset from file starting point to the d
2689 unsigned int ptddd;
2690 ifile.read((char *)ctmp, 4);
2691 convertEndian(ctmp, ptddd);
2692 kPointerToDoseDistData.push_back(ptddd);
2693
2694 // offset from file starting point to the R
2695 ifile.read((char *)ctmp, 4);
2696 convertEndian(ctmp, kPointerToROIData);
2697
2698 // offset from file starting point to the t
2699 ifile.read((char *)ctmp, 4);
2700 convertEndian(ctmp, kPointerToTrackData);
2701 if(DEBUG || kVerbose > 0) {
2702 G4cout << "Each pointer to data : "
2703 << kPointerToModalityData << ", "
2704 << kPointerToDoseDistData[0] << ", "
2705 << kPointerToROIData << ", "
2706 << kPointerToTrackData << G4endl;
2707 }
2708
2709 if(kPointerToModalityData == 0 && kPointerT
2710 kPointerToROIData == 0 && kPointerToTrac
2711 if(DEBUG || kVerbose > 0) {
2712 G4cout << "No data." << G4endl;
2713 }
2714 return false;
2715 }
2716
2717 // event number
2718 /* ver 1
2719 ifile.read(ctmp, sizeof(int));
2720 convertEndian(ctmp, numberOfEvents);
2721 */
2722
2723 int size[3];
2724 float scale;
2725 double dscale;
2726 short minmax[2];
2727 float fCenter[3];
2728 int iCenter[3];
2729
2730 //----- Modality image -----//
2731 // modality image size
2732 ifile.read(ctmp, 3*sizeof(int));
2733 convertEndian(ctmp, size[0]);
2734 convertEndian(ctmp+sizeof(int), size[1]);
2735 convertEndian(ctmp+2*sizeof(int), size[2]);
2736 if(DEBUG || kVerbose > 0) {
2737 G4cout << "Modality image size : ("
2738 << size[0] << ", "
2739 << size[1] << ", "
2740 << size[2] << ")"
2741 << G4endl;
2742 }
2743 kModality.setSize(size);
2744
2745 // modality image voxel spacing
2746 /*
2747 ifile.read(ctmp, 3*sizeof(float));
2748 convertEndian(ctmp, modalityImageVoxelSpa
2749 convertEndian(ctmp+sizeof(float), modalit
2750 convertEndian(ctmp+2*sizeof(float), modal
2751 */
2752
2753 if(kPointerToModalityData != 0) {
2754
2755 // modality density max. & min.
2756 ifile.read((char *)ctmp, 4);
2757 convertEndian(ctmp, minmax[0]);
2758 convertEndian(ctmp+2, minmax[1]);
2759 kModality.setMinMax(minmax);
2760
2761 // modality density scale
2762 ifile.read((char *)ctmp, 4);
2763 convertEndian(ctmp, scale);
2764 kModality.setScale(dscale = scale);
2765 if(DEBUG || kVerbose > 0) {
2766 G4cout << "Modality image min., max., s
2767 << minmax[0] << ", "
2768 << minmax[1] << ", "
2769 << scale << G4endl;
2770 }
2771
2772 // modality density
2773 int psize = size[0]*size[1];
2774 if(DEBUG || kVerbose > 0) G4cout << "Moda
2775 char * cimage = new char[psize*sizeof(sho
2776 for(int i = 0; i < size[2]; i++) {
2777 ifile.read((char *)cimage, psize*sizeof
2778 short * mimage = new short[psize];
2779 for(int j = 0; j < psize; j++) {
2780 convertEndian(cimage+j*sizeof(short), mimag
2781 }
2782 kModality.addImage(mimage);
2783
2784 if(DEBUG || kVerbose > 0) G4cout << "["
2785 }
2786 if(DEBUG || kVerbose > 0) G4cout << G4end
2787 delete [] cimage;
2788
2789 // modality desity map for CT value
2790 size_t msize = minmax[1]-minmax[0]+1;
2791 if(DEBUG || kVerbose > 0) G4cout << "msiz
2792 char * pdmap = new char[msize*sizeof(floa
2793 ifile.read((char *)pdmap, msize*sizeof(fl
2794 float ftmp;
2795 for(int i = 0; i < (int)msize; i++) {
2796 convertEndian(pdmap+i*sizeof(float), ft
2797 kModalityImageDensityMap.push_back(ftmp
2798 }
2799 delete [] pdmap;
2800 if(DEBUG || kVerbose > 0) {
2801 G4cout << "density map : " << std::ends
2802 for(int i = 0; i < 10; i++)
2803 G4cout <<kModalityImageDensityMap[i] << ",
2804 G4cout << G4endl;
2805 for(int i = 0; i < 10; i++) G4cout << "
2806 G4cout << G4endl;
2807 for(size_t i =kModalityImageDensityMap.
2808 G4cout <<kModalityImageDensityMap[i] << ",
2809 G4cout << G4endl;
2810 }
2811
2812 }
2813
2814
2815 //----- dose distribution image -----//
2816 if(kPointerToDoseDistData[0] != 0) {
2817
2818 newDoseDist();
2819
2820 // dose distrbution image size
2821 ifile.read((char *)ctmp, 3*sizeof(int));
2822 convertEndian(ctmp, size[0]);
2823 convertEndian(ctmp+sizeof(int), size[1]);
2824 convertEndian(ctmp+2*sizeof(int), size[2]
2825 if(DEBUG || kVerbose > 0) {
2826 G4cout << "Dose dist. image size : ("
2827 << size[0] << ", "
2828 << size[1] << ", "
2829 << size[2] << ")"
2830 << G4endl;
2831 }
2832 kDose[0].setSize(size);
2833
2834 // dose distribution max. & min.
2835 ifile.read((char *)ctmp, sizeof(short)*2)
2836 convertEndian(ctmp, minmax[0]);
2837 convertEndian(ctmp+2, minmax[1]);
2838 // dose distribution scaling
2839 ifile.read((char *)ctmp, sizeof(float));
2840 convertEndian(ctmp, scale);
2841 kDose[0].setScale(dscale = scale);
2842
2843 double dminmax[2];
2844 for(int i = 0; i < 2; i++) dminmax[i] = m
2845 kDose[0].setMinMax(dminmax);
2846
2847 if(DEBUG || kVerbose > 0) {
2848 G4cout << "Dose dist. image min., max.,
2849 << dminmax[0] << ", "
2850 << dminmax[1] << ", "
2851 << scale << G4endl;
2852 }
2853
2854 // dose distribution image
2855 int dsize = size[0]*size[1];
2856 if(DEBUG || kVerbose > 0) G4cout << "Dose
2857 char * di = new char[dsize*sizeof(short)]
2858 short * shimage = new short[dsize];
2859 for(int z = 0; z < size[2]; z++) {
2860 ifile.read((char *)di, dsize*sizeof(sho
2861 double * dimage = new double[dsize];
2862 for(int xy = 0; xy < dsize; xy++) {
2863 convertEndian(di+xy*sizeof(short), shimage[
2864 dimage[xy] = shimage[xy]*dscale;
2865 }
2866 kDose[0].addImage(dimage);
2867
2868 if(DEBUG || kVerbose > 0) G4cout << "["
2869
2870 if(DEBUG || kVerbose > 0) {
2871 for(int j = 0; j < dsize; j++) {
2872 if(dimage[j] < 0)
2873 G4cout << "[" << j << "," << z << "]"
2874 << dimage[j] << ", ";
2875 }
2876 }
2877 }
2878 delete [] shimage;
2879 delete [] di;
2880 if(DEBUG || kVerbose > 0) G4cout << G4end
2881
2882 /* ver 1
2883 float doseDist;
2884 int dosePid;
2885 double * doseData = new double[numDose
2886 for(int i = 0; i < numDose; i++) {
2887 ifile.read(ctmp, sizeof(int));
2888 convertEndian(ctmp, dosePid);
2889 for(int j = 0; j < numDoseImageVoxels;
2890 ifile.read(ctmp, sizeof(float));
2891 convertEndian(ctmp, doseDist);
2892 doseData[j] = doseDist;
2893 }
2894 setDose(dosePid, doseData);
2895 }
2896 delete [] doseData;
2897 if(totalDose == NULL) totalDose = new
2898 for(int i = 0; i < numDoseImageVoxels;
2899 ifile.read(ctmp, sizeof(float));
2900 convertEndian(ctmp, doseDist);
2901 totalDose[i] = doseDist;
2902 }
2903 */
2904
2905 /* ver 1
2906 // relative location between the two imag
2907 ifile.read(ctmp, 3*sizeof(float));
2908 convertEndian(ctmp, relativeLocation[0]);
2909 convertEndian(ctmp+sizeof(float), relativ
2910 convertEndian(ctmp+2*sizeof(float), relat
2911 */
2912
2913 // relative location of the dose distribu
2914 // the modality image
2915 //ofile.write((char *)relativeLocation, 3
2916 ifile.read((char *)ctmp, 3*sizeof(int));
2917 convertEndian(ctmp, iCenter[0]);
2918 convertEndian(ctmp+sizeof(int), iCenter[1
2919 convertEndian(ctmp+2*sizeof(int), iCenter
2920 for(int i = 0; i < 3; i++) fCenter[i] = (
2921 kDose[0].setCenterPosition(fCenter);
2922
2923 if(DEBUG || kVerbose > 0) {
2924 G4cout << "Dose dist. image relative lo
2925 << fCenter[0] << ", "
2926 << fCenter[1] << ", "
2927 << fCenter[2] << ")" << G4endl;
2928 }
2929
2930
2931 }
2932
2933 //----- ROI image -----//
2934 if(kPointerToROIData != 0) {
2935
2936 newROI();
2937
2938 // ROI image size
2939 ifile.read((char *)ctmp, 3*sizeof(int));
2940 convertEndian(ctmp, size[0]);
2941 convertEndian(ctmp+sizeof(int), size[1]);
2942 convertEndian(ctmp+2*sizeof(int), size[2]
2943 kRoi[0].setSize(size);
2944 if(DEBUG || kVerbose > 0) {
2945 G4cout << "ROI image size : ("
2946 << size[0] << ", "
2947 << size[1] << ", "
2948 << size[2] << ")"
2949 << G4endl;
2950 }
2951
2952 // ROI max. & min.
2953 ifile.read((char *)ctmp, sizeof(short)*2)
2954 convertEndian(ctmp, minmax[0]);
2955 convertEndian(ctmp+sizeof(short), minmax[
2956 kRoi[0].setMinMax(minmax);
2957
2958 // ROI distribution scaling
2959 ifile.read((char *)ctmp, sizeof(float));
2960 convertEndian(ctmp, scale);
2961 kRoi[0].setScale(dscale = scale);
2962 if(DEBUG || kVerbose > 0) {
2963 G4cout << "ROI image min., max., scale
2964 << minmax[0] << ", "
2965 << minmax[1] << ", "
2966 << scale << G4endl;
2967 }
2968
2969 // ROI image
2970 int rsize = size[0]*size[1];
2971 char * ri = new char[rsize*sizeof(short)]
2972 for(int i = 0; i < size[2]; i++) {
2973 ifile.read((char *)ri, rsize*sizeof(sho
2974 short * rimage = new short[rsize];
2975 for(int j = 0; j < rsize; j++) {
2976 convertEndian(ri+j*sizeof(short), rimage[j]
2977 }
2978 kRoi[0].addImage(rimage);
2979
2980 }
2981 delete [] ri;
2982
2983 // ROI relative location
2984 ifile.read((char *)ctmp, 3*sizeof(int));
2985 convertEndian(ctmp, iCenter[0]);
2986 convertEndian(ctmp+sizeof(int), iCenter[1
2987 convertEndian(ctmp+2*sizeof(int), iCenter
2988 for(int i = 0; i < 3; i++) fCenter[i] = i
2989 kRoi[0].setCenterPosition(fCenter);
2990 if(DEBUG || kVerbose > 0) {
2991 G4cout << "ROI image relative location
2992 << fCenter[0] << ", "
2993 << fCenter[1] << ", "
2994 << fCenter[2] << ")" << G4endl;
2995 }
2996
2997 }
2998
2999 //----- track information -----//
3000 if(kPointerToTrackData != 0) {
3001
3002 // track
3003 ifile.read((char *)ctmp, sizeof(int));
3004 int ntrk;
3005 convertEndian(ctmp, ntrk);
3006 if(DEBUG || kVerbose > 0) {
3007 G4cout << "# of tracks: " << ntrk << G4
3008 }
3009
3010 //v4
3011 unsigned char trkcolorv4[3] = {255, 0, 0}
3012
3013 for(int i = 0; i < ntrk; i++) {
3014 float * tp = new float[6];
3015 // v4
3016 std::vector<float *> trkv4;
3017
3018 ifile.read((char *)ctmp, sizeof(float)*
3019 if(DEBUG || kVerbose > 0) if(i < 10) G4
3020 for(int j = 0; j < 3; j++) {
3021 convertEndian(ctmp+j*sizeof(float), tp[j]);
3022 if(DEBUG || kVerbose > 0) if(i < 10) G4cout
3023 }
3024
3025 ifile.read((char *)ctmp, sizeof(float)*
3026 for(int j = 0; j < 3; j++) {
3027 convertEndian(ctmp+j*sizeof(float), tp[j+3]
3028 if(DEBUG || kVerbose > 0) if(i < 10) G4cout
3029 }
3030
3031 kSteps.push_back(tp);
3032 // v4
3033 trkv4.push_back(tp);
3034 addTrack(trkv4, trkcolorv4);
3035
3036 if(DEBUG || kVerbose > 0) if(i < 10) G4
3037 }
3038
3039 }
3040
3041 /* ver 1
3042 // track
3043 int ntracks;
3044 ifile.read(ctmp, sizeof(int));
3045 convertEndian(ctmp, ntracks);
3046 // track displacement
3047 ifile.read(ctmp, 3*sizeof(float));
3048 convertEndian(ctmp, trackDisplacement[0]);
3049 convertEndian(ctmp+sizeof(float), trackDisp
3050 convertEndian(ctmp+2*sizeof(float), trackDi
3051 //
3052 //for(int i = 0; i < ntracks && i < 100; i+
3053 for(int i = 0; i < ntracks; i++) {
3054 DicomDoseTrack trk;
3055 short trackid, parentid, pid;
3056 int npoints;
3057 ifile.read(ctmp, sizeof(short));
3058 convertEndian(ctmp, trackid);
3059 trk.setID(trackid);
3060 ifile.read(ctmp, sizeof(short));
3061 convertEndian(ctmp, parentid);
3062 trk.setParentID(parentid);
3063 ifile.read(ctmp, sizeof(short));
3064 convertEndian(ctmp, pid);
3065 trk.setPID(pid);
3066 ifile.read(ctmp, sizeof(int));
3067 convertEndian(ctmp, npoints);
3068 for(int i = 0; i < npoints; i++) {
3069 ifile.read(ctmp, 3*sizeof(float));
3070 // storing only start and end points
3071 //if(i == 0 || i == npoints - 1) {
3072 float * point = new float[3];
3073 convertEndian(ctmp, point[0]);
3074 convertEndian(ctmp+sizeof(float), point[1])
3075 convertEndian(ctmp+2*sizeof(float), point[2
3076 trk.addPoint(point);
3077 //}
3078 }
3079 track.push_back(trk);
3080 }
3081 */
3082
3083 ifile.close();
3084
3085 return true;
3086 }
3087
3088 bool G4GMocrenIO::retrieveData2(char * _filen
3089 kFileName = _filename;
3090 return retrieveData();
3091 }
3092
3093 void G4GMocrenIO::setID() {
3094 time_t t;
3095 time(&t);
3096
3097 tm * ti;
3098 ti = localtime(&t);
3099
3100 char cmonth[12][4] = {"Jan", "Feb", "Mar",
3101 "May", "Jun", "Jul", "Aug",
3102 "Sep", "Oct", "Nov", "Dec"};
3103 std::stringstream ss;
3104 ss << std::setfill('0')
3105 << std::setw(2)
3106 << ti->tm_hour << ":"
3107 << std::setw(2)
3108 << ti->tm_min << ":"
3109 << std::setw(2)
3110 << ti->tm_sec << ","
3111 << cmonth[ti->tm_mon] << "."
3112 << std::setw(2)
3113 << ti->tm_mday << ","
3114 << ti->tm_year+1900;
3115
3116 kId = ss.str();
3117 }
3118
3119 // get & set the file version
3120 std::string & G4GMocrenIO::getVersion() {retu
3121 void G4GMocrenIO::setVersion(std::string & _v
3122
3123 // set endians of input/output data
3124 void G4GMocrenIO::setLittleEndianInput(bool _
3125 void G4GMocrenIO::setLittleEndianOutput(bool
3126
3127 // voxel spacing
3128 void G4GMocrenIO::setVoxelSpacing(float _spac
3129 for(int i = 0; i < 3; i++) kVoxelSpacing[i]
3130 }
3131 void G4GMocrenIO::getVoxelSpacing(float _spac
3132 for(int i = 0; i < 3; i++) _spacing[i] = kV
3133 }
3134
3135 // get & set number of events
3136 int & G4GMocrenIO::getNumberOfEvents() {
3137 return kNumberOfEvents;
3138 }
3139 void G4GMocrenIO::setNumberOfEvents(int & _nu
3140 kNumberOfEvents = _numberOfEvents;
3141 }
3142 void G4GMocrenIO::addOneEvent() {
3143 kNumberOfEvents++;
3144 }
3145
3146 // set/get pointer the modality image data
3147 void G4GMocrenIO::setPointerToModalityData(un
3148 kPointerToModalityData = _pointer;
3149 }
3150 unsigned int G4GMocrenIO::getPointerToModalit
3151 return kPointerToModalityData;
3152 }
3153 // set/get pointer the dose distribution imag
3154 void G4GMocrenIO::addPointerToDoseDistData(un
3155 kPointerToDoseDistData.push_back(_pointer);
3156 }
3157 unsigned int G4GMocrenIO::getPointerToDoseDis
3158 if(kPointerToDoseDistData.size() == 0 ||
3159 kPointerToDoseDistData.size() < (size_t)
3160 return 0;
3161 else
3162 return kPointerToDoseDistData[_elem];
3163 }
3164
3165 // set/get pointer the ROI image data
3166 void G4GMocrenIO::setPointerToROIData(unsigne
3167 kPointerToROIData = _pointer;
3168 }
3169 unsigned int G4GMocrenIO::getPointerToROIData
3170 return kPointerToROIData;
3171 }
3172 // set/get pointer the track data
3173 void G4GMocrenIO::setPointerToTrackData(unsig
3174 kPointerToTrackData = _pointer;
3175 }
3176 unsigned int G4GMocrenIO::getPointerToTrackDa
3177 return kPointerToTrackData;
3178 }
3179
3180 // calculate pointers for version 4
3181 void G4GMocrenIO::calcPointers4() {
3182
3183 // pointer to modality data
3184 unsigned int pointer = 1070; // up to "poin
3185 int nDoseDist = getNumDoseDist();
3186 pointer += nDoseDist*4;
3187
3188 setPointerToModalityData(pointer);
3189
3190 // pointer to dose data
3191 // ct-density map for modality data
3192 int msize[3];
3193 getModalityImageSize(msize);
3194 short mminmax[2];
3195 getModalityImageMinMax(mminmax);
3196 int pmsize = 2*msize[0]*msize[1]*msize[2];
3197 int pmmap = 4*(mminmax[1] - mminmax[0] + 1)
3198 pointer += 32 + pmsize + pmmap;
3199 //
3200 kPointerToDoseDistData.clear();
3201 if(nDoseDist == 0) {
3202 unsigned int pointer0 = 0;
3203 addPointerToDoseDistData(pointer0);
3204 }
3205 for(int ndose = 0; ndose < nDoseDist; ndose
3206 addPointerToDoseDistData(pointer);
3207 int dsize[3];
3208 getDoseDistSize(dsize);
3209 pointer += 44 + dsize[0]*dsize[1]*dsize[2
3210 }
3211
3212 // pointer to roi data
3213 if(!isROIEmpty()) {
3214 setPointerToROIData(pointer);
3215
3216 int rsize[3];
3217 getROISize(rsize);
3218 int prsize = 2*rsize[0]*rsize[1]*rsize[2]
3219 pointer += 20 + prsize + 12;
3220 } else {
3221 unsigned int pointer0 = 0;
3222 setPointerToROIData(pointer0);
3223 }
3224
3225 // pointer to track data
3226 int ntrk = (int)kTracks.size();
3227 if(ntrk != 0) {
3228 setPointerToTrackData(pointer);
3229
3230 pointer += 4; // # of tracks
3231 for(int nt = 0; nt < ntrk; nt++) {
3232 int nsteps = kTracks[nt].getNumberOfSte
3233 pointer += 4 + 3 + nsteps*(4*6); // # o
3234 }
3235 } else {
3236 unsigned int pointer0 = 0;
3237 setPointerToTrackData(pointer0);
3238 }
3239 if(kVerbose > 0) G4cout << " pointer to the
3240 << kPointerToTrackData << G4endl;
3241
3242 // pointer to detector data
3243 int ndet = (int)kDetectors.size();
3244 if(ndet != 0) {
3245 kPointerToDetectorData = pointer;
3246 } else {
3247 kPointerToDetectorData = 0;
3248 }
3249 if(kVerbose > 0) G4cout << " pointer to the
3250 << kPointerToDetectorData << G4end
3251
3252 }
3253
3254 // calculate pointers for ver.3
3255 void G4GMocrenIO::calcPointers3() {
3256
3257 // pointer to modality data
3258 unsigned int pointer = 1066; // up to "poin
3259 int nDoseDist = getNumDoseDist();
3260 pointer += nDoseDist*4;
3261
3262 setPointerToModalityData(pointer);
3263
3264 // pointer to dose data
3265 // ct-density map for modality data
3266 int msize[3];
3267 getModalityImageSize(msize);
3268 short mminmax[2];
3269 getModalityImageMinMax(mminmax);
3270 int pmsize = 2*msize[0]*msize[1]*msize[2];
3271 int pmmap = 4*(mminmax[1] - mminmax[0] + 1)
3272 pointer += 32 + pmsize + pmmap;
3273 //
3274 kPointerToDoseDistData.clear();
3275 if(nDoseDist == 0) {
3276 unsigned int pointer0 = 0;
3277 addPointerToDoseDistData(pointer0);
3278 }
3279 for(int ndose = 0; ndose < nDoseDist; ndose
3280 addPointerToDoseDistData(pointer);
3281 int dsize[3];
3282 getDoseDistSize(dsize);
3283 pointer += 44 + dsize[0]*dsize[1]*dsize[2
3284 }
3285
3286 // pointer to roi data
3287 if(!isROIEmpty()) {
3288 setPointerToROIData(pointer);
3289
3290 int rsize[3];
3291 getROISize(rsize);
3292 int prsize = 2*rsize[0]*rsize[1]*rsize[2]
3293 pointer += 20 + prsize + 12;
3294 } else {
3295 unsigned int pointer0 = 0;
3296 setPointerToROIData(pointer0);
3297 }
3298
3299 //
3300 if(getNumTracks() != 0)
3301 setPointerToTrackData(pointer);
3302 else {
3303 unsigned int pointer0 = 0;
3304 setPointerToTrackData(pointer0);
3305 }
3306
3307 }
3308
3309 // calculate pointers for ver.2
3310 void G4GMocrenIO::calcPointers2() {
3311
3312 // pointer to modality data
3313 unsigned int pointer = 65;
3314 setPointerToModalityData(pointer);
3315
3316 // pointer to dose data
3317 int msize[3];
3318 getModalityImageSize(msize);
3319 short mminmax[2];
3320 getModalityImageMinMax(mminmax);
3321 int pmsize = 2*msize[0]*msize[1]*msize[2];
3322 int pmmap = 4*(mminmax[1] - mminmax[0] + 1)
3323 pointer += 20 + pmsize + pmmap;
3324 int dsize[3];
3325 getDoseDistSize(dsize);
3326 kPointerToDoseDistData.clear();
3327 if(dsize[0] != 0) {
3328 kPointerToDoseDistData.push_back(pointer)
3329
3330 int pdsize = 2*dsize[0]*dsize[1]*dsize[2]
3331 pointer += 20 + pdsize + 12;
3332 } else {
3333 unsigned int pointer0 = 0;
3334 kPointerToDoseDistData.push_back(pointer0
3335 }
3336
3337 // pointer to roi data
3338 if(!isROIEmpty()) {
3339 int rsize[3];
3340 getROISize(rsize);
3341 setPointerToROIData(pointer);
3342 int prsize = 2*rsize[0]*rsize[1]*rsize[2]
3343 pointer += 20 + prsize + 12;
3344
3345 } else {
3346 unsigned int pointer0 = 0;
3347 setPointerToROIData(pointer0);
3348 }
3349
3350 //
3351 if(getNumTracks() != 0)
3352 setPointerToTrackData(pointer);
3353 else {
3354 unsigned int pointer0 = 0;
3355 setPointerToTrackData(pointer0);
3356 }
3357
3358 }
3359
3360
3361 //----- Modality image -----//
3362 void G4GMocrenIO::getModalityImageSize(int _s
3363
3364 kModality.getSize(_size);
3365 }
3366 void G4GMocrenIO::setModalityImageSize(int _s
3367
3368 kModality.setSize(_size);
3369 }
3370
3371 // get & set the modality image size
3372 void G4GMocrenIO::setModalityImageScale(doubl
3373
3374 kModality.setScale(_scale);
3375 }
3376 double G4GMocrenIO::getModalityImageScale() {
3377
3378 return kModality.getScale();
3379 }
3380
3381 // set the modality image in CT
3382 void G4GMocrenIO::setModalityImage(short * _i
3383
3384 kModality.addImage(_image);
3385 }
3386 short * G4GMocrenIO::getModalityImage(int _z)
3387
3388 return kModality.getImage(_z);
3389 }
3390 void G4GMocrenIO::clearModalityImage() {
3391
3392 kModality.clearImage();
3393 }
3394 // set/get the modality image density map
3395 void G4GMocrenIO::setModalityImageDensityMap(
3396 kModalityImageDensityMap = _map;
3397 }
3398 std::vector<float> & G4GMocrenIO::getModality
3399 return kModalityImageDensityMap;
3400 }
3401 // set the modality image min./max.
3402 void G4GMocrenIO::setModalityImageMinMax(shor
3403
3404 kModality.setMinMax(_minmax);
3405 }
3406 // get the modality image min./max.
3407 void G4GMocrenIO::getModalityImageMinMax(shor
3408
3409 short minmax[2];
3410 kModality.getMinMax(minmax);
3411 for(int i = 0; i < 2; i++) _minmax[i] = min
3412 }
3413 short G4GMocrenIO::getModalityImageMax() {
3414
3415 short minmax[2];
3416 kModality.getMinMax(minmax);
3417 return minmax[1];
3418 }
3419 short G4GMocrenIO::getModalityImageMin() {
3420
3421 short minmax[2];
3422 kModality.getMinMax(minmax);
3423 return minmax[0];
3424 }
3425 // set/get position of the modality image cen
3426 void G4GMocrenIO::setModalityCenterPosition(f
3427
3428 kModality.setCenterPosition(_center);
3429 }
3430 void G4GMocrenIO::getModalityCenterPosition(f
3431
3432 if(isROIEmpty())
3433 for(int i = 0; i < 3; i++) _center[i] = 0
3434 else
3435 kModality.getCenterPosition(_center);
3436 }
3437 // get & set the modality image unit
3438 std::string G4GMocrenIO::getModalityImageUnit
3439 return kModalityUnit;
3440 }
3441 void G4GMocrenIO::setModalityImageUnit(std::s
3442 kModalityUnit = _unit;
3443 }
3444 //
3445 short G4GMocrenIO::convertDensityToHU(float &
3446 short rval = -1024; // default: air
3447 int nmap = (int)kModalityImageDensityMap.si
3448 if(nmap != 0) {
3449 short minmax[2];
3450 kModality.getMinMax(minmax);
3451 rval = minmax[1];
3452 for(int i = 0; i < nmap; i++) {
3453 //G4cout << kModalityImageDensityMap[i]
3454 if(_dens <= kModalityImageDensityMap[i]
3455 rval = i + minmax[0];
3456 break;
3457 }
3458 }
3459 }
3460 return rval;
3461 }
3462
3463
3464 //----- Dose distribution -----//
3465 //
3466 void G4GMocrenIO::newDoseDist() {
3467 GMocrenDataPrimitive<double> doseData;
3468 kDose.push_back(doseData);
3469 }
3470 int G4GMocrenIO::getNumDoseDist() {
3471 return (int)kDose.size();
3472 }
3473
3474 // get & set the dose distribution unit
3475 std::string G4GMocrenIO::getDoseDistUnit(int
3476 // to avoid a warning in the compile proces
3477 if(kDoseUnit.size() > static_cast<size_t>(_
3478
3479 return kDoseUnit;
3480 }
3481 void G4GMocrenIO::setDoseDistUnit(std::string
3482 // to avoid a warning in the compile proces
3483 if(_unit.size() > static_cast<size_t>(_num)
3484
3485 //char unit[13];
3486 //std::strncpy(unit, _unit.c_str(), 12);
3487 //doseUnit = unit;
3488 kDoseUnit = _unit;
3489 }
3490 //
3491 void G4GMocrenIO::getDoseDistSize(int _size[3
3492 if(isDoseEmpty())
3493 for(int i = 0; i < 3; i++) _size[i] = 0;
3494 else
3495 kDose[_num].getSize(_size);
3496 }
3497 void G4GMocrenIO::setDoseDistSize(int _size[3
3498
3499 kDose[_num].setSize(_size);
3500
3501 //resetDose();
3502 }
3503
3504 void G4GMocrenIO::setDoseDistMinMax(short _mi
3505
3506 double minmax[2];
3507 double scale = kDose[_num].getScale();
3508 for(int i = 0; i < 2; i++) minmax[i] = (dou
3509 kDose[_num].setMinMax(minmax);
3510 }
3511 void G4GMocrenIO::getDoseDistMinMax(short _mi
3512
3513 if(isDoseEmpty())
3514 for(int i = 0; i < 2; i++) _minmax[i] = 0
3515 else {
3516 double minmax[2];
3517 kDose[_num].getMinMax(minmax);
3518 double scale = kDose[_num].getScale();
3519 for(int i = 0; i < 2; i++) _minmax[i] = (
3520 }
3521 }
3522 void G4GMocrenIO::setDoseDistMinMax(double _m
3523
3524 kDose[_num].setMinMax(_minmax);
3525 }
3526 void G4GMocrenIO::getDoseDistMinMax(double _m
3527
3528 if(isDoseEmpty())
3529 for(int i = 0; i < 2; i++) _minmax[i] = 0
3530 else
3531 kDose[_num].getMinMax(_minmax);
3532 }
3533
3534 // get & set the dose distribution image scal
3535 void G4GMocrenIO::setDoseDistScale(double & _
3536
3537 kDose[_num].setScale(_scale);
3538 }
3539 double G4GMocrenIO::getDoseDistScale(int _num
3540
3541 if(isDoseEmpty())
3542 return 0.;
3543 else
3544 return kDose[_num].getScale();
3545 }
3546
3547 /*
3548 void G4GMocrenIO::initializeShortDoseDist()
3549 ;
3550 }
3551 void G4GMocrenIO::finalizeShortDoseDist() {
3552 ;
3553 }
3554 */
3555 // set the dose distribution image
3556 void G4GMocrenIO::setShortDoseDist(short * _i
3557
3558 int size[3];
3559 kDose[_num].getSize(size);
3560 int dsize = size[0]*size[1];
3561 double * ddata = new double[dsize];
3562 double scale = kDose[_num].getScale();
3563 double minmax[2];
3564 kDose[_num].getMinMax(minmax);
3565 for(int xy = 0; xy < dsize; xy++) {
3566 ddata[xy] = _image[xy]*scale;
3567 if(ddata[xy] < minmax[0]) minmax[0] = dda
3568 if(ddata[xy] > minmax[1]) minmax[1] = dda
3569 }
3570 kDose[_num].addImage(ddata);
3571
3572 // set min./max.
3573 kDose[_num].setMinMax(minmax);
3574 }
3575 void G4GMocrenIO::getShortDoseDist(short * _d
3576
3577 if(_data == NULL) {
3578 if (G4VisManager::GetVerbosity() >= G4Vis
3579 G4cout << "In G4GMocrenIO::getShortDose
3580 << "first argument is NULL pointer. "
3581 << "The argument must be allocated array.
3582 << G4endl;
3583 G4Exception("G4GMocrenIO::getShortDoseDis
3584 "gMocren2002", FatalException
3585 "Error.");
3586 return;
3587 }
3588
3589 int size[3];
3590 kDose[_num].getSize(size);
3591 //short * shdata = new short[size[0]*size[1
3592 double * ddata = kDose[_num].getImage(_z);
3593 double scale = kDose[_num].getScale();
3594 for(int xy = 0; xy < size[0]*size[1]; xy++)
3595 _data[xy] = (short)(ddata[xy]/scale+0.5);
3596 }
3597 }
3598 void G4GMocrenIO::getShortDoseDistMinMax(shor
3599 double scale = kDose[_num].getScale();
3600 double minmax[2];
3601 kDose[_num].getMinMax(minmax);
3602 for(int i = 0; i < 2; i++)
3603 _minmax[i] = (short)(minmax[i]/scale+0.5)
3604 }
3605 //
3606 void G4GMocrenIO::setDoseDist(double * _image
3607
3608 kDose[_num].addImage(_image);
3609 }
3610 double * G4GMocrenIO::getDoseDist(int _z, int
3611
3612 double * image;
3613 if(isDoseEmpty()) {
3614 image = 0;
3615 } else {
3616 image = kDose[_num].getImage(_z);
3617 }
3618 return image;
3619 }
3620 /*
3621 void G4GMocrenIO::getDoseDist(double * & _i
3622
3623 G4cout << " <" << (void*)_image << "> ";
3624 if(isDoseEmpty()) {
3625 _image = 0;
3626 } else {
3627 _image = kDose[_num].getImage(_z);
3628 G4cout << " <" << (void*)_image << "> ";
3629 G4cout << _image[100] << " ";
3630 }
3631 }
3632 */
3633 bool G4GMocrenIO::addDoseDist(std::vector<dou
3634
3635 int size[3];
3636 getDoseDistSize(size, _num);
3637 std::vector<double *> dosedist = kDose[_num
3638
3639 int nimg = size[0]*size[1];
3640 for(int z = 0; z < size[2]; z++) {
3641 for(int xy = 0; xy < nimg; xy++) {
3642 dosedist[z][xy] += _image[z][xy];
3643 }
3644 }
3645
3646 return true;
3647 }
3648 //void setDoseDistDensityMap(float * _map) {d
3649 // set the dose distribution image displaceme
3650 void G4GMocrenIO::setDoseDistCenterPosition(f
3651
3652 kDose[_num].setCenterPosition(_center);
3653 }
3654 void G4GMocrenIO::getDoseDistCenterPosition(f
3655
3656 if(isDoseEmpty())
3657 for(int i = 0; i < 3; i++) _center[i] = 0
3658 else
3659 kDose[_num].getCenterPosition(_center);
3660 }
3661 // set & get name of dose distribution
3662 void G4GMocrenIO::setDoseDistName(std::string
3663
3664 kDose[_num].setName(_name);
3665 }
3666 std::string G4GMocrenIO::getDoseDistName(int
3667
3668 std::string name;
3669 if(isDoseEmpty())
3670 return name;
3671 else
3672 return kDose[_num].getName();
3673 }
3674 // copy dose distributions
3675 void G4GMocrenIO::copyDoseDist(std::vector<cl
3676 std::vector<class GMocrenDataPrimitive<doub
3677 for(itr = kDose.begin(); itr != kDose.end()
3678 _dose.push_back(*itr);
3679 }
3680 }
3681 // merge two dose distributions
3682 bool G4GMocrenIO::mergeDoseDist(std::vector<c
3683 if(kDose.size() != _dose.size()) {
3684 if (G4VisManager::GetVerbosity() >= G4Vis
3685 G4cout << "G4GMocrenIO::mergeDoseDist()
3686 G4cout << " Unable to merge the dose
3687 G4cout << " because of different size
3688 }
3689 return false;
3690 }
3691
3692 int num = (int)kDose.size();
3693 std::vector<class GMocrenDataPrimitive<doub
3694 std::vector<class GMocrenDataPrimitive<doub
3695 for(int i = 0; i < num; i++, itr1++, itr2++
3696 if (G4VisManager::GetVerbosity() >= G4Vis
3697 if(kVerbose > 0)
3698 G4cout << "merged dose distribution [" << i
3699 *itr1 += *itr2;
3700 }
3701
3702 return true;
3703 }
3704 //
3705 void G4GMocrenIO::clearDoseDistAll() {
3706
3707 if(!isDoseEmpty()) {
3708 for(int i = 0; i < getNumDoseDist(); i++)
3709 kDose[i].clear();
3710 }
3711 kDose.clear();
3712 }
3713 }
3714 //
3715 bool G4GMocrenIO::isDoseEmpty() {
3716 if(kDose.empty()) {
3717 //if (G4VisManager::GetVerbosity() >= G4V
3718 // G4cout << "!!! dose distribution data
3719 return true;
3720 } else {
3721 return false;
3722 }
3723 }
3724
3725 //
3726 void G4GMocrenIO::calcDoseDistScale() {
3727
3728 double scale;
3729 double minmax[2];
3730
3731 for(int i = 0; i < (int)kDose.size(); i++)
3732 kDose[i].getMinMax(minmax);
3733 scale = minmax[1]/DOSERANGE;
3734 kDose[i].setScale(scale);
3735 }
3736 }
3737
3738
3739 //----- RoI -----//
3740
3741 // add one RoI data
3742 void G4GMocrenIO::newROI() {
3743 GMocrenDataPrimitive<short> roiData;
3744 kRoi.push_back(roiData);
3745 }
3746 int G4GMocrenIO::getNumROI() {
3747 return (int)kRoi.size();
3748 }
3749
3750 // set/get the ROI image scale
3751 void G4GMocrenIO::setROIScale(double & _scale
3752
3753 kRoi[_num].setScale(_scale);
3754 }
3755 double G4GMocrenIO::getROIScale(int _num) {
3756
3757 if(isROIEmpty())
3758 return 0.;
3759 else
3760 return kRoi[_num].getScale();
3761 }
3762 // set the ROI image
3763 void G4GMocrenIO::setROI(short * _image, int
3764
3765 kRoi[_num].addImage(_image);
3766 }
3767 short * G4GMocrenIO::getROI(int _z, int _num)
3768
3769 if(isROIEmpty())
3770 return 0;
3771 else
3772 return kRoi[_num].getImage(_z);
3773 }
3774 // set/get the ROI image size
3775 void G4GMocrenIO::setROISize(int _size[3], in
3776
3777 return kRoi[_num].setSize(_size);
3778 }
3779 void G4GMocrenIO::getROISize(int _size[3], in
3780
3781 if(isROIEmpty())
3782 for(int i = 0; i < 3; i++) _size[i] = 0;
3783 else
3784 return kRoi[_num].getSize(_size);
3785 }
3786 // set/get the ROI image min. and max.
3787 void G4GMocrenIO::setROIMinMax(short _minmax[
3788
3789 kRoi[_num].setMinMax(_minmax);
3790 }
3791 void G4GMocrenIO::getROIMinMax(short _minmax[
3792
3793 if(isROIEmpty())
3794 for(int i = 0; i < 2; i++) _minmax[i] = 0
3795 else
3796 kRoi[_num].getMinMax(_minmax);
3797 }
3798 // set/get the ROI image displacement
3799 void G4GMocrenIO::setROICenterPosition(float
3800
3801 kRoi[_num].setCenterPosition(_center);
3802 }
3803 void G4GMocrenIO::getROICenterPosition(float
3804
3805 if(isROIEmpty())
3806 for(int i = 0; i < 3; i++) _center[i] = 0
3807 else
3808 kRoi[_num].getCenterPosition(_center);
3809 }
3810 //
3811 void G4GMocrenIO::clearROIAll() {
3812
3813 if(!isROIEmpty()) {
3814 for(int i = 0; i < getNumROI(); i++) {
3815 kRoi[i].clear();
3816 }
3817 kRoi.clear();
3818 }
3819 }
3820 //
3821 bool G4GMocrenIO::isROIEmpty() {
3822 if(kRoi.empty()) {
3823 //if (G4VisManager::GetVerbosity() >= G4V
3824 // G4cout << "!!! ROI data is empty." <<
3825 return true;
3826 } else {
3827 return false;
3828 }
3829 }
3830
3831
3832
3833 //----- Track information -----//
3834
3835 int G4GMocrenIO::getNumTracks() {
3836 return (int)kSteps.size();
3837 }
3838 int G4GMocrenIO::getNumTracks4() {
3839 return (int)kTracks.size();
3840 }
3841 void G4GMocrenIO::addTrack(float * _tracks) {
3842 kSteps.push_back(_tracks);
3843 }
3844 void G4GMocrenIO::setTracks(std::vector<float
3845 kSteps = _tracks;
3846 }
3847 std::vector<float *> & G4GMocrenIO::getTracks
3848 return kSteps;
3849 }
3850 void G4GMocrenIO::addTrackColor(unsigned char
3851 kStepColors.push_back(_colors);
3852 }
3853 void G4GMocrenIO::setTrackColors(std::vector<
3854 kStepColors = _trackColors;
3855 }
3856 std::vector<unsigned char *> & G4GMocrenIO::g
3857 return kStepColors;
3858 }
3859 void G4GMocrenIO::copyTracks(std::vector<floa
3860 std::vector<unsigned char *> & _
3861 std::vector<float *>::iterator titr;
3862 for(titr = kSteps.begin(); titr != kSteps.e
3863 float * pts = new float[6];
3864 for(int i = 0; i < 6; i++) {
3865 pts[i] = (*titr)[i];
3866 }
3867 _tracks.push_back(pts);
3868 }
3869
3870 std::vector<unsigned char *>::iterator citr
3871 for(citr = kStepColors.begin(); citr != kSt
3872 unsigned char * pts = new unsigned char[3
3873 for(int i = 0; i < 3; i++) {
3874 pts[i] = (*citr)[i];
3875 }
3876 _colors.push_back(pts);
3877 }
3878 }
3879 void G4GMocrenIO::mergeTracks(std::vector<flo
3880 std::vector<unsigned char *> & _color
3881 std::vector<float *>::iterator titr;
3882 for(titr = _tracks.begin(); titr != _tracks
3883 addTrack(*titr);
3884 }
3885
3886 std::vector<unsigned char *>::iterator citr
3887 for(citr = _colors.begin(); citr != _colors
3888 addTrackColor(*citr);
3889 }
3890 }
3891 void G4GMocrenIO::addTrack(std::vector<float
3892
3893 std::vector<float *>::iterator itr = _steps
3894 std::vector<struct GMocrenTrack::Step> st
3895 for(; itr != _steps.end(); itr++) {
3896 struct GMocrenTrack::Step step;
3897 for(int i = 0; i < 3; i++) {
3898 step.startPoint[i] = (*itr)[i];
3899 step.endPoint[i] = (*itr)[i+3];
3900 }
3901 steps.push_back(step);
3902 }
3903 GMocrenTrack track;
3904 track.setTrack(steps);
3905 track.setColor(_color);
3906 kTracks.push_back(track);
3907
3908 }
3909 void G4GMocrenIO::getTrack(int _num, std::vec
3910 std::vector<unsigned char *> & _co
3911
3912 if(_num > (int)kTracks.size()) {
3913 if (G4VisManager::GetVerbosity() >= G4Vis
3914 G4cout << "ERROR in getTrack() : " << G
3915 G4Exception("G4GMocrenIO::getTrack()",
3916 "gMocren2003", FatalException
3917 "Error.");
3918 }
3919 unsigned char * color = new unsigned char[3
3920 kTracks[_num].getColor(color);
3921 _color.push_back(color);
3922
3923 // steps
3924 int nsteps = kTracks[_num].getNumberOfSteps
3925 for(int isteps = 0; isteps < nsteps; isteps
3926 float * stepPoints = new float[6];
3927 kTracks[_num].getStep(stepPoints[0], step
3928 stepPoints[3], stepPoints[4], stepPoi
3929 isteps);
3930 _steps.push_back(stepPoints);
3931 }
3932 }
3933
3934 void G4GMocrenIO::translateTracks(std::vector
3935 std::vector<class GMocrenTrack>::iterator i
3936 for(; itr != kTracks.end(); itr++) {
3937 itr->translate(_translate);
3938 }
3939 }
3940
3941
3942
3943
3944 //----- Detector information -----//
3945 int G4GMocrenIO::getNumberOfDetectors() {
3946 return (int)kDetectors.size();
3947 }
3948 void G4GMocrenIO::addDetector(std::string & _
3949 std::vector<float *> & _det,
3950 unsigned char _color[3]) {
3951
3952 std::vector<float *>::iterator itr = _det
3953 std::vector<struct GMocrenDetector::Edge>
3954 for(; itr != _det.end(); itr++) {
3955 struct GMocrenDetector::Edge edge;
3956 for(int i = 0; i < 3; i++) {
3957 edge.startPoint[i] = (*itr)[i];
3958 edge.endPoint[i] = (*itr)[i+3];
3959 }
3960 edges.push_back(edge);
3961 }
3962 GMocrenDetector detector;
3963 detector.setDetector(edges);
3964 detector.setColor(_color);
3965 detector.setName(_name);
3966 kDetectors.push_back(detector);
3967
3968 }
3969
3970 void G4GMocrenIO::getDetector(int _num, std::
3971 std::vector<unsigned char *> & _color
3972 std::string & _detName) {
3973
3974 if(_num > (int)kDetectors.size()) {
3975 if (G4VisManager::GetVerbosity() >= G4Vis
3976 G4cout << "ERROR in getDetector() : " <
3977
3978 G4Exception("G4GMocrenIO::getDetector()",
3979 "gMocren2004", FatalException
3980 "Error.");
3981 }
3982
3983 _detName = kDetectors[_num].getName();
3984
3985 unsigned char * color = new unsigned char[3
3986 kDetectors[_num].getColor(color);
3987 _color.push_back(color);
3988
3989 // edges
3990 int nedges = kDetectors[_num].getNumberOfEd
3991 for(int ne = 0; ne < nedges; ne++) {
3992 float * edgePoints = new float[6];
3993 kDetectors[_num].getEdge(edgePoints[0], e
3994 edgePoints[3], edgePoints[4], edge
3995 ne);
3996 _edges.push_back(edgePoints);
3997 }
3998 }
3999
4000 void G4GMocrenIO::translateDetector(std::vect
4001 std::vector<class GMocrenDetector>::iterato
4002 for(; itr != kDetectors.end(); itr++) {
4003 itr->translate(_translate);
4004 }
4005 }
4006
4007 // endian conversion
4008 template <typename T>
4009 void G4GMocrenIO::convertEndian(char * _val,
4010
4011 if((kLittleEndianOutput && !kLittleEndianIn
4012 (!kLittleEndianOutput && kLittleEndianIn
4013
4014 const int SIZE = sizeof(_rval);
4015 char ctemp;
4016 for(int i = 0; i < SIZE/2; i++) {
4017 ctemp = _val[i];
4018 _val[i] = _val[SIZE - 1 - i];
4019 _val[SIZE - 1 - i] = ctemp;
4020 }
4021 }
4022 _rval = *(T *)_val;
4023 }
4024
4025 // inversion of byte order
4026 template <typename T>
4027 void G4GMocrenIO::invertByteOrder(char * _val
4028
4029 const int SIZE = sizeof(_rval);
4030 //char * cval = new char[SIZE];
4031 union {
4032 char cu[16];
4033 T tu;
4034 } uni;
4035 for(int i = 0; i < SIZE; i++) {
4036 uni.cu[i] = _val[SIZE-1-i];
4037 //cval[i] = _val[SIZE-i-1];
4038 }
4039 //_rval = *(T *)cval;
4040 _rval = uni.tu;
4041 //delete [] cval;
4042 }
4043
4044 //----- kVerbose information -----//
4045 void G4GMocrenIO::setVerboseLevel(int _level)
4046 kVerbose = _level;
4047 }
4048
4049