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