ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4GMocrenIO.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4GMocrenIO.cc
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 //
27 //
28 //
29 // File I/O manager class for writing or reading calcuated dose
30 // distribution and some event information
31 //
32 // Created: Mar. 31, 2009 Akinori Kimura : release for the gMocrenFile driver
33 //
34 // Akinori Kimura
35 // gMocren home page:
36 // http://geant4.kek.jp/gMocren/
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 GMocrenDataIO class-----//
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.begin();
66  for(; itr != image.end(); itr++) {
67  delete [] *itr;
68  }
69  */
70 }
71 
72 template <typename T> GMocrenDataPrimitive<T> &
73 GMocrenDataPrimitive<T>::operator = (const GMocrenDataPrimitive<T> & _right) {
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] = _right.kMinmax[i];
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] =_right.kImage[z][i];
86  kImage.push_back(img);
87  }
88  return *this;
89 }
90 
91 template <typename T> GMocrenDataPrimitive<T> &
92 GMocrenDataPrimitive<T>::operator + (const GMocrenDataPrimitive<T> & _right) {
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 = false;
98  if(kCenter[i] != _right.kCenter[i]) stat = false;
99  }
100  if(!stat) {
101  if (G4VisManager::GetVerbosity() >= G4VisManager::errors)
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] = _right.minmax[0];
113  //if(mms[1] < _right.minmax[1]) mms[1] = _right.minmax[1];
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[z][xy];
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 GMocrenDataPrimitive<T> & _right) {
135 
136  bool stat = true;
137  for(int i = 0; i < 3; i++) {
138  if(kSize[i] != _right.kSize[i]) stat = false;
139  if(kCenter[i] != _right.kCenter[i]) stat = false;
140  }
141  if(!stat) {
142  if (G4VisManager::GetVerbosity() >= G4VisManager::errors)
143  G4cout << "Warning: operator += " << G4endl
144  << " Cannot do the operator +="
145  << G4endl;
146  return *this;
147  }
148 
149  if(kMinmax[0] > _right.kMinmax[0]) kMinmax[0] = _right.kMinmax[0];
150  if(kMinmax[1] < _right.kMinmax[1]) kMinmax[1] = _right.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] = kImage[z][xy];
157  if(kMinmax[1] < kImage[z][xy]) kMinmax[1] = kImage[z][xy];
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(); itr++) {
182  delete [] *itr;
183  }
184  kImage.clear();
185 }
186 template <typename T>
187 void GMocrenDataPrimitive<T>::setSize(int _size[3]) {
188  for(int i = 0; i < 3; i++) kSize[i] = _size[i];
189 }
190 template <typename T>
191 void GMocrenDataPrimitive<T>::getSize(int _size[3]) {
192  for(int i = 0; i < 3; i++) _size[i] = kSize[i];
193 }
194 template <typename T>
195 void GMocrenDataPrimitive<T>::setScale(double & _scale) {
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 _minmax[2]) {
204  for(int i = 0; i < 2; i++) kMinmax[i] = _minmax[i];
205 }
206 template <typename T>
207 void GMocrenDataPrimitive<T>::getMinMax(T _minmax[2]) {
208  for(int i = 0; i < 2; i++) _minmax[i] = kMinmax[i];
209 
210 }
211 template <typename T>
212 void GMocrenDataPrimitive<T>::setImage(std::vector<T *> & _image) {
213  kImage = _image;
214 }
215 template <typename T>
216 void GMocrenDataPrimitive<T>::addImage(T * _image) {
217  kImage.push_back(_image);
218 }
219 template <typename T>
220 std::vector<T *> & GMocrenDataPrimitive<T>::getImage() {
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>::setCenterPosition(float _center[3]) {
230  for(int i = 0; i < 3; i++) kCenter[i] = _center[i];
231 }
232 template <typename T>
233 void GMocrenDataPrimitive<T>::getCenterPosition(float _center[3]) {
234  for(int i = 0; i < 3; i++) _center[i] = kCenter[i];
235 }
236 template <typename T>
237 void GMocrenDataPrimitive<T>::setName(std::string & _name) {
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, float _starty, float _startz,
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, float & _starty, float & _startz,
266  float & _endx, float & _endy, float & _endz,
267  int _num) {
268  if(_num >= (int)kTrack.size()) {
269  if (G4VisManager::GetVerbosity() >= G4VisManager::errors)
270  G4cout << "GMocrenTrack::getStep(...) Error: "
271  << "invalid step # : " << _num << G4endl;
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> & _translate) {
283  std::vector<struct Step>::iterator itr = kTrack.begin();
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, float _starty, float _startz,
306  float _endx, float _endy, float _endz) {
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, float & _starty, float & _startz,
317  float & _endx, float & _endy, float & _endz,
318  int _num) {
319  if(_num >= (int)kDetector.size()) {
320  if (G4VisManager::GetVerbosity() >= G4VisManager::errors)
321  G4cout << "GMocrenDetector::getEdge(...) Error: "
322  << "invalid edge # : " << _num << G4endl;
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<float> & _translate) {
334  std::vector<struct Edge>::iterator itr = kDetector.begin();
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; // Big endian
361 #endif
362 std::string G4GMocrenIO::kComment;
363 //
364 std::string G4GMocrenIO::kFileName = "dose.gdd";
365 
366 //
367 unsigned int G4GMocrenIO::kPointerToModalityData = 0;
368 std::vector<unsigned int> G4GMocrenIO::kPointerToDoseDistData;
369 unsigned int G4GMocrenIO::kPointerToROIData = 0;
370 unsigned int G4GMocrenIO::kPointerToTrackData = 0;
371 unsigned int G4GMocrenIO::kPointerToDetectorData = 0;
372 
373 // modality
374 float G4GMocrenIO::kVoxelSpacing[3] = {0., 0., 0.};
375 class GMocrenDataPrimitive<short> G4GMocrenIO::kModality;
376 std::vector<float> G4GMocrenIO::kModalityImageDensityMap;
377 std::string G4GMocrenIO::kModalityUnit = "g/cm3 "; // 12 Bytes
378 
379 // dose
380 std::vector<class GMocrenDataPrimitive<double> > G4GMocrenIO::kDose;
381 std::string G4GMocrenIO::kDoseUnit = "keV "; // 12 Bytes
382 
383 // ROI
384 std::vector<class GMocrenDataPrimitive<short> > G4GMocrenIO::kRoi;
385 
386 // track
387 std::vector<float *> G4GMocrenIO::kSteps;
388 std::vector<unsigned char *> G4GMocrenIO::kStepColors;
389 std::vector<class GMocrenTrack> G4GMocrenIO::kTracks;
390 
391 // detector
392 std::vector<class GMocrenDetector> G4GMocrenIO::kDetectors;
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] = 0.;
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(); itr++) delete [] *itr;
442  kSteps.clear();
443  std::vector<unsigned char *>::iterator citr;
444  for(citr = kStepColors.begin(); citr != kStepColors.end(); citr++)
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) <<<<<<<" << G4endl;
467  if(DEBUG || kVerbose > 0)
468  G4cout << " " << kFileName << G4endl;
469 
470  // output file open
471  std::ofstream ofile(kFileName.c_str(),
472  std::ios_base::out|std::ios_base::binary);
473  if(DEBUG || kVerbose > 0)
474  G4cout << " file open status: " << ofile.rdbuf() << G4endl;
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, sizeof(char));
485  char littleEndian = 0x01;
486  ofile.write((char *)&littleEndian, sizeof(char));
487  if(DEBUG || kVerbose > 0) {
488  //G4cout << "Endian: " << (int)kLittleEndianOutput << G4endl;
489  G4cout << "Endian: " << (int)littleEndian << G4endl;
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, itmp[0]);
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], ftmp[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 = 0;
534 
535  // offset from file starting point to the modality image data
536  if(kLittleEndianOutput) {
537  ofile.write((char *)&kPointerToModalityData, 4);
538  } else {
539  invertByteOrder((char *)&kPointerToModalityData, itmp[0]);
540  ofile.write((char *)itmp, 4);
541  }
542 
543  // # of dose distributions
544  //int nDoseDist = (int)pointerToDoseDistData.size();
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 dose image data
554  if(kLittleEndianOutput) {
555  for(int i = 0; i < nDoseDist; i++) {
556  ofile.write((char *)&kPointerToDoseDistData[i], 4);
557  }
558  } else {
559  for(int i = 0; i < nDoseDist; i++) {
560  invertByteOrder((char *)&kPointerToDoseDistData[i], itmp[0]);
561  ofile.write((char *)itmp, 4);
562  }
563  }
564 
565  // offset from file starting point to the ROI image data
566  if(kLittleEndianOutput) {
567  ofile.write((char *)&kPointerToROIData, 4);
568  } else {
569  invertByteOrder((char *)&kPointerToROIData, itmp[0]);
570  ofile.write((char *)itmp, 4);
571  }
572 
573  // offset from file starting point to the track data
574  if(kLittleEndianOutput) {
575  ofile.write((char *)&kPointerToTrackData, 4);
576  } else {
577  invertByteOrder((char *)&kPointerToTrackData, itmp[0]);
578  ofile.write((char *)itmp, 4);
579  }
580 
581  // offset from file starting point to the detector data
582  if(kLittleEndianOutput) {
583  ofile.write((char *)&kPointerToDetectorData, 4);
584  } else {
585  invertByteOrder((char *)&kPointerToDetectorData, itmp[0]);
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[j]);
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., scale : "
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 << "Modality image : ";
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(short));
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 << "[" << i << "]" << image[(size_t)(psize*0.55)] << ", ";
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] << ", " << minmax[1] << G4endl;
678  float * pdmap = new float[msize];
679  for(int i = 0; i < (int)msize; i++) pdmap[i] =kModalityImageDensityMap[i];
680 
681  if(kLittleEndianOutput) {
682  ofile.write((char *)pdmap, msize*sizeof(float));
683  } else {
684  for(int j = 0; j < (int)msize; j++) {
685  invertByteOrder((char *)&pdmap[j], ftmp[0]);
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; ndose++) {
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 << "] image size : ("
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(shorft)*2
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(), 12);
736  cdunit[12] = '\0';
737  ofile.write(cdunit, 12);
738  if(DEBUG || kVerbose > 0) {
739  G4cout << "Dose dist. unit : " << kDoseUnit << G4endl;
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); //sizeof(short)
767  } else {
768  for(int j = 0; j < dsize; j++) {
769  invertByteOrder((char *)&dimage[j], stmp[0]);
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 << G4endl;
783  delete [] dimage;
784 
785  // relative location of the dose distribution image for
786  // the modality image
787  getDoseDistCenterPosition(fCenter, ndose);
788  for(int i = 0; i < 3; i++) iCenter[i] = (int)fCenter[i];
789  if(kLittleEndianOutput) {
790  ofile.write((char *)iCenter, 3*4); // 3*sizeof(int)
791  } else {
792  for(int j = 0; j < 3; j++)
793  invertByteOrder((char *)&iCenter[j], itmp[j]);
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 << G4endl;
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)*2);
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(short));
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] = (int)fCenter[i];
878  if(kLittleEndianOutput) {
879  ofile.write((char *)iCenter, 3*sizeof(int));
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 = 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].getNumberOfSteps();
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 << G4endl;
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; isteps++) {
931  kTracks[nt].getStep(stepPoints[0], stepPoints[1], stepPoints[2],
932  stepPoints[3], stepPoints[4], stepPoints[5],
933  isteps);
934 
935  if(kLittleEndianOutput) {
936  ofile.write((char *)stepPoints, sizeof(float)*6);
937  } else {
938  for(int j = 0; j < 6; j++)
939  invertByteOrder((char *)&stepPoints[j], ftmp[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 = 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].getNumberOfEdges();
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 : " << nedges << G4endl;
973  }
974 
975  // edges
976  float edgePoints[6];
977  for(int ne = 0; ne < nedges; ne++) {
978  kDetectors[nd].getEdge(edgePoints[0], edgePoints[1], edgePoints[2],
979  edgePoints[3], edgePoints[4], edgePoints[5],
980  ne);
981 
982  if(kLittleEndianOutput) {
983  ofile.write((char *)edgePoints, sizeof(float)*6);
984  } else {
985  for(int j = 0; j < 6; j++)
986  invertByteOrder((char *)&edgePoints[j], ftmp[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].getName();
1014  dname.resize(80);
1015  ofile.write((char *)dname.c_str(), 80);
1016  if(DEBUG || kVerbose > 0) {
1017  G4cout << " detector name : " << dname << G4endl;
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: " << kFileName << G4endl;
1029 
1030  return true;
1031 }
1032 bool G4GMocrenIO::storeData3() {
1033 
1034  if(kVerbose > 0) G4cout << ">>>>>>> store data (ver.3) <<<<<<<" << G4endl;
1035  if(kVerbose > 0) G4cout << " " << kFileName << G4endl;
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::binary);
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, sizeof(char));
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 modality image data
1079  ofile.write((char *)&kPointerToModalityData, 4);
1080 
1081  // # of dose distributions
1082  //int nDoseDist = (int)pointerToDoseDistData.size();
1083  int nDoseDist = getNumDoseDist();
1084  ofile.write((char *)&nDoseDist, 4);
1085 
1086  // offset from file starting point to the dose image data
1087  for(int i = 0; i < nDoseDist; i++) {
1088  ofile.write((char *)&kPointerToDoseDistData[i], 4);
1089  }
1090 
1091  // offset from file starting point to the ROI image data
1092  ofile.write((char *)&kPointerToROIData, 4);
1093 
1094  // offset from file starting point to the track data
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., scale : "
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 << "Modality image : ";
1145  for(int i = 0; i < size[2]; i++) {
1146  short * image = kModality.getImage(i);
1147  ofile.write((char *)image, psize*sizeof(short));
1148 
1149  if(DEBUG || kVerbose > 0) G4cout << "[" << i << "]" << image[(size_t)(psize*0.55)] << ", ";
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] =kModalityImageDensityMap[i];
1157  ofile.write((char *)pdmap, msize*sizeof(float));
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; ndose++) {
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 << "] image size : ("
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); // sizeof(shorft)*2
1188 
1189  // dose distribution unit
1190  ofile.write((char *)kDoseUnit.c_str(), 12);
1191  if(DEBUG || kVerbose > 0) {
1192  G4cout << "Dose dist. unit : " << kDoseUnit << G4endl;
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); //sizeof(short)
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 << G4endl;
1224  delete [] dimage;
1225 
1226  // relative location of the dose distribution image for
1227  // the modality image
1228  getDoseDistCenterPosition(fCenter, ndose);
1229  for(int i = 0; i < 3; i++) iCenter[i] = (int)fCenter[i];
1230  ofile.write((char *)iCenter, 3*4); // 3*sizeof(int)
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)*2);
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*sizeof(short));
1273 
1274  }
1275 
1276  // ROI relative location
1277  kRoi[0].getCenterPosition(fCenter);
1278  for(int i = 0; i < 3; i++) iCenter[i] = (int)fCenter[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 = 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() >= G4VisManager::errors)
1305  G4cout << "# of track color information must be the same as # of tracks."
1306  << G4endl;
1307  unsigned char white[3] = {255,255,255}; // default color
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 data (ver.2) <<<<<<<" << G4endl;
1334  if(kVerbose > 0) G4cout << " " << kFileName << G4endl;
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::binary);
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, sizeof(char));
1354 
1355  /*
1356  // event number
1357  ofile.write((char *)&numberOfEvents, sizeof(int));
1358  float imageSpacing[3];
1359  imageSpacing[0] = modalityImageVoxelSpacing[0];
1360  imageSpacing[1] = modalityImageVoxelSpacing[1];
1361  imageSpacing[2] = modalityImageVoxelSpacing[2];
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 modality image data
1378  ofile.write((char *)&kPointerToModalityData, 4);
1379 
1380  // offset from file starting point to the dose image data
1381  ofile.write((char *)&kPointerToDoseDistData[0], 4);
1382 
1383  // offset from file starting point to the ROI image data
1384  ofile.write((char *)&kPointerToROIData, 4);
1385 
1386  // offset from file starting point to the track data
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., scale : "
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 << "Modality image : ";
1435  for(int i = 0; i < size[2]; i++) {
1436  short * image =kModality.getImage(i);
1437  ofile.write((char *)image, psize*sizeof(short));
1438 
1439  if(DEBUG || kVerbose > 0) G4cout << "[" << i << "]" << image[(size_t)(psize*0.55)] << ", ";
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] =kModalityImageDensityMap[i];
1447  ofile.write((char *)pdmap, msize*sizeof(float));
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)*2);
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., scale : "
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*sizeof(short));
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 << G4endl;
1503  delete [] dimage;
1504 
1505  // relative location of the dose distribution image for
1506  // the modality image
1507  kDose[0].getCenterPosition(fCenter);
1508  for(int i = 0; i < 3; i++) iCenter[i] = (int)fCenter[i];
1509  ofile.write((char *)iCenter, 3*sizeof(int));
1510  if(DEBUG || kVerbose > 0) {
1511  G4cout << "Dose dist. image relative location : ("
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)*2);
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*sizeof(short));
1551 
1552  }
1553 
1554  // ROI relative location
1555  kRoi[0].getCenterPosition(fCenter);
1556  for(int i = 0; i < 3; i++) iCenter[i] = (int)fCenter[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 = 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::ios_base::in|std::ios_base::binary);
1598  if(!ifile) {
1599  if (G4VisManager::GetVerbosity() >= G4VisManager::errors)
1600  G4cout << "Cannot open file: " << kFileName
1601  << " in G4GMocrenIO::retrieveData()." << G4endl;
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.3) <<<<<<<" << G4endl;
1616  G4cout << " " << kFileName << G4endl;
1617  retrieveData3();
1618  } else if (ver == 0x04) {
1619  G4cout << ">>>>>>> retrieve data (ver.4) <<<<<<<" << G4endl;
1620  G4cout << " " << kFileName << G4endl;
1621  retrieveData4();
1622  } else {
1623  if (G4VisManager::GetVerbosity() >= G4VisManager::errors) {
1624  G4cout << "Error -- invalid file version : " << (int)ver
1625  << G4endl;
1626  G4cout << " " << kFileName << G4endl;
1627  }
1628  G4Exception("G4GMocrenIO::retrieveDadta()",
1629  "gMocren2001", FatalException,
1630  "Error.");
1631 
1632  }
1633  } else if(std::strncmp(verid, "GRAPE", 5) == 0) {
1634  G4cout << ">>>>>>> retrieve data (ver.2) <<<<<<<" << G4endl;
1635  G4cout << " " << kFileName << G4endl;
1636  retrieveData2();
1637  } else {
1638  if (G4VisManager::GetVerbosity() >= G4VisManager::errors)
1639  G4cout << kFileName << " was not gdd file." << G4endl;
1640  return false;
1641  }
1642 
1643  return true;
1644 }
1645 
1646 bool G4GMocrenIO::retrieveData(char * _filename) {
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::ios_base::in|std::ios_base::binary);
1658  if(!ifile) {
1659  if (G4VisManager::GetVerbosity() >= G4VisManager::errors)
1660  G4cout << "Cannot open file: " << kFileName
1661  << " in G4GMocrenIO::retrieveData3()." << G4endl;
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 version : " << kVersion << G4endl;
1679 
1680  // endian
1681  ifile.read((char *)&kLittleEndianInput, sizeof(char));
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 modality image data
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 dose image data
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 ROI image data
1738  ifile.read((char *)ctmp, 4);
1739  convertEndian(ctmp, kPointerToROIData);
1740 
1741  // offset from file starting point to the track data
1742  ifile.read((char *)ctmp, 4);
1743  convertEndian(ctmp, kPointerToTrackData);
1744 
1745  // offset from file starting point to the detector data
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 && kPointerToDoseDistData.size() == 0 &&
1763  kPointerToROIData == 0 && kPointerToTrackData == 0) {
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, modalityImageVoxelSpacing[0]);
1802  convertEndian(ctmp+sizeof(float), modalityImageVoxelSpacing[1]);
1803  convertEndian(ctmp+2*sizeof(float), modalityImageVoxelSpacing[2]);
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., scale : "
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 << "Modality image (" << psize << "): ";
1835  char * cimage = new char[psize*sizeof(short)];
1836  for(int i = 0; i < size[2]; i++) {
1837  ifile.read((char *)cimage, psize*sizeof(short));
1838  short * mimage = new short[psize];
1839  for(int j = 0; j < psize; j++) {
1840  convertEndian(cimage+j*sizeof(short), mimage[j]);
1841  }
1842  kModality.addImage(mimage);
1843 
1844  if(DEBUG || kVerbose > 0) G4cout << "[" << i << "]" << mimage[(size_t)(psize*0.55)] << ", ";
1845  }
1846  if(DEBUG || kVerbose > 0) G4cout << G4endl;
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 << "msize: " << msize << G4endl;
1852  char * pdmap = new char[msize*sizeof(float)];
1853  ifile.read((char *)pdmap, msize*sizeof(float));
1854  float ftmp;
1855  for(int i = 0; i < (int)msize; i++) {
1856  convertEndian(pdmap+i*sizeof(float), ftmp);
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.size() - 10; i <kModalityImageDensityMap.size(); i++)
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 : " << kDoseUnit << G4endl;
1908  }
1909 
1910  // dose distribution scaling
1911  ifile.read((char *)ctmp, 4); // sizeof(float)
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] = minmax[i]*dscale;
1917  kDose[ndose].setMinMax(dminmax);
1918 
1919  if(DEBUG || kVerbose > 0) {
1920  G4cout << "Dose dist. image min., max., scale : "
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 dist. (" << dsize << "): ";
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(short));
1933  double * dimage = new double[dsize];
1934  for(int xy = 0; xy < dsize; xy++) {
1935  convertEndian(di+xy*sizeof(short), shimage[xy]);
1936  dimage[xy] = shimage[xy]*dscale;
1937  }
1938  kDose[ndose].addImage(dimage);
1939 
1940  if(DEBUG || kVerbose > 0) G4cout << "[" << z << "]" << dimage[(size_t)(dsize*0.55)] << ", ";
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 << G4endl;
1953 
1954  ifile.read((char *)ctmp, 3*4); // 3*sizeof(int)
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] = (float)iCenter[i];
1959  kDose[ndose].setCenterPosition(fCenter);
1960 
1961  if(DEBUG || kVerbose > 0) {
1962  G4cout << "Dose dist. image relative location : ("
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 : " << dosename << G4endl;
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[1]);
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(short));
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[2]);
2035  for(int i = 0; i < 3; i++) fCenter[i] = iCenter[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 << G4endl;
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), steppoint[k]);
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 << steps[0][j] << " ";
2091  int nstp = steps.size();
2092  G4cout << "<-> ";
2093  for(int j = 3; j < 6; j++) G4cout << steps[nstp-1][j] << " ";
2094  G4cout << " rgb( ";
2095  for(int j = 0; j < 3; j++) G4cout << (int)rgb[j] << " ";
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 : " << nedges << G4endl;
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(&cftmp[sizeof(float)*j], edgePoints[j]);
2136  detector.push_back(edgePoints);
2137 
2138  }
2139 
2140  if(DEBUG || kVerbose > 0) {
2141  G4cout << " first edge : (" << detector[0][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 << G4endl;
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 * _filename) {
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::ios_base::in|std::ios_base::binary);
2191  if(!ifile) {
2192  if (G4VisManager::GetVerbosity() >= G4VisManager::errors)
2193  G4cout << "Cannot open file: " << kFileName
2194  << " in G4GMocrenIO::retrieveData3()." << G4endl;
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 version : " << kVersion << G4endl;
2212 
2213  // endian
2214  ifile.read((char *)&kLittleEndianInput, sizeof(char));
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 modality image data
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 dose image data
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 ROI image data
2270  ifile.read((char *)ctmp, 4);
2271  convertEndian(ctmp, kPointerToROIData);
2272 
2273  // offset from file starting point to the track data
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 && kPointerToDoseDistData.size() == 0 &&
2286  kPointerToROIData == 0 && kPointerToTrackData == 0) {
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, modalityImageVoxelSpacing[0]);
2325  convertEndian(ctmp+sizeof(float), modalityImageVoxelSpacing[1]);
2326  convertEndian(ctmp+2*sizeof(float), modalityImageVoxelSpacing[2]);
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., scale : "
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 << "Modality image (" << psize << "): ";
2357  char * cimage = new char[psize*sizeof(short)];
2358  for(int i = 0; i < size[2]; i++) {
2359  ifile.read((char *)cimage, psize*sizeof(short));
2360  short * mimage = new short[psize];
2361  for(int j = 0; j < psize; j++) {
2362  convertEndian(cimage+j*sizeof(short), mimage[j]);
2363  }
2364  kModality.addImage(mimage);
2365 
2366  if(DEBUG || kVerbose > 0) G4cout << "[" << i << "]" << mimage[(size_t)(psize*0.55)] << ", ";
2367  }
2368  if(DEBUG || kVerbose > 0) G4cout << G4endl;
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 << "msize: " << msize << G4endl;
2374  char * pdmap = new char[msize*sizeof(float)];
2375  ifile.read((char *)pdmap, msize*sizeof(float));
2376  float ftmp;
2377  for(int i = 0; i < (int)msize; i++) {
2378  convertEndian(pdmap+i*sizeof(float), ftmp);
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.size() - 10; i <kModalityImageDensityMap.size(); i++)
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 : " << kDoseUnit << G4endl;
2428  }
2429 
2430  // dose distribution scaling
2431  ifile.read((char *)ctmp, 4); // sizeof(float)
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] = minmax[i]*dscale;
2437  kDose[ndose].setMinMax(dminmax);
2438 
2439  if(DEBUG || kVerbose > 0) {
2440  G4cout << "Dose dist. image min., max., scale : "
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 dist. (" << dsize << "): ";
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(short));
2453  double * dimage = new double[dsize];
2454  for(int xy = 0; xy < dsize; xy++) {
2455  convertEndian(di+xy*sizeof(short), shimage[xy]);
2456  dimage[xy] = shimage[xy]*dscale;
2457  }
2458  kDose[ndose].addImage(dimage);
2459 
2460  if(DEBUG || kVerbose > 0) G4cout << "[" << z << "]" << dimage[(size_t)(dsize*0.55)] << ", ";
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 << G4endl;
2473 
2474  ifile.read((char *)ctmp, 3*4); // 3*sizeof(int)
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] = (float)iCenter[i];
2479  kDose[ndose].setCenterPosition(fCenter);
2480 
2481  if(DEBUG || kVerbose > 0) {
2482  G4cout << "Dose dist. image relative location : ("
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[1]);
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(short));
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[2]);
2546  for(int i = 0; i < 3; i++) fCenter[i] = iCenter[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 << G4endl;
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)*3);
2576  if(DEBUG || kVerbose > 0) if(i < 10) G4cout << i << ": " ;
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 << tp[j] << ", ";
2580  }
2581 
2582  ifile.read((char *)ctmp, sizeof(float)*3);
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 << tp[j+3] << ", ";
2586  }
2587  addTrack(tp);
2588  if(DEBUG || kVerbose > 0) if(i < 10) G4cout << G4endl;
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[3];
2600  ifile.read((char *)rgb, 3);
2601  addTrackColor(rgb);
2602 
2603  // v4
2604  for(int j = 0; j < 3; j++) trkcolorv4[j] = rgb[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 * _filename) {
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::ios_base::in|std::ios_base::binary);
2629  if(!ifile) {
2630  if (G4VisManager::GetVerbosity() >= G4VisManager::errors)
2631  G4cout << "Cannot open file: " << kFileName
2632  << " in G4GMocrenIO::retrieveData2()." << G4endl;
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 version : " << kVersion << G4endl;
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, sizeof(char));
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 modality image data
2685  ifile.read((char *)ctmp, 4);
2686  convertEndian(ctmp, kPointerToModalityData);
2687 
2688  // offset from file starting point to the dose image data
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 ROI image data
2695  ifile.read((char *)ctmp, 4);
2696  convertEndian(ctmp, kPointerToROIData);
2697 
2698  // offset from file starting point to the track data
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 && kPointerToDoseDistData.size() == 0 &&
2710  kPointerToROIData == 0 && kPointerToTrackData == 0) {
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, modalityImageVoxelSpacing[0]);
2749  convertEndian(ctmp+sizeof(float), modalityImageVoxelSpacing[1]);
2750  convertEndian(ctmp+2*sizeof(float), modalityImageVoxelSpacing[2]);
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., scale : "
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 << "Modality image (" << psize << "): ";
2775  char * cimage = new char[psize*sizeof(short)];
2776  for(int i = 0; i < size[2]; i++) {
2777  ifile.read((char *)cimage, psize*sizeof(short));
2778  short * mimage = new short[psize];
2779  for(int j = 0; j < psize; j++) {
2780  convertEndian(cimage+j*sizeof(short), mimage[j]);
2781  }
2782  kModality.addImage(mimage);
2783 
2784  if(DEBUG || kVerbose > 0) G4cout << "[" << i << "]" << mimage[(size_t)(psize*0.55)] << ", ";
2785  }
2786  if(DEBUG || kVerbose > 0) G4cout << G4endl;
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 << "msize: " << msize << G4endl;
2792  char * pdmap = new char[msize*sizeof(float)];
2793  ifile.read((char *)pdmap, msize*sizeof(float));
2794  float ftmp;
2795  for(int i = 0; i < (int)msize; i++) {
2796  convertEndian(pdmap+i*sizeof(float), ftmp);
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.size() - 10; i <kModalityImageDensityMap.size(); i++)
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] = minmax[i]*dscale;
2845  kDose[0].setMinMax(dminmax);
2846 
2847  if(DEBUG || kVerbose > 0) {
2848  G4cout << "Dose dist. image min., max., scale : "
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 dist. (" << dsize << "): ";
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(short));
2861  double * dimage = new double[dsize];
2862  for(int xy = 0; xy < dsize; xy++) {
2863  convertEndian(di+xy*sizeof(short), shimage[xy]);
2864  dimage[xy] = shimage[xy]*dscale;
2865  }
2866  kDose[0].addImage(dimage);
2867 
2868  if(DEBUG || kVerbose > 0) G4cout << "[" << z << "]" << dimage[(size_t)(dsize*0.55)] << ", ";
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 << G4endl;
2881 
2882  /* ver 1
2883  float doseDist;
2884  int dosePid;
2885  double * doseData = new double[numDoseImageVoxels];
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; j++) {
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 double[numDoseImageVoxels];
2898  for(int i = 0; i < numDoseImageVoxels; i++) {
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 images
2907  ifile.read(ctmp, 3*sizeof(float));
2908  convertEndian(ctmp, relativeLocation[0]);
2909  convertEndian(ctmp+sizeof(float), relativeLocation[1]);
2910  convertEndian(ctmp+2*sizeof(float), relativeLocation[2]);
2911  */
2912 
2913  // relative location of the dose distribution image for
2914  // the modality image
2915  //ofile.write((char *)relativeLocation, 3*sizeof(float));
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[2]);
2920  for(int i = 0; i < 3; i++) fCenter[i] = (float)iCenter[i];
2921  kDose[0].setCenterPosition(fCenter);
2922 
2923  if(DEBUG || kVerbose > 0) {
2924  G4cout << "Dose dist. image relative location : ("
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[1]);
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(short));
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[2]);
2988  for(int i = 0; i < 3; i++) fCenter[i] = iCenter[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 << G4endl;
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)*3);
3019  if(DEBUG || kVerbose > 0) if(i < 10) G4cout << i << ": " ;
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 << tp[j] << ", ";
3023  }
3024 
3025  ifile.read((char *)ctmp, sizeof(float)*3);
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 << tp[j+3] << ", ";
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) G4cout << G4endl;
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), trackDisplacement[2]); // exchanged with [1]
3050  convertEndian(ctmp+2*sizeof(float), trackDisplacement[1]);
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 * _filename) {
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", "Apr",
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() {return kVersion;}
3121 void G4GMocrenIO::setVersion(std::string & _version) {kVersion = _version;}
3122 
3123 // set endians of input/output data
3124 void G4GMocrenIO::setLittleEndianInput(bool _little) {kLittleEndianInput = _little;}
3125 void G4GMocrenIO::setLittleEndianOutput(bool _little) {kLittleEndianOutput = _little;}
3126 
3127 // voxel spacing
3128 void G4GMocrenIO::setVoxelSpacing(float _spacing[3]) {
3129  for(int i = 0; i < 3; i++) kVoxelSpacing[i] = _spacing[i];
3130 }
3131 void G4GMocrenIO::getVoxelSpacing(float _spacing[3]) {
3132  for(int i = 0; i < 3; i++) _spacing[i] = kVoxelSpacing[i];
3133 }
3134 
3135 // get & set number of events
3136 int & G4GMocrenIO::getNumberOfEvents() {
3137  return kNumberOfEvents;
3138 }
3139 void G4GMocrenIO::setNumberOfEvents(int & _numberOfEvents) {
3140  kNumberOfEvents = _numberOfEvents;
3141 }
3142 void G4GMocrenIO::addOneEvent() {
3143  kNumberOfEvents++;
3144 }
3145 
3146 // set/get pointer the modality image data
3147 void G4GMocrenIO::setPointerToModalityData(unsigned int & _pointer) {
3148  kPointerToModalityData = _pointer;
3149 }
3150 unsigned int G4GMocrenIO::getPointerToModalityData() {
3151  return kPointerToModalityData;
3152 }
3153 // set/get pointer the dose distribution image data
3154 void G4GMocrenIO::addPointerToDoseDistData(unsigned int & _pointer) {
3155  kPointerToDoseDistData.push_back(_pointer);
3156 }
3157 unsigned int G4GMocrenIO::getPointerToDoseDistData(int _elem) {
3158  if(kPointerToDoseDistData.size() == 0 ||
3159  kPointerToDoseDistData.size() < (size_t)_elem)
3160  return 0;
3161  else
3162  return kPointerToDoseDistData[_elem];
3163 }
3164 
3165 // set/get pointer the ROI image data
3166 void G4GMocrenIO::setPointerToROIData(unsigned int & _pointer) {
3167  kPointerToROIData = _pointer;
3168 }
3169 unsigned int G4GMocrenIO::getPointerToROIData() {
3170  return kPointerToROIData;
3171 }
3172 // set/get pointer the track data
3173 void G4GMocrenIO::setPointerToTrackData(unsigned int & _pointer) {
3174  kPointerToTrackData = _pointer;
3175 }
3176 unsigned int G4GMocrenIO::getPointerToTrackData() {
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 "pointer to the detector data" except for "pointer to the dose dist data"
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]*2 + 80;
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 = 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].getNumberOfSteps();
3233  pointer += 4 + 3 + nsteps*(4*6); // # of steps + color + steps(float*6)
3234  }
3235  } else {
3236  unsigned int pointer0 = 0;
3237  setPointerToTrackData(pointer0);
3238  }
3239  if(kVerbose > 0) G4cout << " pointer to the track data :"
3240  << kPointerToTrackData << G4endl;
3241 
3242  // pointer to detector data
3243  int ndet = kDetectors.size();
3244  if(ndet != 0) {
3245  kPointerToDetectorData = pointer;
3246  } else {
3247  kPointerToDetectorData = 0;
3248  }
3249  if(kVerbose > 0) G4cout << " pointer to the detector data :"
3250  << kPointerToDetectorData << G4endl;
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 "pointer to the track data" except for "pointer to the dose dist data"
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]*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 _size[3]) {
3363 
3364  kModality.getSize(_size);
3365 }
3366 void G4GMocrenIO::setModalityImageSize(int _size[3]) {
3367 
3368  kModality.setSize(_size);
3369 }
3370 
3371 // get & set the modality image size
3372 void G4GMocrenIO::setModalityImageScale(double & _scale) {
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 * _image) {
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(std::vector<float> & _map) {
3396  kModalityImageDensityMap = _map;
3397 }
3398 std::vector<float> & G4GMocrenIO::getModalityImageDensityMap() {
3399  return kModalityImageDensityMap;
3400 }
3401 // set the modality image min./max.
3402 void G4GMocrenIO::setModalityImageMinMax(short _minmax[2]) {
3403 
3404  kModality.setMinMax(_minmax);
3405 }
3406 // get the modality image min./max.
3407 void G4GMocrenIO::getModalityImageMinMax(short _minmax[2]) {
3408 
3409  short minmax[2];
3410  kModality.getMinMax(minmax);
3411  for(int i = 0; i < 2; i++) _minmax[i] = minmax[i];
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 center
3426 void G4GMocrenIO::setModalityCenterPosition(float _center[3]) {
3427 
3428  kModality.setCenterPosition(_center);
3429 }
3430 void G4GMocrenIO::getModalityCenterPosition(float _center[3]) {
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::string & _unit) {
3442  kModalityUnit = _unit;
3443 }
3444 //
3445 short G4GMocrenIO::convertDensityToHU(float & _dens) {
3446  short rval = -1024; // default: air
3447  int nmap = (int)kModalityImageDensityMap.size();
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] << G4endl;
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 _num) {
3476  // to avoid a warning in the compile process
3477  if(kDoseUnit.size() > static_cast<size_t>(_num)) return kDoseUnit;
3478 
3479  return kDoseUnit;
3480 }
3481 void G4GMocrenIO::setDoseDistUnit(std::string & _unit, int _num) {
3482  // to avoid a warning in the compile process
3483  if(_unit.size() > static_cast<size_t>(_num)) kDoseUnit = _unit;
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], int _num) {
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], int _num) {
3498 
3499  kDose[_num].setSize(_size);
3500 
3501  //resetDose();
3502 }
3503 
3504 void G4GMocrenIO::setDoseDistMinMax(short _minmax[2], int _num) {
3505 
3506  double minmax[2];
3507  double scale = kDose[_num].getScale();
3508  for(int i = 0; i < 2; i++) minmax[i] = (double)_minmax[i]*scale;
3509  kDose[_num].setMinMax(minmax);
3510 }
3511 void G4GMocrenIO::getDoseDistMinMax(short _minmax[2], int _num) {
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] = (short)(minmax[i]/scale+0.5);
3520  }
3521 }
3522 void G4GMocrenIO::setDoseDistMinMax(double _minmax[2], int _num) {
3523 
3524  kDose[_num].setMinMax(_minmax);
3525 }
3526 void G4GMocrenIO::getDoseDistMinMax(double _minmax[2], int _num) {
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 scale
3535 void G4GMocrenIO::setDoseDistScale(double & _scale, int _num) {
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 * _image, int _num) {
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] = ddata[xy];
3568  if(ddata[xy] > minmax[1]) minmax[1] = ddata[xy];
3569  }
3570  kDose[_num].addImage(ddata);
3571 
3572  // set min./max.
3573  kDose[_num].setMinMax(minmax);
3574 }
3575 void G4GMocrenIO::getShortDoseDist(short * _data, int _z, int _num) {
3576 
3577  if(_data == NULL) {
3578  if (G4VisManager::GetVerbosity() >= G4VisManager::errors)
3579  G4cout << "In G4GMocrenIO::getShortDoseDist(), "
3580  << "first argument is NULL pointer. "
3581  << "The argument must be allocated array."
3582  << G4endl;
3583  G4Exception("G4GMocrenIO::getShortDoseDist()",
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); //there is never negative value
3596  }
3597 }
3598 void G4GMocrenIO::getShortDoseDistMinMax(short _minmax[2], int _num) {
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, int _num) {
3607 
3608  kDose[_num].addImage(_image);
3609 }
3610 double * G4GMocrenIO::getDoseDist(int _z, int _num) {
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 * & _image, int _z, int _num) {
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<double *> & _image, int _num) {
3634 
3635  int size[3];
3636  getDoseDistSize(size, _num);
3637  std::vector<double *> dosedist = kDose[_num].getImage();
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) {doseImageDensityMap = _map;};
3649 // set the dose distribution image displacement
3650 void G4GMocrenIO::setDoseDistCenterPosition(float _center[3], int _num) {
3651 
3652  kDose[_num].setCenterPosition(_center);
3653 }
3654 void G4GMocrenIO::getDoseDistCenterPosition(float _center[3], int _num) {
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 _name, int _num) {
3663 
3664  kDose[_num].setName(_name);
3665 }
3666 std::string G4GMocrenIO::getDoseDistName(int _num) {
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<class GMocrenDataPrimitive<double> > & _dose) {
3676  std::vector<class GMocrenDataPrimitive<double> >::iterator itr;
3677  for(itr = kDose.begin(); itr != kDose.end(); itr++) {
3678  _dose.push_back(*itr);
3679  }
3680 }
3681 // merge two dose distributions
3682 bool G4GMocrenIO::mergeDoseDist(std::vector<class GMocrenDataPrimitive<double> > & _dose) {
3683  if(kDose.size() != _dose.size()) {
3684  if (G4VisManager::GetVerbosity() >= G4VisManager::errors) {
3685  G4cout << "G4GMocrenIO::mergeDoseDist() : Error" << G4endl;
3686  G4cout << " Unable to merge the dose distributions,"<< G4endl;
3687  G4cout << " because of different size of dose maps."<< G4endl;
3688  }
3689  return false;
3690  }
3691 
3692  int num = kDose.size();
3693  std::vector<class GMocrenDataPrimitive<double> >::iterator itr1 = kDose.begin();
3694  std::vector<class GMocrenDataPrimitive<double> >::iterator itr2 = _dose.begin();
3695  for(int i = 0; i < num; i++, itr1++, itr2++) {
3696  if (G4VisManager::GetVerbosity() >= G4VisManager::errors)
3697  if(kVerbose > 0)
3698  G4cout << "merged dose distribution [" << i << "]" << G4endl;
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() >= G4VisManager::errors)
3718  // G4cout << "!!! dose distribution data is empty." << G4endl;
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, int _num) {
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 _num) {
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], int _num) {
3776 
3777  return kRoi[_num].setSize(_size);
3778 }
3779 void G4GMocrenIO::getROISize(int _size[3], int _num) {
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[2], int _num) {
3788 
3789  kRoi[_num].setMinMax(_minmax);
3790 }
3791 void G4GMocrenIO::getROIMinMax(short _minmax[2], int _num) {
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 _center[3], int _num) {
3800 
3801  kRoi[_num].setCenterPosition(_center);
3802 }
3803 void G4GMocrenIO::getROICenterPosition(float _center[3], int _num) {
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() >= G4VisManager::errors)
3824  // G4cout << "!!! ROI data is empty." << G4endl;
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 *> & _tracks) {
3845  kSteps = _tracks;
3846 }
3847 std::vector<float *> & G4GMocrenIO::getTracks() {
3848  return kSteps;
3849 }
3850 void G4GMocrenIO::addTrackColor(unsigned char * _colors) {
3851  kStepColors.push_back(_colors);
3852 }
3853 void G4GMocrenIO::setTrackColors(std::vector<unsigned char *> & _trackColors) {
3854  kStepColors = _trackColors;
3855 }
3856 std::vector<unsigned char *> & G4GMocrenIO::getTrackColors() {
3857  return kStepColors;
3858 }
3859 void G4GMocrenIO::copyTracks(std::vector<float *> & _tracks,
3860  std::vector<unsigned char *> & _colors) {
3861  std::vector<float *>::iterator titr;
3862  for(titr = kSteps.begin(); titr != kSteps.end(); titr++) {
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 != kStepColors.end(); citr++) {
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<float *> & _tracks,
3880  std::vector<unsigned char *> & _colors) {
3881  std::vector<float *>::iterator titr;
3882  for(titr = _tracks.begin(); titr != _tracks.end(); titr++) {
3883  addTrack(*titr);
3884  }
3885 
3886  std::vector<unsigned char *>::iterator citr;
3887  for(citr = _colors.begin(); citr != _colors.end(); citr++) {
3888  addTrackColor(*citr);
3889  }
3890 }
3891 void G4GMocrenIO::addTrack(std::vector<float *> & _steps, unsigned char _color[3]) {
3892 
3893  std::vector<float *>::iterator itr = _steps.begin();
3894  std::vector<struct GMocrenTrack::Step> steps;
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::vector<float *> & _steps,
3910  std::vector<unsigned char *> & _color) {
3911 
3912  if(_num > (int)kTracks.size()) {
3913  if (G4VisManager::GetVerbosity() >= G4VisManager::errors)
3914  G4cout << "ERROR in getTrack() : " << G4endl;
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], stepPoints[1], stepPoints[2],
3928  stepPoints[3], stepPoints[4], stepPoints[5],
3929  isteps);
3930  _steps.push_back(stepPoints);
3931  }
3932 }
3933 
3934 void G4GMocrenIO::translateTracks(std::vector<float> & _translate) {
3935  std::vector<class GMocrenTrack>::iterator itr = kTracks.begin();
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 & _name,
3949  std::vector<float *> & _det,
3950  unsigned char _color[3]) {
3951 
3952  std::vector<float *>::iterator itr = _det.begin();
3953  std::vector<struct GMocrenDetector::Edge> edges;
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::vector<float *> & _edges,
3971  std::vector<unsigned char *> & _color,
3972  std::string & _detName) {
3973 
3974  if(_num > (int)kDetectors.size()) {
3975  if (G4VisManager::GetVerbosity() >= G4VisManager::errors)
3976  G4cout << "ERROR in getDetector() : " << G4endl;
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].getNumberOfEdges();
3991  for(int ne = 0; ne < nedges; ne++) {
3992  float * edgePoints = new float[6];
3993  kDetectors[_num].getEdge(edgePoints[0], edgePoints[1], edgePoints[2],
3994  edgePoints[3], edgePoints[4], edgePoints[5],
3995  ne);
3996  _edges.push_back(edgePoints);
3997  }
3998 }
3999 
4000 void G4GMocrenIO::translateDetector(std::vector<float> & _translate) {
4001  std::vector<class GMocrenDetector>::iterator itr = kDetectors.begin();
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, T & _rval) {
4010 
4011  if((kLittleEndianOutput && !kLittleEndianInput) || // big endian
4012  (!kLittleEndianOutput && kLittleEndianInput)) { // little endian
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, T & _rval) {
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