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

Last change on this file since 1346 was 1337, checked in by garnier, 15 years ago

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