Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
DicomVFileImage.cc
Go to the documentation of this file.
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 #include "DicomFileStructure.hh"
28 #include "DicomROI.hh"
29 
30 #include "G4GeometryTolerance.hh"
31 
32 #include "dcmtk/dcmdata/dcfilefo.h"
33 #include "dcmtk/dcmdata/dcdeftag.h"
34 #include "dcmtk/dcmdata/dcpixel.h"
35 #include "dcmtk/dcmdata/dcpxitem.h"
36 #include "dcmtk/dcmdata/dcpixseq.h"
37 #include "dcmtk/dcmrt/drtimage.h"
38 
39 #include <set>
40 
41 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
43 {
45 }
46 
47 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
49 {
51 }
52 
53 
54 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
56 {
57  std::vector<double> dImagePositionPatient = Read1Data(theDataset, DCM_ImagePositionPatient,3);
58  fLocation = dImagePositionPatient[2];
59  std::vector<double> dSliceThickness = Read1Data(theDataset, DCM_SliceThickness, 1);
60  std::vector<double> dPixelSpacing = Read1Data(theDataset, DCM_PixelSpacing, 2);
61 
62  std::vector<double> dRows = Read1Data(theDataset, DCM_Rows, 1);
63  std::vector<double> dColumns = Read1Data(theDataset, DCM_Columns, 1);
64  fNoVoxelY = dRows[0];
65  fNoVoxelX = dColumns[0];
66  fNoVoxelZ = 1;
67 
68  fMinX = dImagePositionPatient[0]; // center of upper corner of pixel?
69  fMaxX = dImagePositionPatient[0]+dColumns[0]*dPixelSpacing[0];
70 
71  fMinY = dImagePositionPatient[1];
72  fMaxY = dImagePositionPatient[1]+dRows[0]*dPixelSpacing[1];
73 
74  fMinZ = dImagePositionPatient[2]-dSliceThickness[0]/2.;
75  fMaxZ = dImagePositionPatient[2]+dSliceThickness[0]/2.;
76  fVoxelDimX = dPixelSpacing[0];
77  fVoxelDimY = dPixelSpacing[1];
78  fVoxelDimZ = dSliceThickness[0];
79 
80  if( DicomFileMgr::verbose >= debugVerb ) G4cout << " DicomVFileImage::ReadData: fNoVoxel "
81  << fNoVoxelX << " " << fNoVoxelY << " " << fNoVoxelZ << G4endl;
82  if( DicomFileMgr::verbose >= debugVerb ) G4cout << " DicomVFileImage::ReadData: fMin "
83  << fMinX << " " << fMinY << " " << fMinZ << G4endl;
84  if( DicomFileMgr::verbose >= debugVerb ) G4cout << " DicomVFileImage::ReadData: fMax "
85  << fMaxX << " " << fMaxY << " " << fMaxZ << G4endl;
86  if( DicomFileMgr::verbose >= debugVerb ) G4cout << " DicomVFileImage::ReadData: fVoxelDim "
87  << fVoxelDimX << " " << fVoxelDimY << " " << fVoxelDimZ << G4endl;
88 
89  std::vector<double> dImageOrientationPatient =
90  Read1Data(theDataset, DCM_ImageOrientationPatient,6);
91  fOrientationRows = G4ThreeVector(dImageOrientationPatient[0],dImageOrientationPatient[1],
92  dImageOrientationPatient[2]);
93  fOrientationColumns = G4ThreeVector(dImageOrientationPatient[3],dImageOrientationPatient[4],
94  dImageOrientationPatient[5]);
95 
96  if( fOrientationRows != G4ThreeVector(1,0,0)
97  || fOrientationColumns != G4ThreeVector(0,1,0) ) {
98  G4cerr << " OrientationRows " << fOrientationRows << " OrientationColumns "
100  G4Exception("DicomVFileImage::ReadData",
101  "DFCT0002",
102  JustWarning,
103  "OrientationRows must be (1,0,0) and OrientationColumns (0,1,0), please contact GAMOS authors");
104  }
105  fBitAllocated = Read1Data(theDataset, DCM_BitsAllocated, 1)[0];
106  if( DicomFileMgr::verbose >= 4 ) G4cout << " BIT ALLOCATED " << fBitAllocated << G4endl;
107 
108  std::vector<double> dRescaleSlope = Read1Data(theDataset, DCM_RescaleSlope, 1);
109  if( dRescaleSlope.size() == 1 ) {
110  fRescaleSlope = dRescaleSlope[0];
111  } else {
112  fRescaleSlope = 1;
113  }
114  std::vector<double> dRescaleIntercept = Read1Data(theDataset, DCM_RescaleIntercept, 1);
115  if( dRescaleIntercept.size() == 1 ) {
116  fRescaleIntercept = dRescaleIntercept[0];
117  } else {
118  fRescaleIntercept = 1;
119  }
120 
121  ReadPixelData();
122 }
123 
124 
125 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
126 void DicomVFileImage::ReadPixelData()
127 {
128  // READING THE PIXELS :
129  OFCondition result = EC_Normal;
130  //---- CHECK IF DATA IS COMPRESSED
131  DcmElement* element = NULL;
132  result = theDataset->findAndGetElement(DCM_PixelData, element);
133  if (result.bad() || element == NULL) {
134  G4Exception("ReadData",
135  "findAndGetElement(DCM_PixelData, ",
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()",
154  "DFCT004",
156  "Compressed pixel data is not supported");
157 
159  << " DicomVFileImage::ReadData: result == EC_Normal Reading compressed data " << std::endl;
160  DcmPixelItem* pixitem = NULL;
161  // Access first frame (skipping offset table)
162  for( int ii = 1; ii < 2; ii++ ) {
163  OFCondition cond = dseq->getItem(pixitem, ii);
164  if( !cond.good()) break;
165  G4cout << ii << " PIX LENGTH " << pixitem->getLength() << G4endl;
166  }
167  if (pixitem == NULL) {
168  G4Exception("ReadData",
169  "dseq->getItem()",
171  "No DcmPixelItem in DcmPixelSequence");
172  }
173  Uint8* pixData = NULL;
174  // Get the length of this pixel item
175  // (i.e. fragment, i.e. most of the time, the lenght of the frame)
176  Uint32 length = pixitem->getLength();
177  if (length == 0) {
178  G4Exception("ReadData",
179  "pixitem->getLength()",
181  "PixelData empty");
182  }
183 
185  << " DicomVFileImage::ReadData: number of pixels " << length << G4endl;
186  // Finally, get the compressed data for this pixel item
187  result = pixitem->getUint8Array(pixData);
188  } else { // UNCOMPRESSED DATA
189  if(fBitAllocated == 8) { // Case 8 bits :
190  Uint8* pixData = NULL;
191  if(! (element->getUint8Array(pixData)).good() ) {
192  G4Exception("ReadData",
193  "getUint8Array pixData, ",
195  ("PixelData not found: " + G4String(result.text())).c_str());
196  }
197  for( int ir = 0; ir < fNoVoxelY; ir++ ) {
198  for( int ic = 0; ic < fNoVoxelX; ic++ ) {
199  fHounsfieldV.push_back(pixData[ic+ir*fNoVoxelX]*fRescaleSlope + fRescaleIntercept);
200  }
201  }
202  } else if(fBitAllocated == 16) { // Case 16 bits :
203  Uint16* pixData = NULL;
204  if(! (element->getUint16Array(pixData)).good() ) {
205  G4Exception("ReadData",
206  "getUint16Array pixData, ",
208  ("PixelData not found: " + G4String(result.text())).c_str());
209  }
210  for( int ir = 0; ir < fNoVoxelY; ir++ ) {
211  for( int ic = 0; ic < fNoVoxelX; ic++ ) {
212  fHounsfieldV.push_back(pixData[ic+ir*fNoVoxelX]*fRescaleSlope + fRescaleIntercept);
213  }
214  }
215  } else if(fBitAllocated == 32) { // Case 32 bits :
216  Uint32* pixData = NULL;
217  if(! (element->getUint32Array(pixData)).good() ) {
218  G4Exception("ReadData",
219  "getUint32Array pixData, ",
221  ("PixelData not found: " + G4String(result.text())).c_str());
222  }
223  for( int ir = 0; ir < fNoVoxelY; ir++ ) {
224  for( int ic = 0; ic < fNoVoxelX; ic++ ) {
225  fHounsfieldV.push_back(pixData[ic+ir*fNoVoxelX]*fRescaleSlope + fRescaleIntercept);
226  }
227  }
228  }
229  }
230 
231 }
232 
233 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
235 {
236  *this = *this + rhs;
237 }
238 
239 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
241 {
242  //----- Check that both slices has the same dimensions
243  if( fNoVoxelX != rhs.GetNoVoxelX()
244  || fNoVoxelY != rhs.GetNoVoxelY() ) {
245  G4cerr << "DicomVFileImage error adding two slice headers:\
246  !!! Different number of voxels: "
247  << " X= " << fNoVoxelX << " =? " << rhs.GetNoVoxelX()
248  << " Y= " << fNoVoxelY << " =? " << rhs.GetNoVoxelY()
249  << " Z= " << fNoVoxelZ << " =? " << rhs.GetNoVoxelZ()
250  << G4endl;
251  G4Exception("DicomVFileImage::DicomVFileImage",
252  "",FatalErrorInArgument,"");
253  }
254  //----- Check that both slices has the same extensions
255  if( fMinX != rhs.GetMinX() || fMaxX != rhs.GetMaxX()
256  || fMinY != rhs.GetMinY() || fMaxY != rhs.GetMaxY() ) {
257  G4cerr << "DicomVFileImage error adding two slice headers:\
258  !!! Different extensions: "
259  << " Xmin= " << fMinX << " =? " << rhs.GetMinX()
260  << " Xmax= " << fMaxX << " =? " << rhs.GetMaxX()
261  << " Ymin= " << fMinY << " =? " << rhs.GetMinY()
262  << " Ymax= " << fMaxY << " =? " << rhs.GetMaxY()
263  << G4endl;
264  G4Exception("DicomVFileImage::operator+","",
266  }
267 
268  //----- Check that both slices has the same orientations
269  if( fOrientationRows != rhs.GetOrientationRows() ||
271  G4cerr << "DicomVFileImage error adding two slice headers: !!!\
272  Slices have different orientations "
273  << " Orientation Rows = " << fOrientationRows << " & " << rhs.GetOrientationRows()
274  << " Orientation Columns " << fOrientationColumns << " & "
275  << rhs.GetOrientationColumns() << G4endl;
276  G4Exception("DicomVFileImage::operator+","",
278  }
279 
280  //----- Check that the slices are contiguous in Z
281  if( std::fabs( fMinZ - rhs.GetMaxZ() ) >
283  std::fabs( fMaxZ - rhs.GetMinZ() ) >
285  G4cerr << "DicomVFileImage error adding two slice headers: !!!\
286  Slices are not contiguous in Z "
287  << " Zmin= " << fMinZ << " & " << rhs.GetMinZ()
288  << " Zmax= " << fMaxZ << " & " << rhs.GetMaxZ()
289  << G4endl;
290  G4Exception("DicomVFileImage::operator+","",
291  JustWarning,"");
292  }
293 
294  //----- Build slice header copying first one
295  DicomVFileImage temp( *this );
296 
297  //----- Add data from second slice header
298  temp.SetMinZ( std::min( fMinZ, rhs.GetMinZ() ) );
299  temp.SetMaxZ( std::max( fMaxZ, rhs.GetMaxZ() ) );
300  temp.SetNoVoxelZ( fNoVoxelZ + rhs.GetNoVoxelZ() );
301 
302  return temp;
303 }
304 
305 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
306 void DicomVFileImage::DumpHeaderToTextFile(std::ofstream& fout)
307 {
308  if( DicomFileMgr::verbose >= warningVerb ) G4cout << fLocation << " DumpHeaderToTextFile "
309  << fFileName << " " << fHounsfieldV.size() << G4endl;
310 
311  G4String fName = fFileName.substr(0,fFileName.length()-3) + "g4dcm";
312  std::ofstream out(fName.c_str());
313 
315  << "### DicomVFileImage::Dumping Z Slice header to Text file " << G4endl;
316 
317  G4int fCompress = theFileMgr->GetCompression();
318  fout << fNoVoxelX/fCompress << " " << fNoVoxelY/fCompress << " " << fNoVoxelZ << std::endl;
319  fout << fMinX << " " << fMaxX << std::endl;
320  fout << fMinY << " " << fMaxY << std::endl;
321  fout << fMinZ << " " << fMaxZ << std::endl;
322 
323 }
324 
325 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
326 void DicomVFileImage::Print(std::ostream& out )
327 {
328  G4int fCompress = theFileMgr->GetCompression();
329  out << "@@ CT Slice " << fLocation << G4endl;
330 
331  out << "@ NoVoxels " << fNoVoxelX/fCompress << " " << fNoVoxelY/fCompress << " "
332  << fNoVoxelZ << G4endl;
333  out << "@ DIM X: " << fMinX << " " << fMaxX
334  << " Y: " << fMinY << " " << fMaxY
335  << " Z: " << fMinZ << " " << fMaxZ << G4endl;
336 }
337 
G4double G4ParticleHPJENDLHEData::G4double result
G4double fRescaleIntercept
CLHEP::Hep3Vector G4ThreeVector
G4String fName
Definition: G4AttUtils.hh:55
void SetMinZ(const G4double &val)
G4double GetMinX() const
G4int GetCompression() const
Definition: DicomFileMgr.hh:91
G4double fRescaleSlope
G4double GetMinZ() const
virtual std::vector< G4double > Read1Data(DcmDataset *dset, DcmTagKey tagKey, G4int nData)
Definition: DicomVFile.cc:40
G4ThreeVector GetOrientationRows() const
G4int GetNoVoxelY() const
int G4int
Definition: G4Types.hh:78
G4double fBitAllocated
void operator+=(const DicomVFileImage &rhs)
G4String fFileName
Definition: DicomVFile.hh:55
G4int GetNoVoxelZ() const
void SetMaxZ(const G4double &val)
G4GLOB_DLL std::ostream G4cout
DicomVFileImage operator+(const DicomVFileImage &rhs)
G4double GetRadialTolerance() const
std::vector< int > fHounsfieldV
G4double GetMaxZ() const
G4double GetMaxX() const
G4double GetMinY() const
DcmDataset * theDataset
Definition: DicomVFile.hh:52
static int verbose
G4ThreeVector fOrientationRows
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
DicomFileMgr * theFileMgr
T max(const T t1, const T t2)
brief Return the largest of the two arguments
void DumpHeaderToTextFile(std::ofstream &fout)
G4ThreeVector fOrientationColumns
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
G4double GetMaxY() const
G4ThreeVector GetOrientationColumns() const
#define G4endl
Definition: G4ios.hh:61
void SetNoVoxelZ(const G4int &val)
G4int GetNoVoxelX() const
virtual void ReadData()
static G4GeometryTolerance * GetInstance()
static DicomFileMgr * GetInstance()
Definition: DicomFileMgr.cc:41
G4GLOB_DLL std::ostream G4cerr