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