| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 #include "DicomVFileImage.hh" 26 #include "DicomVFileImage.hh"
27 <<
28 #include "DicomFileStructure.hh" 27 #include "DicomFileStructure.hh"
29 #include "DicomROI.hh" 28 #include "DicomROI.hh"
30 #include "dcmtk/dcmdata/dcdeftag.h" << 29
>> 30 #include "G4GeometryTolerance.hh"
>> 31
31 #include "dcmtk/dcmdata/dcfilefo.h" 32 #include "dcmtk/dcmdata/dcfilefo.h"
>> 33 #include "dcmtk/dcmdata/dcdeftag.h"
32 #include "dcmtk/dcmdata/dcpixel.h" 34 #include "dcmtk/dcmdata/dcpixel.h"
33 #include "dcmtk/dcmdata/dcpixseq.h" <<
34 #include "dcmtk/dcmdata/dcpxitem.h" 35 #include "dcmtk/dcmdata/dcpxitem.h"
>> 36 #include "dcmtk/dcmdata/dcpixseq.h"
35 #include "dcmtk/dcmrt/drtimage.h" 37 #include "dcmtk/dcmrt/drtimage.h"
36 38
37 #include "G4GeometryTolerance.hh" <<
38 <<
39 #include <set> 39 #include <set>
40 40
41 //....oooOO0OOooo........oooOO0OOooo........oo 41 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
42 DicomVFileImage::DicomVFileImage() 42 DicomVFileImage::DicomVFileImage()
43 { 43 {
44 theFileMgr = DicomFileMgr::GetInstance(); 44 theFileMgr = DicomFileMgr::GetInstance();
45 } 45 }
46 46
47 //....oooOO0OOooo........oooOO0OOooo........oo 47 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
48 DicomVFileImage::DicomVFileImage(DcmDataset* d 48 DicomVFileImage::DicomVFileImage(DcmDataset* dset) : DicomVFile(dset)
49 { 49 {
50 theFileMgr = DicomFileMgr::GetInstance(); 50 theFileMgr = DicomFileMgr::GetInstance();
51 } 51 }
52 52
>> 53
53 //....oooOO0OOooo........oooOO0OOooo........oo 54 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
54 void DicomVFileImage::ReadData() 55 void DicomVFileImage::ReadData()
55 { 56 {
56 std::vector<double> dImagePositionPatient = << 57 std::vector<double> dImagePositionPatient = Read1Data(theDataset, DCM_ImagePositionPatient,3);
57 fLocation = dImagePositionPatient[2]; 58 fLocation = dImagePositionPatient[2];
58 std::vector<double> dSliceThickness = Read1D 59 std::vector<double> dSliceThickness = Read1Data(theDataset, DCM_SliceThickness, 1);
59 std::vector<double> dPixelSpacing = Read1Dat 60 std::vector<double> dPixelSpacing = Read1Data(theDataset, DCM_PixelSpacing, 2);
60 61
61 std::vector<double> dRows = Read1Data(theDat 62 std::vector<double> dRows = Read1Data(theDataset, DCM_Rows, 1);
62 std::vector<double> dColumns = Read1Data(the 63 std::vector<double> dColumns = Read1Data(theDataset, DCM_Columns, 1);
63 fNoVoxelsY = dRows[0]; << 64 fNoVoxelY = dRows[0];
64 fNoVoxelsX = dColumns[0]; << 65 fNoVoxelX = dColumns[0];
65 fNoVoxelsZ = 1; << 66 fNoVoxelZ = 1;
66 << 67
67 fMinX = dImagePositionPatient[0]; // center << 68 fMinX = dImagePositionPatient[0]; // center of upper corner of pixel?
68 fMaxX = dImagePositionPatient[0] + dColumns[ << 69 fMaxX = dImagePositionPatient[0]+dColumns[0]*dPixelSpacing[0];
69 70
70 fMinY = dImagePositionPatient[1]; 71 fMinY = dImagePositionPatient[1];
71 fMaxY = dImagePositionPatient[1] + dRows[0] << 72 fMaxY = dImagePositionPatient[1]+dRows[0]*dPixelSpacing[1];
72 << 73
73 fMinZ = dImagePositionPatient[2] - dSliceThi << 74 fMinZ = dImagePositionPatient[2]-dSliceThickness[0]/2.;
74 fMaxZ = dImagePositionPatient[2] + dSliceThi << 75 fMaxZ = dImagePositionPatient[2]+dSliceThickness[0]/2.;
75 fVoxelDimX = dPixelSpacing[0]; 76 fVoxelDimX = dPixelSpacing[0];
76 fVoxelDimY = dPixelSpacing[1]; 77 fVoxelDimY = dPixelSpacing[1];
77 fVoxelDimZ = dSliceThickness[0]; 78 fVoxelDimZ = dSliceThickness[0];
78 79
79 if (DicomFileMgr::verbose >= debugVerb) << 80 if( DicomFileMgr::verbose >= debugVerb ) G4cout << " DicomVFileImage::ReadData: fNoVoxel "
80 G4cout << " DicomVFileImage::ReadData: fN << 81 << fNoVoxelX << " " << fNoVoxelY << " " << fNoVoxelZ << G4endl;
81 << fNoVoxelsZ << G4endl; << 82 if( DicomFileMgr::verbose >= debugVerb ) G4cout << " DicomVFileImage::ReadData: fMin "
82 if (DicomFileMgr::verbose >= debugVerb) << 83 << fMinX << " " << fMinY << " " << fMinZ << G4endl;
83 G4cout << " DicomVFileImage::ReadData: fM << 84 if( DicomFileMgr::verbose >= debugVerb ) G4cout << " DicomVFileImage::ReadData: fMax "
84 << G4endl; << 85 << fMaxX << " " << fMaxY << " " << fMaxZ << G4endl;
85 if (DicomFileMgr::verbose >= debugVerb) << 86 if( DicomFileMgr::verbose >= debugVerb ) G4cout << " DicomVFileImage::ReadData: fVoxelDim "
86 G4cout << " DicomVFileImage::ReadData: fM << 87 << fVoxelDimX << " " << fVoxelDimY << " " << fVoxelDimZ << G4endl;
87 << G4endl; << 88
88 if (DicomFileMgr::verbose >= debugVerb) <<
89 G4cout << " DicomVFileImage::ReadData: fV <<
90 << fVoxelDimZ << G4endl; <<
91 <<
92 std::vector<double> dImageOrientationPatient 89 std::vector<double> dImageOrientationPatient =
93 Read1Data(theDataset, DCM_ImageOrientation << 90 Read1Data(theDataset, DCM_ImageOrientationPatient,6);
94 fOrientationRows = G4ThreeVector(dImageOrien << 91 fOrientationRows = G4ThreeVector(dImageOrientationPatient[0],dImageOrientationPatient[1],
95 dImageOrien << 92 dImageOrientationPatient[2]);
96 fOrientationColumns = G4ThreeVector(dImageOr << 93 fOrientationColumns = G4ThreeVector(dImageOrientationPatient[3],dImageOrientationPatient[4],
97 dImageOr << 94 dImageOrientationPatient[5]);
98 << 95
99 if (fOrientationRows != G4ThreeVector(1, 0, << 96 if( fOrientationRows != G4ThreeVector(1,0,0)
100 G4cerr << " OrientationRows " << fOrientat << 97 || fOrientationColumns != G4ThreeVector(0,1,0) ) {
>> 98 G4cerr << " OrientationRows " << fOrientationRows << " OrientationColumns "
101 << fOrientationColumns << G4endl; 99 << fOrientationColumns << G4endl;
102 G4Exception("DicomVFileImage::ReadData", " << 100 G4Exception("DicomVFileImage::ReadData",
103 "OrientationRows must be (1,0, << 101 "DFCT0002",
104 "GAMOS authors"); << 102 JustWarning,
>> 103 "OrientationRows must be (1,0,0) and OrientationColumns (0,1,0), please contact GAMOS authors");
105 } 104 }
106 fBitAllocated = Read1Data(theDataset, DCM_Bi 105 fBitAllocated = Read1Data(theDataset, DCM_BitsAllocated, 1)[0];
107 if (DicomFileMgr::verbose >= 4) G4cout << " << 106 if( DicomFileMgr::verbose >= 4 ) G4cout << " BIT ALLOCATED " << fBitAllocated << G4endl;
108 107
109 std::vector<double> dRescaleSlope = Read1Dat 108 std::vector<double> dRescaleSlope = Read1Data(theDataset, DCM_RescaleSlope, 1);
110 if (dRescaleSlope.size() == 1) { << 109 if( dRescaleSlope.size() == 1 ) {
111 fRescaleSlope = dRescaleSlope[0]; 110 fRescaleSlope = dRescaleSlope[0];
112 } << 111 } else {
113 else { <<
114 fRescaleSlope = 1; 112 fRescaleSlope = 1;
115 } 113 }
116 std::vector<double> dRescaleIntercept = Read 114 std::vector<double> dRescaleIntercept = Read1Data(theDataset, DCM_RescaleIntercept, 1);
117 if (dRescaleIntercept.size() == 1) { << 115 if( dRescaleIntercept.size() == 1 ) {
118 fRescaleIntercept = dRescaleIntercept[0]; 116 fRescaleIntercept = dRescaleIntercept[0];
119 } << 117 } else {
120 else { <<
121 fRescaleIntercept = 1; 118 fRescaleIntercept = 1;
122 } 119 }
123 << 120
124 ReadPixelData(); 121 ReadPixelData();
125 } 122 }
126 123
>> 124
127 //....oooOO0OOooo........oooOO0OOooo........oo 125 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
128 void DicomVFileImage::ReadPixelData() 126 void DicomVFileImage::ReadPixelData()
129 { 127 {
130 // READING THE PIXELS : 128 // READING THE PIXELS :
131 OFCondition result = EC_Normal; 129 OFCondition result = EC_Normal;
132 //---- CHECK IF DATA IS COMPRESSED 130 //---- CHECK IF DATA IS COMPRESSED
133 DcmElement* element = NULL; 131 DcmElement* element = NULL;
134 result = theDataset->findAndGetElement(DCM_P 132 result = theDataset->findAndGetElement(DCM_PixelData, element);
135 if (result.bad() || element == NULL) { 133 if (result.bad() || element == NULL) {
136 G4Exception("ReadData", "findAndGetElement << 134 G4Exception("ReadData",
>> 135 "findAndGetElement(DCM_PixelData, ",
>> 136 FatalException,
137 ("Element PixelData not found: 137 ("Element PixelData not found: " + G4String(result.text())).c_str());
138 } 138 }
139 DcmPixelData* dpix = NULL; << 139 DcmPixelData *dpix = NULL;
140 dpix = OFstatic_cast(DcmPixelData*, element) 140 dpix = OFstatic_cast(DcmPixelData*, element);
141 // If we have compressed data, we must utili 141 // If we have compressed data, we must utilize DcmPixelSequence
142 // in order to access it in raw format, e. 142 // in order to access it in raw format, e. g. for decompressing it
143 // with an external library. 143 // with an external library.
144 DcmPixelSequence* dseq = NULL; << 144 DcmPixelSequence *dseq = NULL;
145 E_TransferSyntax xferSyntax = EXS_Unknown; 145 E_TransferSyntax xferSyntax = EXS_Unknown;
146 const DcmRepresentationParameter* rep = NULL << 146 const DcmRepresentationParameter *rep = NULL;
147 // Find the key that is needed to access the 147 // Find the key that is needed to access the right representation of the data within DCMTK
148 dpix->getOriginalRepresentationKey(xferSynta 148 dpix->getOriginalRepresentationKey(xferSyntax, rep);
149 // Access original data representation and g 149 // Access original data representation and get result within pixel sequence
150 result = dpix->getEncapsulatedRepresentation 150 result = dpix->getEncapsulatedRepresentation(xferSyntax, rep, dseq);
151 if (result == EC_Normal) // COMPRESSED DATA << 151 if ( result == EC_Normal ) // COMPRESSED DATA
152 { 152 {
153 G4Exception("DicomVFileImage::ReadData()", << 153 G4Exception("DicomVFileImage::ReadData()",
>> 154 "DFCT004",
>> 155 FatalException,
154 "Compressed pixel data is not 156 "Compressed pixel data is not supported");
155 << 157
156 if (DicomFileMgr::verbose >= debugVerb) << 158 if( DicomFileMgr::verbose >= debugVerb ) G4cout
157 G4cout << " DicomVFileImage::ReadData: << 159 << " DicomVFileImage::ReadData: result == EC_Normal Reading compressed data " << std::endl;
158 << std::endl; <<
159 DcmPixelItem* pixitem = NULL; 160 DcmPixelItem* pixitem = NULL;
160 // Access first frame (skipping offset tab 161 // Access first frame (skipping offset table)
161 for (int ii = 1; ii < 2; ii++) { << 162 for( int ii = 1; ii < 2; ii++ ) {
162 OFCondition cond = dseq->getItem(pixitem << 163 OFCondition cond = dseq->getItem(pixitem, ii);
163 if (!cond.good()) break; << 164 if( !cond.good()) break;
164 G4cout << ii << " PIX LENGTH " << pixite 165 G4cout << ii << " PIX LENGTH " << pixitem->getLength() << G4endl;
165 } 166 }
166 if (pixitem == NULL) { 167 if (pixitem == NULL) {
167 G4Exception("ReadData", "dseq->getItem() << 168 G4Exception("ReadData",
168 "No DcmPixelItem in DcmPixel << 169 "dseq->getItem()",
>> 170 FatalException,
>> 171 "No DcmPixelItem in DcmPixelSequence");
169 } 172 }
170 Uint8* pixData = NULL; 173 Uint8* pixData = NULL;
171 // Get the length of this pixel item 174 // Get the length of this pixel item
172 // (i.e. fragment, i.e. most of the time, 175 // (i.e. fragment, i.e. most of the time, the lenght of the frame)
173 Uint32 length = pixitem->getLength(); 176 Uint32 length = pixitem->getLength();
174 if (length == 0) { 177 if (length == 0) {
175 G4Exception("ReadData", "pixitem->getLen << 178 G4Exception("ReadData",
>> 179 "pixitem->getLength()",
>> 180 FatalException,
>> 181 "PixelData empty");
176 } 182 }
177 183
178 if (DicomFileMgr::verbose >= debugVerb) << 184 if( DicomFileMgr::verbose >= debugVerb ) G4cout
179 G4cout << " DicomVFileImage::ReadData: << 185 << " DicomVFileImage::ReadData: number of pixels " << length << G4endl;
180 // Finally, get the compressed data for th 186 // Finally, get the compressed data for this pixel item
181 result = pixitem->getUint8Array(pixData); 187 result = pixitem->getUint8Array(pixData);
182 } << 188 } else { // UNCOMPRESSED DATA
183 else { // UNCOMPRESSED DATA << 189 if(fBitAllocated == 8) { // Case 8 bits :
184 if (fBitAllocated == 8) { // Case 8 bits <<
185 Uint8* pixData = NULL; 190 Uint8* pixData = NULL;
186 if (!(element->getUint8Array(pixData)).g << 191 if(! (element->getUint8Array(pixData)).good() ) {
187 G4Exception("ReadData", "getUint8Array << 192 G4Exception("ReadData",
>> 193 "getUint8Array pixData, ",
>> 194 FatalException,
188 ("PixelData not found: " + 195 ("PixelData not found: " + G4String(result.text())).c_str());
189 } 196 }
190 for (int ir = 0; ir < fNoVoxelsY; ir++) << 197 for( int ir = 0; ir < fNoVoxelY; ir++ ) {
191 for (int ic = 0; ic < fNoVoxelsX; ic++ << 198 for( int ic = 0; ic < fNoVoxelX; ic++ ) {
192 fHounsfieldV.push_back(pixData[ic + << 199 fHounsfieldV.push_back(pixData[ic+ir*fNoVoxelX]*fRescaleSlope + fRescaleIntercept);
193 } 200 }
194 } 201 }
195 } << 202 } else if(fBitAllocated == 16) { // Case 16 bits :
196 else if (fBitAllocated == 16) { // Case 1 <<
197 Uint16* pixData = NULL; 203 Uint16* pixData = NULL;
198 if (!(element->getUint16Array(pixData)). << 204 if(! (element->getUint16Array(pixData)).good() ) {
199 G4Exception("ReadData", "getUint16Arra << 205 G4Exception("ReadData",
>> 206 "getUint16Array pixData, ",
>> 207 FatalException,
200 ("PixelData not found: " + 208 ("PixelData not found: " + G4String(result.text())).c_str());
201 } 209 }
202 for (int ir = 0; ir < fNoVoxelsY; ir++) << 210 for( int ir = 0; ir < fNoVoxelY; ir++ ) {
203 for (int ic = 0; ic < fNoVoxelsX; ic++ << 211 for( int ic = 0; ic < fNoVoxelX; ic++ ) {
204 fHounsfieldV.push_back(pixData[ic + << 212 fHounsfieldV.push_back(pixData[ic+ir*fNoVoxelX]*fRescaleSlope + fRescaleIntercept);
205 } 213 }
206 } 214 }
207 } << 215 } else if(fBitAllocated == 32) { // Case 32 bits :
208 else if (fBitAllocated == 32) { // Case 3 <<
209 Uint32* pixData = NULL; 216 Uint32* pixData = NULL;
210 if (!(element->getUint32Array(pixData)). << 217 if(! (element->getUint32Array(pixData)).good() ) {
211 G4Exception("ReadData", "getUint32Arra << 218 G4Exception("ReadData",
>> 219 "getUint32Array pixData, ",
>> 220 FatalException,
212 ("PixelData not found: " + 221 ("PixelData not found: " + G4String(result.text())).c_str());
213 } 222 }
214 for (int ir = 0; ir < fNoVoxelsY; ir++) << 223 for( int ir = 0; ir < fNoVoxelY; ir++ ) {
215 for (int ic = 0; ic < fNoVoxelsX; ic++ << 224 for( int ic = 0; ic < fNoVoxelX; ic++ ) {
216 fHounsfieldV.push_back(pixData[ic + << 225 fHounsfieldV.push_back(pixData[ic+ir*fNoVoxelX]*fRescaleSlope + fRescaleIntercept);
217 } 226 }
218 } 227 }
219 } 228 }
220 } << 229 }
>> 230
221 } 231 }
222 232
223 //....oooOO0OOooo........oooOO0OOooo........oo 233 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
224 void DicomVFileImage::operator+=(const DicomVF << 234 void DicomVFileImage::operator+=( const DicomVFileImage& rhs )
225 { 235 {
226 *this = *this + rhs; 236 *this = *this + rhs;
227 } 237 }
228 238
229 //....oooOO0OOooo........oooOO0OOooo........oo 239 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
230 DicomVFileImage DicomVFileImage::operator+(con << 240 DicomVFileImage DicomVFileImage::operator+( const DicomVFileImage& rhs )
231 { 241 {
232 //----- Check that both slices has the same 242 //----- Check that both slices has the same dimensions
233 if (fNoVoxelsX != rhs.GetNoVoxelsX() || fNoV << 243 if( fNoVoxelX != rhs.GetNoVoxelX()
>> 244 || fNoVoxelY != rhs.GetNoVoxelY() ) {
234 G4cerr << "DicomVFileImage error adding tw 245 G4cerr << "DicomVFileImage error adding two slice headers:\
235 !!! Different number of voxels: " 246 !!! Different number of voxels: "
236 << " X= " << fNoVoxelsX << " =? " << 247 << " X= " << fNoVoxelX << " =? " << rhs.GetNoVoxelX()
237 << " =? " << rhs.GetNoVoxelsY() << << 248 << " Y= " << fNoVoxelY << " =? " << rhs.GetNoVoxelY()
>> 249 << " Z= " << fNoVoxelZ << " =? " << rhs.GetNoVoxelZ()
238 << G4endl; 250 << G4endl;
239 G4Exception("DicomVFileImage::DicomVFileIm << 251 G4Exception("DicomVFileImage::DicomVFileImage",
>> 252 "",FatalErrorInArgument,"");
240 } 253 }
241 //----- Check that both slices has the same 254 //----- Check that both slices has the same extensions
242 if (fMinX != rhs.GetMinX() || fMaxX != rhs.G << 255 if( fMinX != rhs.GetMinX() || fMaxX != rhs.GetMaxX()
243 || fMaxY != rhs.GetMaxY()) << 256 || fMinY != rhs.GetMinY() || fMaxY != rhs.GetMaxY() ) {
244 { <<
245 G4cerr << "DicomVFileImage error adding tw 257 G4cerr << "DicomVFileImage error adding two slice headers:\
246 !!! Different extensions: " 258 !!! Different extensions: "
247 << " Xmin= " << fMinX << " =? " << << 259 << " Xmin= " << fMinX << " =? " << rhs.GetMinX()
248 << rhs.GetMaxX() << " Ymin= " << f << 260 << " Xmax= " << fMaxX << " =? " << rhs.GetMaxX()
249 << " =? " << rhs.GetMaxY() << G4end << 261 << " Ymin= " << fMinY << " =? " << rhs.GetMinY()
250 G4Exception("DicomVFileImage::operator+", << 262 << " Ymax= " << fMaxY << " =? " << rhs.GetMaxY()
>> 263 << G4endl;
>> 264 G4Exception("DicomVFileImage::operator+","",
>> 265 FatalErrorInArgument,"");
251 } 266 }
252 267
253 //----- Check that both slices has the same 268 //----- Check that both slices has the same orientations
254 if (fOrientationRows != rhs.GetOrientationRo << 269 if( fOrientationRows != rhs.GetOrientationRows() ||
255 || fOrientationColumns != rhs.GetOrienta << 270 fOrientationColumns != rhs.GetOrientationColumns() ) {
256 { <<
257 G4cerr << "DicomVFileImage error adding tw 271 G4cerr << "DicomVFileImage error adding two slice headers: !!!\
258 Slices have different orientations " 272 Slices have different orientations "
259 << " Orientation Rows = " << fOrie << 273 << " Orientation Rows = " << fOrientationRows << " & " << rhs.GetOrientationRows()
260 << " Orientation Columns " << fOri << 274 << " Orientation Columns " << fOrientationColumns << " & "
261 << rhs.GetOrientationColumns() << G 275 << rhs.GetOrientationColumns() << G4endl;
262 G4Exception("DicomVFileImage::operator+", << 276 G4Exception("DicomVFileImage::operator+","",
>> 277 FatalErrorInArgument,"");
263 } 278 }
264 << 279
265 //----- Check that the slices are contiguous 280 //----- Check that the slices are contiguous in Z
266 if (std::fabs(fMinZ - rhs.GetMaxZ()) > G4Geo << 281 if( std::fabs( fMinZ - rhs.GetMaxZ() ) >
267 && std::fabs(fMaxZ - rhs.GetMinZ()) << 282 G4GeometryTolerance::GetInstance()->GetRadialTolerance() &&
268 > G4GeometryTolerance::GetInstance( << 283 std::fabs( fMaxZ - rhs.GetMinZ() ) >
269 { << 284 G4GeometryTolerance::GetInstance()->GetRadialTolerance() ){
270 G4cerr << "DicomVFileImage error adding tw 285 G4cerr << "DicomVFileImage error adding two slice headers: !!!\
271 Slices are not contiguous in Z " 286 Slices are not contiguous in Z "
272 << " Zmin= " << fMinZ << " & " << << 287 << " Zmin= " << fMinZ << " & " << rhs.GetMinZ()
273 << rhs.GetMaxZ() << G4endl; << 288 << " Zmax= " << fMaxZ << " & " << rhs.GetMaxZ()
274 G4Exception("DicomVFileImage::operator+", << 289 << G4endl;
>> 290 G4Exception("DicomVFileImage::operator+","",
>> 291 JustWarning,"");
275 } 292 }
276 << 293
277 //----- Build slice header copying first one 294 //----- Build slice header copying first one
278 DicomVFileImage temp(*this); << 295 DicomVFileImage temp( *this );
279 << 296
280 //----- Add data from second slice header 297 //----- Add data from second slice header
281 temp.SetMinZ(std::min(fMinZ, rhs.GetMinZ())) << 298 temp.SetMinZ( std::min( fMinZ, rhs.GetMinZ() ) );
282 temp.SetMaxZ(std::max(fMaxZ, rhs.GetMaxZ())) << 299 temp.SetMaxZ( std::max( fMaxZ, rhs.GetMaxZ() ) );
283 temp.SetNoVoxelsZ(fNoVoxelsZ + rhs.GetNoVoxe << 300 temp.SetNoVoxelZ( fNoVoxelZ + rhs.GetNoVoxelZ() );
284 << 301
285 return temp; 302 return temp;
286 } 303 }
287 304
288 //....oooOO0OOooo........oooOO0OOooo........oo 305 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
289 void DicomVFileImage::DumpHeaderToTextFile(std 306 void DicomVFileImage::DumpHeaderToTextFile(std::ofstream& fout)
290 { 307 {
291 if (DicomFileMgr::verbose >= warningVerb) << 308 if( DicomFileMgr::verbose >= warningVerb ) G4cout << fLocation << " DumpHeaderToTextFile "
292 G4cout << fLocation << " DumpHeaderToTextF << 309 << fFileName << " " << fHounsfieldV.size() << G4endl;
293 << G4endl; <<
294 310
295 G4String fName = fFileName.substr(0, fFileNa << 311 G4String fName = fFileName.substr(0,fFileName.length()-3) + "g4dcm";
296 std::ofstream out(fName.c_str()); 312 std::ofstream out(fName.c_str());
297 313
298 if (DicomFileMgr::verbose >= warningVerb) << 314 if( DicomFileMgr::verbose >= warningVerb ) G4cout
299 G4cout << "### DicomVFileImage::Dumping Z << 315 << "### DicomVFileImage::Dumping Z Slice header to Text file " << G4endl;
300 316
301 G4int fCompress = theFileMgr->GetCompression 317 G4int fCompress = theFileMgr->GetCompression();
302 fout << fNoVoxelsX / fCompress << " " << fNo << 318 fout << fNoVoxelX/fCompress << " " << fNoVoxelY/fCompress << " " << fNoVoxelZ << std::endl;
303 fout << fMinX << " " << fMaxX << std::endl; 319 fout << fMinX << " " << fMaxX << std::endl;
304 fout << fMinY << " " << fMaxY << std::endl; 320 fout << fMinY << " " << fMaxY << std::endl;
305 fout << fMinZ << " " << fMaxZ << std::endl; 321 fout << fMinZ << " " << fMaxZ << std::endl;
>> 322
306 } 323 }
307 324
308 //....oooOO0OOooo........oooOO0OOooo........oo 325 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
309 void DicomVFileImage::Print(std::ostream& out) << 326 void DicomVFileImage::Print(std::ostream& out )
310 { 327 {
311 G4int fCompress = theFileMgr->GetCompression 328 G4int fCompress = theFileMgr->GetCompression();
312 out << "@@ CT Slice " << fLocation << G4endl 329 out << "@@ CT Slice " << fLocation << G4endl;
313 330
314 out << "@ NoVoxels " << fNoVoxelsX / fCompre << 331 out << "@ NoVoxels " << fNoVoxelX/fCompress << " " << fNoVoxelY/fCompress << " "
315 << fNoVoxelsZ << G4endl; << 332 << fNoVoxelZ << G4endl;
316 out << "@ DIM X: " << fMinX << " " << fMaxX << 333 out << "@ DIM X: " << fMinX << " " << fMaxX
317 << " " << fMaxZ << G4endl; << 334 << " Y: " << fMinY << " " << fMaxY
>> 335 << " Z: " << fMinZ << " " << fMaxZ << G4endl;
318 } 336 }
>> 337
319 338