CustusX  2023.01.05-dev+develop.0da12
An IGT application
cxDicomImageReader.cpp
Go to the documentation of this file.
1 /*=========================================================================
2 This file is part of CustusX, an Image Guided Therapy Application.
3 
4 Copyright (c) SINTEF Department of Medical Technology.
5 All rights reserved.
6 
7 CustusX is released under a BSD 3-Clause license.
8 
9 See Lisence.txt (https://github.com/SINTEFMedtek/CustusX/blob/master/License.txt) for details.
10 =========================================================================*/
11 
12 #include "cxDicomImageReader.h"
13 
14 #include "cxVolumeHelpers.h"
15 #include "dcvrpn.h"
16 #include "cxLogger.h"
17 
18 namespace cx
19 {
20 
22 {
23  DicomImageReaderPtr retval(new DicomImageReader);
24  if (retval->loadFile(filename))
25  return retval;
26  else
27  return DicomImageReaderPtr();
28 }
29 
30 DicomImageReader::DicomImageReader() :
31  mDataset(NULL)
32 {
33 }
34 
35 bool DicomImageReader::loadFile(QString filename)
36 {
37  mFilename = filename;
38  OFCondition status = mFileFormat.loadFile(filename.toLatin1().data());
39  if( !status.good() )
40  {
41  return false;
42  }
43 
44  mDataset = mFileFormat.getDataset();
45  return true;
46 }
47 
49 {
50  return this->wrapInCTK(mDataset);
51 }
52 
53 double DicomImageReader::getDouble(const DcmTagKey& tag, const unsigned long pos, const OFBool searchIntoSub) const
54 {
55 // return this->getDouble(mDataset, tag, pos, searchIntoSub);
56  double retval = 0;
57  OFCondition condition;
58  condition = mDataset->findAndGetFloat64(tag, retval, pos, searchIntoSub);
59  if (!condition.good())
60  {
61  QString tagName = this->item()->TagDescription(tag);
62  this->error(QString("Failed to get tag %1/%2").arg(tagName).arg(pos));
63  }
64  return retval;
65 }
66 //double DicomImageReader::getDouble(DcmObject* dcmObject, const DcmTagKey& tag, const unsigned long pos, const OFBool searchIntoSub) const
67 //{
68 // DcmStack stack;
69 // dcmObject->search(tag, stack);
70 
71 // DcmElement* element = dynamic_cast<DcmElement*>(stack.top());
72 
73 // if(!element)
74 // {
75 // QString tagName = this->item()->TagDescription(tag);
76 // this->error(QString("Failed to get DcmAttributeTag with tag %1/%2").arg(tagName).arg(pos));
77 // return 0;
78 // }
79 // else
80 // {
81 // double val;
82 // element->getFloat64(val);
83 // }
84 
85 // double retval = 0;
86 // OFCondition condition;
87 // condition = element->getFloat64(retval, pos);
88 // if (!condition.good())
89 // {
90 // QString tagName = this->item()->TagDescription(tag);
91 // this->error(QString("Failed to get tag %1/%2").arg(tagName).arg(pos));
92 // }
93 // return retval;
94 //}
95 
96 DicomImageReader::WindowLevel DicomImageReader::getWindowLevel() const
97 {
98  WindowLevel retval;
99  retval.center = this->getDouble(DCM_WindowCenter, 0, OFTrue);
100  retval.width = this->getDouble(DCM_WindowWidth, 0, OFTrue);
101  return retval;
102 }
103 
105 {
106  //DICOM standard PS3.3 section C.7.6.1.1.2 Image Type
107  //http://dicom.nema.org/medical/dicom/current/output/html/part03.html#sect_C.7.6.1.1.2
108  bool retval = false;
109 
110  OFCondition condition;
111  OFString value;
112  condition = mDataset->findAndGetOFString(DCM_ImageType, value, 2, OFTrue);
113  if (condition.good())
114  {
115  QString imageSpesialization(value.c_str());
116  if (imageSpesialization.compare("LOCALIZER", Qt::CaseSensitive) == 0)
117  retval = true;
118  }
119  return retval;
120 }
121 
123 {
124  int numberOfFrames = this->item()->GetElementAsInteger(DCM_NumberOfFrames);
125  if (numberOfFrames==0)
126  {
127  unsigned short rows = 0;
128  unsigned short columns = 0;
129  mDataset->findAndGetUint16(DCM_Rows, rows, 0, OFTrue);
130  mDataset->findAndGetUint16(DCM_Columns, columns, 0, OFTrue);
131  if (rows*columns > 0)
132  numberOfFrames = 1; // seems like we have a 2D image
133  }
134  return numberOfFrames;
135 }
136 
138 {
139  Vector3D pos;
140  Vector3D e_x;
141  Vector3D e_y;
142 
143  for (int i=0; i<3; ++i)
144  {
145  e_x[i] = this->getDouble(DCM_ImageOrientationPatient, i, OFTrue);
146  e_y[i] = this->getDouble(DCM_ImageOrientationPatient, i+3, OFTrue);
147  pos[i] = this->getDouble(DCM_ImagePositionPatient, i, OFTrue);
148  }
149 
150  Vector3D zero_vec(0,0,0);
151  if( similar(e_x,zero_vec) && similar(e_y,zero_vec)) // Zero matrix
152  {
153  report("Set transform matrix to identity");
154  e_x[0]=1;
155  e_y[1]=1;
156  }
157 
158  Transform3D retval = cx::createTransformIJC(e_x, e_y, pos);
159  return retval;
160 }
161 ctkDICOMItemPtr DicomImageReader::wrapInCTK(DcmItem* item) const
162 {
163  if (!item)
164  return ctkDICOMItemPtr();
165  ctkDICOMItemPtr retval(new ctkDICOMItem);
166  retval->InitializeFromItem(item);
167  return retval;
168 }
169 
170 void DicomImageReader::error(QString message) const
171 {
172  reportError(QString("Dicom convert: [%1] in %2").arg(message).arg(mFilename));
173 }
174 
176 {
177  //TODO: Use DicomImage::createMonochromeImage() to get a monochrome copy for convenience
178 
179  DicomImage dicomImage(mFilename.toLatin1().data()); //, CIF_MayDetachPixelData );
180  const DiPixel *pixels = dicomImage.getInterData();
181  if (!pixels)
182  {
183  this->error("Found no pixel data");
184  return vtkImageDataPtr();
185  }
186 
187  vtkImageDataPtr data = vtkImageDataPtr::New();
188 
189  data->SetSpacing(this->getSpacing().data());
190 
191  Eigen::Array3i dim = this->getDim(dicomImage);
192  data->SetExtent(0, dim[0]-1, 0, dim[1]-1, 0, dim[2]-1);
193 
194  int samplesPerPixel = pixels->getPlanes();
195  int scalarSize = dim.prod() * samplesPerPixel;
196  int pixelDepth = dicomImage.getDepth();
197 
198  switch (pixels->getRepresentation())
199  {
200  case EPR_Uint8:
201 // std::cout << " VTK_UNSIGNED_CHAR" << std::endl;
202  data->AllocateScalars(VTK_UNSIGNED_CHAR, samplesPerPixel);
203  break;
204  case EPR_Uint16:
205 // std::cout << " VTK_UNSIGNED_SHORT" << std::endl;
206  data->AllocateScalars(VTK_UNSIGNED_SHORT, samplesPerPixel);
207  break;
208  case EPR_Uint32:
209 // std::cout << " VTK_UNSIGNED_INT" << std::endl;
210  data->AllocateScalars(VTK_UNSIGNED_INT, samplesPerPixel);
211  break;
212  case EPR_Sint8:
213 // std::cout << " VTK_CHAR" << std::endl;
214  data->AllocateScalars(VTK_CHAR, samplesPerPixel);
215  break;
216  case EPR_Sint16:
217 // std::cout << " VTK_SHORT" << std::endl;
218  data->AllocateScalars(VTK_SHORT, samplesPerPixel);
219  break;
220  case EPR_Sint32:
221 // std::cout << " VTK_INT" << std::endl;
222  data->AllocateScalars(VTK_INT, samplesPerPixel);
223  break;
224  }
225 
226  int bytesPerPixel = data->GetScalarSize() * samplesPerPixel;
227 
228  memcpy(data->GetScalarPointer(), pixels->getData(), pixels->getCount()*bytesPerPixel);
229  if (pixels->getCount()!=scalarSize)
230  this->error("Mismatch in pixel counts");
231  setDeepModified(data);
232  return data;
233 }
234 
235 Eigen::Array3d DicomImageReader::getSpacing() const
236 {
237  Eigen::Array3d spacing;
238  spacing[0] = this->getDouble(DCM_PixelSpacing, 0, OFTrue);
239  spacing[1] = this->getDouble(DCM_PixelSpacing, 1, OFTrue);
240  spacing[2] = 0;
241 
242  if(this->isMultiFrameImage())
243  {
244  double sliceSpacing = this->getSliceSpacing();
245  if(similar(sliceSpacing, 0))
246  CX_LOG_WARNING() << "Cannot get slice spacing. Spacing set to 0. ";
247  else
248  spacing[2] = sliceSpacing;
249  }
250 
251  return spacing;
252 }
253 
254 bool DicomImageReader::isMultiFrameImage() const
255 {
256  //For now, just use number of z positions as indicator
257  QVector<double> zPos = this->getPositions(2);
258  if(zPos.size() < 2)
259  return false;
260  return true;
261 }
262 
263 double DicomImageReader::calculateMultiFrameSpacing(int frameIndex) const
264 {
265  static QVector<double> xPos = this->getPositions(0);
266  static QVector<double> yPos = this->getPositions(1);
267  static QVector<double> zPos = this->getPositions(2);
268 
269  if (frameIndex < 1 || frameIndex > (zPos.length()-1))
270  return 0;
271 
272  double dist = std::sqrt( pow(xPos[frameIndex] - xPos[frameIndex-1],2)
273  + pow(yPos[frameIndex] - yPos[frameIndex-1],2) + pow(zPos[frameIndex] - zPos[frameIndex-1],2) );
274  return dist;
275 }
276 
277 double DicomImageReader::getSliceSpacing() const
278 {
279  QVector<double> zPos = this->getPositions(2);
280  if(zPos.size() < 2)
281  return 0;
282 
283  double retval = calculateMultiFrameSpacing(1);
284 
285  double tolerance = retval/10000;
286 
287  for(int i = 2; i < zPos.size(); ++i)
288  {
289  double dist = calculateMultiFrameSpacing(i);
290  if(!similar(dist, retval,tolerance))
291  CX_LOG_WARNING() << "Distance between frame: " << i - 1 << " and " << i << " is: "
292  << dist << " != " << "Dist between frame 0 and 1: " << retval;
293  }
294  return retval;
295 }
296 
297 QVector<double> DicomImageReader::getPositions(int cIndex) const
298 {
299  QVector<double> retval;
300  DcmStack cleanStack;
301  DcmElement* stackElement;
302  OFCondition condition;
303  int i = 0;
304  do
305  {
306  do
307  condition = mDataset->nextObject(cleanStack, OFTrue);
308  while(condition.good() && cleanStack.top()->getTag() != DCM_ImagePositionPatient);
309 
310  ++i;
311  if(condition.good())
312  {
313  stackElement = dynamic_cast<DcmElement*>(cleanStack.top());
314  double val;
315  condition = stackElement->getFloat64(val, cIndex);
316  if(condition.bad())
317  {
318  CX_LOG_WARNING() << "Cannot get pos for frame " << i;
319  return retval;
320  }
321  retval << val;
322  // CX_LOG_DEBUG() << "frame " << i << " pos: " << val;
323  }
324  }
325  while(condition.good());
326  return retval;
327 }
328 
329 //Transform3D DicomImageReader::getTransform(DcmItem* dcmItem)
330 //{
331 
332 // e_x[i] = this->getDouble(DCM_ImageOrientationPatient, i, OFTrue);
333 // e_y[i] = this->getDouble(DCM_ImageOrientationPatient, i+3, OFTrue);
334 // pos[i] = this->getDouble(DCM_ImagePositionPatient, i, OFTrue);
335 
336 
337 // condition = stackElement->findAndGetOFString(DCM_ImagePositionPatient, value, 2, OFTrue);
338 //}
339 
340 Eigen::Array3i DicomImageReader::getDim(const DicomImage& dicomImage) const
341 {
342  Eigen::Array3i dim;
343  dim[0] = dicomImage.getWidth();
344  dim[1] = dicomImage.getHeight();
345  dim[2] = dicomImage.getFrameCount();
346  return dim;
347 }
348 
349 QString DicomImageReader::getPatientName() const
350 {
351  QString rawName = this->item()->GetElementAsString(DCM_PatientName);
352  return this->formatPatientName(rawName);
353 }
354 
355 QString DicomImageReader::formatPatientName(QString rawName) const
356 {
357  // ripped from ctkDICOMModel
358 
359  OFString dicomName = rawName.toStdString().c_str();
360  OFString formattedName;
361  OFString lastName, firstName, middleName, namePrefix, nameSuffix;
362  OFCondition l_error = DcmPersonName::getNameComponentsFromString(dicomName,
363  lastName, firstName, middleName, namePrefix, nameSuffix);
364  if (l_error.good())
365  {
366  formattedName.clear();
367  /* concatenate name components per this convention
368  * Last, First Middle, Suffix (Prefix)
369  * */
370  if (!lastName.empty())
371  {
372  formattedName += lastName;
373  if ( !(firstName.empty() && middleName.empty()) )
374  {
375  formattedName += ",";
376  }
377  }
378  if (!firstName.empty())
379  {
380  formattedName += " ";
381  formattedName += firstName;
382  }
383  if (!middleName.empty())
384  {
385  formattedName += " ";
386  formattedName += middleName;
387  }
388  if (!nameSuffix.empty())
389  {
390  formattedName += ", ";
391  formattedName += nameSuffix;
392  }
393  if (!namePrefix.empty())
394  {
395  formattedName += " (";
396  formattedName += namePrefix;
397  formattedName += ")";
398  }
399  }
400  return QString(formattedName.c_str());
401 }
402 
403 
404 } // namespace cx
405 
static DicomImageReaderPtr createFromFile(QString filename)
void reportError(QString msg)
Definition: cxLogger.cpp:71
Transform3D Transform3D
Transform3D is a representation of an affine 3D transform.
WindowLevel getWindowLevel() const
ctkDICOMItemPtr item() const
Transform3D getImageTransformPatient() const
vtkImageDataPtr createVtkImageData()
boost::shared_ptr< class ctkDICOMItem > ctkDICOMItemPtr
Transform3D createTransformIJC(const Vector3D &ivec, const Vector3D &jvec, const Vector3D &center)
QString getPatientName() const
boost::shared_ptr< class DicomImageReader > DicomImageReaderPtr
Eigen::Vector3d Vector3D
Vector3D is a representation of a point or vector in 3D.
Definition: cxVector3D.h:42
void report(QString msg)
Definition: cxLogger.cpp:69
void setDeepModified(vtkImageDataPtr image)
bool similar(const CameraInfo &lhs, const CameraInfo &rhs, double tol)
#define CX_LOG_WARNING
Definition: cxLogger.h:98
vtkSmartPointer< class vtkImageData > vtkImageDataPtr
Namespace for all CustusX production code.