Fraxinus  16.5.0-fx-rc1
An IGT application
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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) 2008-2014, SINTEF Department of Medical Technology
5 All rights reserved.
6 
7 Redistribution and use in source and binary forms, with or without
8 modification, are permitted provided that the following conditions are met:
9 
10 1. Redistributions of source code must retain the above copyright notice,
11  this list of conditions and the following disclaimer.
12 
13 2. Redistributions in binary form must reproduce the above copyright notice,
14  this list of conditions and the following disclaimer in the documentation
15  and/or other materials provided with the distribution.
16 
17 3. Neither the name of the copyright holder nor the names of its contributors
18  may be used to endorse or promote products derived from this software
19  without specific prior written permission.
20 
21 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
22 AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
23 IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
24 DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
25 FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
26 DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
27 SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
28 CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
29 OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
30 OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
31 =========================================================================*/
32 
33 #include "cxDicomImageReader.h"
34 
35 #include "cxVolumeHelpers.h"
36 #include "dcvrpn.h"
37 #include "cxLogger.h"
38 
39 namespace cx
40 {
41 
43 {
45  if (retval->loadFile(filename))
46  return retval;
47  else
48  return DicomImageReaderPtr();
49 }
50 
51 DicomImageReader::DicomImageReader() :
52  mDataset(NULL)
53 {
54 }
55 
56 bool DicomImageReader::loadFile(QString filename)
57 {
58  mFilename = filename;
59  OFCondition status = mFileFormat.loadFile(filename.toLatin1().data());
60  if( !status.good() )
61  {
62  return false;
63  }
64 
65  mDataset = mFileFormat.getDataset();
66  return true;
67 }
68 
70 {
71  return this->wrapInCTK(mDataset);
72 }
73 
74 double DicomImageReader::getDouble(const DcmTagKey& tag, const unsigned long pos, const OFBool searchIntoSub) const
75 {
76 // return this->getDouble(mDataset, tag, pos, searchIntoSub);
77  double retval = 0;
78  OFCondition condition;
79  condition = mDataset->findAndGetFloat64(tag, retval, pos, searchIntoSub);
80  if (!condition.good())
81  {
82  QString tagName = this->item()->TagDescription(tag);
83  this->error(QString("Failed to get tag %1/%2").arg(tagName).arg(pos));
84  }
85  return retval;
86 }
87 //double DicomImageReader::getDouble(DcmObject* dcmObject, const DcmTagKey& tag, const unsigned long pos, const OFBool searchIntoSub) const
88 //{
89 // DcmStack stack;
90 // dcmObject->search(tag, stack);
91 
92 // DcmElement* element = dynamic_cast<DcmElement*>(stack.top());
93 
94 // if(!element)
95 // {
96 // QString tagName = this->item()->TagDescription(tag);
97 // this->error(QString("Failed to get DcmAttributeTag with tag %1/%2").arg(tagName).arg(pos));
98 // return 0;
99 // }
100 // else
101 // {
102 // double val;
103 // element->getFloat64(val);
104 // }
105 
106 // double retval = 0;
107 // OFCondition condition;
108 // condition = element->getFloat64(retval, pos);
109 // if (!condition.good())
110 // {
111 // QString tagName = this->item()->TagDescription(tag);
112 // this->error(QString("Failed to get tag %1/%2").arg(tagName).arg(pos));
113 // }
114 // return retval;
115 //}
116 
118 {
119  WindowLevel retval;
120  retval.center = this->getDouble(DCM_WindowCenter, 0, OFTrue);
121  retval.width = this->getDouble(DCM_WindowWidth, 0, OFTrue);
122  return retval;
123 }
124 
126 {
127  //DICOM standard PS3.3 section C.7.6.1.1.2 Image Type
128  //http://dicom.nema.org/medical/dicom/current/output/html/part03.html#sect_C.7.6.1.1.2
129  bool retval = false;
130 
131  OFCondition condition;
132  OFString value;
133  condition = mDataset->findAndGetOFString(DCM_ImageType, value, 2, OFTrue);
134  if (condition.good())
135  {
136  QString imageSpesialization(value.c_str());
137  if (imageSpesialization.compare("LOCALIZER", Qt::CaseSensitive) == 0)
138  retval = true;
139  }
140  return retval;
141 }
142 
144 {
145  int numberOfFrames = this->item()->GetElementAsInteger(DCM_NumberOfFrames);
146  if (numberOfFrames==0)
147  {
148  unsigned short rows = 0;
149  unsigned short columns = 0;
150  mDataset->findAndGetUint16(DCM_Rows, rows, 0, OFTrue);
151  mDataset->findAndGetUint16(DCM_Columns, columns, 0, OFTrue);
152  if (rows*columns > 0)
153  numberOfFrames = 1; // seems like we have a 2D image
154  }
155  return numberOfFrames;
156 }
157 
159 {
160  Vector3D pos;
161  Vector3D e_x;
162  Vector3D e_y;
163 
164  for (int i=0; i<3; ++i)
165  {
166  e_x[i] = this->getDouble(DCM_ImageOrientationPatient, i, OFTrue);
167  e_y[i] = this->getDouble(DCM_ImageOrientationPatient, i+3, OFTrue);
168  pos[i] = this->getDouble(DCM_ImagePositionPatient, i, OFTrue);
169  }
170 
171  Vector3D zero_vec(0,0,0);
172  if( similar(e_x,zero_vec) && similar(e_y,zero_vec)) // Zero matrix
173  {
174  report("Set transform matrix to identity");
175  e_x[0]=1;
176  e_y[1]=1;
177  }
178 
179  Transform3D retval = cx::createTransformIJC(e_x, e_y, pos);
180  return retval;
181 }
182 ctkDICOMItemPtr DicomImageReader::wrapInCTK(DcmItem* item) const
183 {
184  if (!item)
185  return ctkDICOMItemPtr();
186  ctkDICOMItemPtr retval(new ctkDICOMItem);
187  retval->InitializeFromItem(item);
188  return retval;
189 }
190 
191 void DicomImageReader::error(QString message) const
192 {
193  reportError(QString("Dicom convert: [%1] in %2").arg(message).arg(mFilename));
194 }
195 
197 {
198  //TODO: Use DicomImage::createMonochromeImage() to get a monochrome copy for convenience
199 
200  DicomImage dicomImage(mFilename.toLatin1().data()); //, CIF_MayDetachPixelData );
201  const DiPixel *pixels = dicomImage.getInterData();
202  if (!pixels)
203  {
204  this->error("Found no pixel data");
205  return vtkImageDataPtr();
206  }
207 
208  vtkImageDataPtr data = vtkImageDataPtr::New();
209 
210  data->SetSpacing(this->getSpacing().data());
211 
212  Eigen::Array3i dim = this->getDim(dicomImage);
213  data->SetExtent(0, dim[0]-1, 0, dim[1]-1, 0, dim[2]-1);
214 
215  int samplesPerPixel = pixels->getPlanes();
216  int scalarSize = dim.prod() * samplesPerPixel;
217  int pixelDepth = dicomImage.getDepth();
218 
219  switch (pixels->getRepresentation())
220  {
221  case EPR_Uint8:
222 // std::cout << " VTK_UNSIGNED_CHAR" << std::endl;
223  data->AllocateScalars(VTK_UNSIGNED_CHAR, samplesPerPixel);
224  break;
225  case EPR_Uint16:
226 // std::cout << " VTK_UNSIGNED_SHORT" << std::endl;
227  data->AllocateScalars(VTK_UNSIGNED_SHORT, samplesPerPixel);
228  break;
229  case EPR_Uint32:
230  this->error("DICOM EPR_Uint32 not supported");
231  return vtkImageDataPtr();
232  break;
233  case EPR_Sint8:
234 // std::cout << " VTK_CHAR" << std::endl;
235  data->AllocateScalars(VTK_CHAR, samplesPerPixel);
236  break;
237  case EPR_Sint16:
238 // std::cout << " VTK_SHORT" << std::endl;
239  data->AllocateScalars(VTK_SHORT, samplesPerPixel);
240  break;
241  case EPR_Sint32:
242  this->error("DICOM EPR_Sint32 not supported");
243  return vtkImageDataPtr();
244  break;
245  }
246 
247  int bytesPerPixel = data->GetScalarSize() * samplesPerPixel;
248 
249  memcpy(data->GetScalarPointer(), pixels->getData(), pixels->getCount()*bytesPerPixel);
250  if (pixels->getCount()!=scalarSize)
251  this->error("Mismatch in pixel counts");
252  setDeepModified(data);
253  return data;
254 }
255 
256 Eigen::Array3d DicomImageReader::getSpacing() const
257 {
258  Eigen::Array3d spacing;
259  spacing[0] = this->getDouble(DCM_PixelSpacing, 0, OFTrue);
260  spacing[1] = this->getDouble(DCM_PixelSpacing, 1, OFTrue);
261 
262  double sliceThickness = this->getDouble(DCM_SliceThickness, 0, OFTrue);
263  spacing[2] = sliceThickness;
264 
265  if(this->isMultiFrameImage())
266  {
267  double sliceSpacing = this->getSliceSpacing();
268  if(similar(sliceSpacing, 0))
269  CX_LOG_WARNING() << "Cannot get slice spacing. Using slice thickness instead: " << sliceThickness;
270  else
271  spacing[2] = sliceSpacing;
272  }
273 
274 // double spacingBetweenSlices = this->getDouble(DCM_SpacingBetweenSlices, 0, OFTrue);//Usually only for MR
275 // std::cout << "DCM_SpacingBetweenSlices: " << spacingBetweenSlices << std::endl;
276 
277 // std::cout << " spacing: " << spacing << std::endl;
278  return spacing;
279 }
280 
281 bool DicomImageReader::isMultiFrameImage() const
282 {
283  //For now, just use number of z positions as indicator
284  QVector<double> zPos = this->getZPositions();
285  if(zPos.size() < 2)
286  return false;
287  return true;
288 }
289 
290 double DicomImageReader::getSliceSpacing() const
291 {
292  double retval;
293 
294  QVector<double> zPos = this->getZPositions();
295  if(zPos.size() < 2)
296  return 0;
297  retval = zPos[1] - zPos[0];
298 
299  for(int i = 2; i < zPos.size(); ++i)
300  {
301  double dist = zPos[i] - zPos[i-1];
302  if(!similar(dist, retval))
303  CX_LOG_WARNING() << "Distance between frame: " << i << " and " << i+1 << " is: " << dist << " != " << "dist between frame 0 and 1: " << retval;
304  }
305  if(retval < 0)
306  retval = zPos[0] - zPos[1];
307  return retval;
308 }
309 
310 QVector<double> DicomImageReader::getZPositions() const
311 {
312  QVector<double> retval;
313  DcmStack cleanStack;
314  DcmElement* stackElement;
315  OFCondition condition;
316  int i = 0;
317  do
318  {
319  do
320  condition = mDataset->nextObject(cleanStack, OFTrue);
321  while(condition.good() && cleanStack.top()->getTag() != DCM_ImagePositionPatient);
322 
323  ++i;
324  if(condition.good())
325  {
326  stackElement = dynamic_cast<DcmElement*>(cleanStack.top());
327  double val;
328  condition = stackElement->getFloat64(val, 2);
329  if(condition.bad())
330  {
331  CX_LOG_WARNING() << "Cannot get z pos for frame " << i;
332  return retval;
333  }
334  retval << val;
335 // CX_LOG_DEBUG() << "frame " << i << " z pos: " << val;
336  }
337  }
338  while(condition.good());
339  return retval;
340 }
341 
342 //Transform3D DicomImageReader::getTransform(DcmItem* dcmItem)
343 //{
344 
345 // e_x[i] = this->getDouble(DCM_ImageOrientationPatient, i, OFTrue);
346 // e_y[i] = this->getDouble(DCM_ImageOrientationPatient, i+3, OFTrue);
347 // pos[i] = this->getDouble(DCM_ImagePositionPatient, i, OFTrue);
348 
349 
350 // condition = stackElement->findAndGetOFString(DCM_ImagePositionPatient, value, 2, OFTrue);
351 //}
352 
353 Eigen::Array3i DicomImageReader::getDim(const DicomImage& dicomImage) const
354 {
355  Eigen::Array3i dim;
356  dim[0] = dicomImage.getWidth();
357  dim[1] = dicomImage.getHeight();
358  dim[2] = dicomImage.getFrameCount();
359  return dim;
360 }
361 
363 {
364  QString rawName = this->item()->GetElementAsString(DCM_PatientName);
365  return this->formatPatientName(rawName);
366 }
367 
368 QString DicomImageReader::formatPatientName(QString rawName) const
369 {
370  // ripped from ctkDICOMModel
371 
372  OFString dicomName = rawName.toStdString().c_str();
373  OFString formattedName;
374  OFString lastName, firstName, middleName, namePrefix, nameSuffix;
375  OFCondition l_error = DcmPersonName::getNameComponentsFromString(dicomName,
376  lastName, firstName, middleName, namePrefix, nameSuffix);
377  if (l_error.good())
378  {
379  formattedName.clear();
380  /* concatenate name components per this convention
381  * Last, First Middle, Suffix (Prefix)
382  * */
383  if (!lastName.empty())
384  {
385  formattedName += lastName;
386  if ( !(firstName.empty() && middleName.empty()) )
387  {
388  formattedName += ",";
389  }
390  }
391  if (!firstName.empty())
392  {
393  formattedName += " ";
394  formattedName += firstName;
395  }
396  if (!middleName.empty())
397  {
398  formattedName += " ";
399  formattedName += middleName;
400  }
401  if (!nameSuffix.empty())
402  {
403  formattedName += ", ";
404  formattedName += nameSuffix;
405  }
406  if (!namePrefix.empty())
407  {
408  formattedName += " (";
409  formattedName += namePrefix;
410  formattedName += ")";
411  }
412  }
413  return QString(formattedName.c_str());
414 }
415 
416 
417 } // namespace cx
418 
static DicomImageReaderPtr createFromFile(QString filename)
void reportError(QString msg)
Definition: cxLogger.cpp:92
Transform3D Transform3D
Transform3D is a representation of an affine 3D transform.
WindowLevel getWindowLevel() const
ctkDICOMItemPtr item() const
Transform3D getImageTransformPatient() const
vtkImageDataPtr createVtkImageData()
bool similar(const DoubleBoundingBox3D &a, const DoubleBoundingBox3D &b, double tol)
Transform3D createTransformIJC(const Vector3D &ivec, const Vector3D &jvec, const Vector3D &center)
QString getPatientName() const
boost::shared_ptr< class ctkDICOMItem > ctkDICOMItemPtr
boost::shared_ptr< class DicomImageReader > DicomImageReaderPtr
Eigen::Vector3d Vector3D
Vector3D is a representation of a point or vector in 3D.
Definition: cxVector3D.h:63
void report(QString msg)
Definition: cxLogger.cpp:90
void setDeepModified(vtkImageDataPtr image)
#define CX_LOG_WARNING
Definition: cxLogger.h:113
vtkSmartPointer< class vtkImageData > vtkImageDataPtr