| Geant4 Cross Reference |
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 #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........oooOO0OOooo........oooOO0OOooo......
42 DicomVFileImage::DicomVFileImage()
43 {
44 theFileMgr = DicomFileMgr::GetInstance();
45 }
46
47 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
48 DicomVFileImage::DicomVFileImage(DcmDataset* dset) : DicomVFile(dset)
49 {
50 theFileMgr = DicomFileMgr::GetInstance();
51 }
52
53 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
54 void DicomVFileImage::ReadData()
55 {
56 std::vector<double> dImagePositionPatient = Read1Data(theDataset, DCM_ImagePositionPatient, 3);
57 fLocation = dImagePositionPatient[2];
58 std::vector<double> dSliceThickness = Read1Data(theDataset, DCM_SliceThickness, 1);
59 std::vector<double> dPixelSpacing = Read1Data(theDataset, DCM_PixelSpacing, 2);
60
61 std::vector<double> dRows = Read1Data(theDataset, DCM_Rows, 1);
62 std::vector<double> dColumns = Read1Data(theDataset, DCM_Columns, 1);
63 fNoVoxelsY = dRows[0];
64 fNoVoxelsX = dColumns[0];
65 fNoVoxelsZ = 1;
66
67 fMinX = dImagePositionPatient[0]; // center of upper corner of pixel?
68 fMaxX = dImagePositionPatient[0] + dColumns[0] * dPixelSpacing[0];
69
70 fMinY = dImagePositionPatient[1];
71 fMaxY = dImagePositionPatient[1] + dRows[0] * dPixelSpacing[1];
72
73 fMinZ = dImagePositionPatient[2] - dSliceThickness[0] / 2.;
74 fMaxZ = dImagePositionPatient[2] + dSliceThickness[0] / 2.;
75 fVoxelDimX = dPixelSpacing[0];
76 fVoxelDimY = dPixelSpacing[1];
77 fVoxelDimZ = dSliceThickness[0];
78
79 if (DicomFileMgr::verbose >= debugVerb)
80 G4cout << " DicomVFileImage::ReadData: fNoVoxels " << fNoVoxelsX << " " << fNoVoxelsY << " "
81 << fNoVoxelsZ << G4endl;
82 if (DicomFileMgr::verbose >= debugVerb)
83 G4cout << " DicomVFileImage::ReadData: fMin " << fMinX << " " << fMinY << " " << fMinZ
84 << G4endl;
85 if (DicomFileMgr::verbose >= debugVerb)
86 G4cout << " DicomVFileImage::ReadData: fMax " << fMaxX << " " << fMaxY << " " << fMaxZ
87 << G4endl;
88 if (DicomFileMgr::verbose >= debugVerb)
89 G4cout << " DicomVFileImage::ReadData: fVoxelDim " << fVoxelDimX << " " << fVoxelDimY << " "
90 << fVoxelDimZ << G4endl;
91
92 std::vector<double> dImageOrientationPatient =
93 Read1Data(theDataset, DCM_ImageOrientationPatient, 6);
94 fOrientationRows = G4ThreeVector(dImageOrientationPatient[0], dImageOrientationPatient[1],
95 dImageOrientationPatient[2]);
96 fOrientationColumns = G4ThreeVector(dImageOrientationPatient[3], dImageOrientationPatient[4],
97 dImageOrientationPatient[5]);
98
99 if (fOrientationRows != G4ThreeVector(1, 0, 0) || fOrientationColumns != G4ThreeVector(0, 1, 0)) {
100 G4cerr << " OrientationRows " << fOrientationRows << " OrientationColumns "
101 << fOrientationColumns << G4endl;
102 G4Exception("DicomVFileImage::ReadData", "DFCT0002", JustWarning,
103 "OrientationRows must be (1,0,0) and OrientationColumns (0,1,0), please contact "
104 "GAMOS authors");
105 }
106 fBitAllocated = Read1Data(theDataset, DCM_BitsAllocated, 1)[0];
107 if (DicomFileMgr::verbose >= 4) G4cout << " BIT ALLOCATED " << fBitAllocated << G4endl;
108
109 std::vector<double> dRescaleSlope = Read1Data(theDataset, DCM_RescaleSlope, 1);
110 if (dRescaleSlope.size() == 1) {
111 fRescaleSlope = dRescaleSlope[0];
112 }
113 else {
114 fRescaleSlope = 1;
115 }
116 std::vector<double> dRescaleIntercept = Read1Data(theDataset, DCM_RescaleIntercept, 1);
117 if (dRescaleIntercept.size() == 1) {
118 fRescaleIntercept = dRescaleIntercept[0];
119 }
120 else {
121 fRescaleIntercept = 1;
122 }
123
124 ReadPixelData();
125 }
126
127 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
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_PixelData, element);
135 if (result.bad() || element == NULL) {
136 G4Exception("ReadData", "findAndGetElement(DCM_PixelData, ", FatalException,
137 ("Element PixelData not found: " + G4String(result.text())).c_str());
138 }
139 DcmPixelData* dpix = NULL;
140 dpix = OFstatic_cast(DcmPixelData*, element);
141 // If we have compressed data, we must utilize DcmPixelSequence
142 // in order to access it in raw format, e. g. for decompressing it
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 right representation of the data within DCMTK
148 dpix->getOriginalRepresentationKey(xferSyntax, rep);
149 // Access original data representation and get result within pixel sequence
150 result = dpix->getEncapsulatedRepresentation(xferSyntax, rep, dseq);
151 if (result == EC_Normal) // COMPRESSED DATA
152 {
153 G4Exception("DicomVFileImage::ReadData()", "DFCT004", FatalException,
154 "Compressed pixel data is not supported");
155
156 if (DicomFileMgr::verbose >= debugVerb)
157 G4cout << " DicomVFileImage::ReadData: result == EC_Normal Reading compressed data "
158 << std::endl;
159 DcmPixelItem* pixitem = NULL;
160 // Access first frame (skipping offset table)
161 for (int ii = 1; ii < 2; ii++) {
162 OFCondition cond = dseq->getItem(pixitem, ii);
163 if (!cond.good()) break;
164 G4cout << ii << " PIX LENGTH " << pixitem->getLength() << G4endl;
165 }
166 if (pixitem == NULL) {
167 G4Exception("ReadData", "dseq->getItem()", FatalException,
168 "No DcmPixelItem in DcmPixelSequence");
169 }
170 Uint8* pixData = NULL;
171 // Get the length of this pixel item
172 // (i.e. fragment, i.e. most of the time, the lenght of the frame)
173 Uint32 length = pixitem->getLength();
174 if (length == 0) {
175 G4Exception("ReadData", "pixitem->getLength()", FatalException, "PixelData empty");
176 }
177
178 if (DicomFileMgr::verbose >= debugVerb)
179 G4cout << " DicomVFileImage::ReadData: number of pixels " << length << G4endl;
180 // Finally, get the compressed data for this pixel item
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)).good()) {
187 G4Exception("ReadData", "getUint8Array pixData, ", FatalException,
188 ("PixelData not found: " + G4String(result.text())).c_str());
189 }
190 for (int ir = 0; ir < fNoVoxelsY; ir++) {
191 for (int ic = 0; ic < fNoVoxelsX; ic++) {
192 fHounsfieldV.push_back(pixData[ic + ir * fNoVoxelsX] * fRescaleSlope + fRescaleIntercept);
193 }
194 }
195 }
196 else if (fBitAllocated == 16) { // Case 16 bits :
197 Uint16* pixData = NULL;
198 if (!(element->getUint16Array(pixData)).good()) {
199 G4Exception("ReadData", "getUint16Array pixData, ", FatalException,
200 ("PixelData not found: " + G4String(result.text())).c_str());
201 }
202 for (int ir = 0; ir < fNoVoxelsY; ir++) {
203 for (int ic = 0; ic < fNoVoxelsX; ic++) {
204 fHounsfieldV.push_back(pixData[ic + ir * fNoVoxelsX] * fRescaleSlope + fRescaleIntercept);
205 }
206 }
207 }
208 else if (fBitAllocated == 32) { // Case 32 bits :
209 Uint32* pixData = NULL;
210 if (!(element->getUint32Array(pixData)).good()) {
211 G4Exception("ReadData", "getUint32Array pixData, ", FatalException,
212 ("PixelData not found: " + G4String(result.text())).c_str());
213 }
214 for (int ir = 0; ir < fNoVoxelsY; ir++) {
215 for (int ic = 0; ic < fNoVoxelsX; ic++) {
216 fHounsfieldV.push_back(pixData[ic + ir * fNoVoxelsX] * fRescaleSlope + fRescaleIntercept);
217 }
218 }
219 }
220 }
221 }
222
223 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
224 void DicomVFileImage::operator+=(const DicomVFileImage& rhs)
225 {
226 *this = *this + rhs;
227 }
228
229 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
230 DicomVFileImage DicomVFileImage::operator+(const DicomVFileImage& rhs)
231 {
232 //----- Check that both slices has the same dimensions
233 if (fNoVoxelsX != rhs.GetNoVoxelsX() || fNoVoxelsY != rhs.GetNoVoxelsY()) {
234 G4cerr << "DicomVFileImage error adding two slice headers:\
235 !!! Different number of voxels: "
236 << " X= " << fNoVoxelsX << " =? " << rhs.GetNoVoxelsX() << " Y= " << fNoVoxelsY
237 << " =? " << rhs.GetNoVoxelsY() << " Z= " << fNoVoxelsZ << " =? " << rhs.GetNoVoxelsZ()
238 << G4endl;
239 G4Exception("DicomVFileImage::DicomVFileImage", "", FatalErrorInArgument, "");
240 }
241 //----- Check that both slices has the same extensions
242 if (fMinX != rhs.GetMinX() || fMaxX != rhs.GetMaxX() || fMinY != rhs.GetMinY()
243 || fMaxY != rhs.GetMaxY())
244 {
245 G4cerr << "DicomVFileImage error adding two slice headers:\
246 !!! Different extensions: "
247 << " Xmin= " << fMinX << " =? " << rhs.GetMinX() << " Xmax= " << fMaxX << " =? "
248 << rhs.GetMaxX() << " Ymin= " << fMinY << " =? " << rhs.GetMinY() << " Ymax= " << fMaxY
249 << " =? " << rhs.GetMaxY() << G4endl;
250 G4Exception("DicomVFileImage::operator+", "", FatalErrorInArgument, "");
251 }
252
253 //----- Check that both slices has the same orientations
254 if (fOrientationRows != rhs.GetOrientationRows()
255 || fOrientationColumns != rhs.GetOrientationColumns())
256 {
257 G4cerr << "DicomVFileImage error adding two slice headers: !!!\
258 Slices have different orientations "
259 << " Orientation Rows = " << fOrientationRows << " & " << rhs.GetOrientationRows()
260 << " Orientation Columns " << fOrientationColumns << " & "
261 << rhs.GetOrientationColumns() << G4endl;
262 G4Exception("DicomVFileImage::operator+", "", FatalErrorInArgument, "");
263 }
264
265 //----- Check that the slices are contiguous in Z
266 if (std::fabs(fMinZ - rhs.GetMaxZ()) > G4GeometryTolerance::GetInstance()->GetRadialTolerance()
267 && std::fabs(fMaxZ - rhs.GetMinZ())
268 > G4GeometryTolerance::GetInstance()->GetRadialTolerance())
269 {
270 G4cerr << "DicomVFileImage error adding two slice headers: !!!\
271 Slices are not contiguous in Z "
272 << " Zmin= " << fMinZ << " & " << rhs.GetMinZ() << " Zmax= " << fMaxZ << " & "
273 << rhs.GetMaxZ() << G4endl;
274 G4Exception("DicomVFileImage::operator+", "", JustWarning, "");
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.GetNoVoxelsZ());
284
285 return temp;
286 }
287
288 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
289 void DicomVFileImage::DumpHeaderToTextFile(std::ofstream& fout)
290 {
291 if (DicomFileMgr::verbose >= warningVerb)
292 G4cout << fLocation << " DumpHeaderToTextFile " << fFileName << " " << fHounsfieldV.size()
293 << G4endl;
294
295 G4String fName = fFileName.substr(0, fFileName.length() - 3) + "g4dcm";
296 std::ofstream out(fName.c_str());
297
298 if (DicomFileMgr::verbose >= warningVerb)
299 G4cout << "### DicomVFileImage::Dumping Z Slice header to Text file " << G4endl;
300
301 G4int fCompress = theFileMgr->GetCompression();
302 fout << fNoVoxelsX / fCompress << " " << fNoVoxelsY / fCompress << " " << fNoVoxelsZ << std::endl;
303 fout << fMinX << " " << fMaxX << std::endl;
304 fout << fMinY << " " << fMaxY << std::endl;
305 fout << fMinZ << " " << fMaxZ << std::endl;
306 }
307
308 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
309 void DicomVFileImage::Print(std::ostream& out)
310 {
311 G4int fCompress = theFileMgr->GetCompression();
312 out << "@@ CT Slice " << fLocation << G4endl;
313
314 out << "@ NoVoxels " << fNoVoxelsX / fCompress << " " << fNoVoxelsY / fCompress << " "
315 << fNoVoxelsZ << G4endl;
316 out << "@ DIM X: " << fMinX << " " << fMaxX << " Y: " << fMinY << " " << fMaxY << " Z: " << fMinZ
317 << " " << fMaxZ << G4endl;
318 }
319