source: trunk/source/visualization/gMocren/src/G4GMocrenIO.cc @ 1347

Last change on this file since 1347 was 1347, checked in by garnier, 13 years ago

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