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 {
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 
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  default:
225  CX_LOG_WARNING() << "Unknown pixel format: " << pixels->getRepresentation();
226  return vtkImageDataPtr();
227  break;
228  }
229 
230  int bytesPerPixel = data->GetScalarSize() * samplesPerPixel;
231 
232  if (pixels->getCount()!=scalarSize)
233  {
234  CX_LOG_WARNING() << "DicomImageReader::createVtkImageData: Mismatch in pixel counts: " << pixels->getCount() << " != " << scalarSize << " Wrong number of bytesPerPixel? = " << bytesPerPixel;
235  return vtkImageDataPtr();//Prevent seg.fault. in memcopy line below
236  }
237  memcpy(data->GetScalarPointer(), pixels->getData(), pixels->getCount()*bytesPerPixel);
238  setDeepModified(data);
239  return data;
240 }
241 
242 Eigen::Array3d DicomImageReader::getSpacing() const
243 {
244  Eigen::Array3d spacing;
245  spacing[0] = this->getDouble(DCM_PixelSpacing, 0, OFTrue);
246  spacing[1] = this->getDouble(DCM_PixelSpacing, 1, OFTrue);
247  spacing[2] = 0;
248 
249  if(this->isMultiFrameImage())
250  {
251  double sliceSpacing = this->getSliceSpacing();
252  if(similar(sliceSpacing, 0))
253  CX_LOG_WARNING() << "Cannot get slice spacing. Spacing set to 0. ";
254  else
255  spacing[2] = sliceSpacing;
256  }
257 
258  return spacing;
259 }
260 
261 bool DicomImageReader::isMultiFrameImage() const
262 {
263  //For now, just use number of z positions as indicator
264  QVector<double> zPos = this->getPositions(2);
265  if(zPos.size() < 2)
266  return false;
267  return true;
268 }
269 
270 double DicomImageReader::calculateMultiFrameSpacing(int frameIndex) const
271 {
272  static QVector<double> xPos = this->getPositions(0);
273  static QVector<double> yPos = this->getPositions(1);
274  static QVector<double> zPos = this->getPositions(2);
275 
276  if (frameIndex < 1 || frameIndex > (zPos.length()-1))
277  return 0;
278 
279  double dist = std::sqrt( pow(xPos[frameIndex] - xPos[frameIndex-1],2)
280  + pow(yPos[frameIndex] - yPos[frameIndex-1],2) + pow(zPos[frameIndex] - zPos[frameIndex-1],2) );
281  return dist;
282 }
283 
284 double DicomImageReader::getSliceSpacing() const
285 {
286  QVector<double> zPos = this->getPositions(2);
287  if(zPos.size() < 2)
288  return 0;
289 
290  double retval = calculateMultiFrameSpacing(1);
291 
292  double tolerance = retval/10000;
293 
294  for(int i = 2; i < zPos.size(); ++i)
295  {
296  double dist = calculateMultiFrameSpacing(i);
297  if(!similar(dist, retval,tolerance))
298  CX_LOG_WARNING() << "Distance between frame: " << i - 1 << " and " << i << " is: "
299  << dist << " != " << "Dist between frame 0 and 1: " << retval;
300  }
301  return retval;
302 }
303 
304 QVector<double> DicomImageReader::getPositions(int cIndex) const
305 {
306  QVector<double> retval;
307  DcmStack cleanStack;
308  DcmElement* stackElement;
309  OFCondition condition;
310  int i = 0;
311  do
312  {
313  do
314  condition = mDataset->nextObject(cleanStack, OFTrue);
315  while(condition.good() && cleanStack.top()->getTag() != DCM_ImagePositionPatient);
316 
317  ++i;
318  if(condition.good())
319  {
320  stackElement = dynamic_cast<DcmElement*>(cleanStack.top());
321  double val;
322  condition = stackElement->getFloat64(val, cIndex);
323  if(condition.bad())
324  {
325  CX_LOG_WARNING() << "Cannot get pos for frame " << i;
326  return retval;
327  }
328  retval << val;
329  // CX_LOG_DEBUG() << "frame " << i << " pos: " << val;
330  }
331  }
332  while(condition.good());
333  return retval;
334 }
335 
336 //Transform3D DicomImageReader::getTransform(DcmItem* dcmItem)
337 //{
338 
339 // e_x[i] = this->getDouble(DCM_ImageOrientationPatient, i, OFTrue);
340 // e_y[i] = this->getDouble(DCM_ImageOrientationPatient, i+3, OFTrue);
341 // pos[i] = this->getDouble(DCM_ImagePositionPatient, i, OFTrue);
342 
343 
344 // condition = stackElement->findAndGetOFString(DCM_ImagePositionPatient, value, 2, OFTrue);
345 //}
346 
347 Eigen::Array3i DicomImageReader::getDim(const DicomImage& dicomImage) const
348 {
349  Eigen::Array3i dim;
350  dim[0] = dicomImage.getWidth();
351  dim[1] = dicomImage.getHeight();
352  dim[2] = dicomImage.getFrameCount();
353  return dim;
354 }
355 
357 {
358  QString rawName = this->item()->GetElementAsString(DCM_PatientName);
359  return this->formatPatientName(rawName);
360 }
361 
362 QString DicomImageReader::formatPatientName(QString rawName) const
363 {
364  // ripped from ctkDICOMModel
365 
366  OFString dicomName = rawName.toStdString().c_str();
367  OFString formattedName;
368  OFString lastName, firstName, middleName, namePrefix, nameSuffix;
369  OFCondition l_error = DcmPersonName::getNameComponentsFromString(dicomName,
370  lastName, firstName, middleName, namePrefix, nameSuffix);
371  if (l_error.good())
372  {
373  formattedName.clear();
374  /* concatenate name components per this convention
375  * Last, First Middle, Suffix (Prefix)
376  * */
377  if (!lastName.empty())
378  {
379  formattedName += lastName;
380  if ( !(firstName.empty() && middleName.empty()) )
381  {
382  formattedName += ",";
383  }
384  }
385  if (!firstName.empty())
386  {
387  formattedName += " ";
388  formattedName += firstName;
389  }
390  if (!middleName.empty())
391  {
392  formattedName += " ";
393  formattedName += middleName;
394  }
395  if (!nameSuffix.empty())
396  {
397  formattedName += ", ";
398  formattedName += nameSuffix;
399  }
400  if (!namePrefix.empty())
401  {
402  formattedName += " (";
403  formattedName += namePrefix;
404  formattedName += ")";
405  }
406  }
407  return QString(formattedName.c_str());
408 }
409 
410 
411 } // namespace cx
412 
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.